*========================================= PROGRAM STATICS *----------------------------------------- * Yuri Mishin, Feb.-March 1996; revised June 1998 *==================================================================== * * Calculation of equilibrium structure and energy of point and * planar defects by "molecular statics" using the embedded atom * method (EAM). Relaxation with respect to atomic displacements * of free atoms is performed by the conjugate gradient method. * For a planar defect in the X-Z plane, rigid-body displacements * of the grains in the Y-directon can be included as an additional * mode of relaxation. The block can be additionally relaxed with * respect to its uniform isotropic deformation (expansion or * contraction) or elastic deformation in Y-direction * * This version of the program includes simulation of a mode I * crack in the block * * Program limitations: * maximum number of chemical species: 4 * maximum number of sublattices: 8 * maximum total number of atoms: 100000 * maximum number of neighbours 150 * *IMPORTANT: For suggestions, comments or bugs contact mishin@vt.edu * I am not responsible for any unauthorized modifications * of the original version of this code. - Y.M. *=================================================================== implicit real*8 (a-h,o-z) real*8 KI, mu, nu logical bcx,bcy,bcz,ipoint,iplane,dfrm character*72 outfile,plotfile,strfile,deffile,crcfile character*1 dummy parameter (nbasmax = 30000, ! max number of basis atoms & maxatom = 100000, ! max number of atoms in block & maxw = 15*maxatom + 2 ! max dimension of array w & , nspmax=4) common /accuracy/ acc common /basis/ xbas(nbasmax),ybas(nbasmax),zbas(nbasmax), & lbas(nbasmax),nbas common /block/ search2,xmn,ymn,zmn,xmx,ymx,zmx, & xtmn,ytmn,ztmn,xtmx,ytmx,ztmx, & xdisp,ydisp,zdisp,mdx,mdy,mdz, & ntot,nb,nfree common /atoms/ lco(maxatom) character*4 element(nspmax) common/chemical/cutoff2,np(nspmax,nspmax),numat,npairs,npr, & element common /cyc/ bcx,bcy,bcz common /defct/ nvac,ilat,ipoint,iplane common /crack/ x0,y0,KI,mu,nu common/minimize/ xdisp1,ydisp1,zdisp1,minmode,iter1D dimension co(3*maxatom),co1(3*maxatom),g(3*maxatom) dimension w(maxw) data ipoint/.false./, iplane/.false./, acc /1.D-4/ open (1, file = 'crystal.dat', status= 'old') read (1,*) outfile open (6, file = outfile, status= 'unknown' ) write (6, 300) write (6, 301) ****************************************************************** * Develop interpolation coefficients for interatomic potentials * ****************************************************************** call INTERPOL ********************************************* * Create the initial atomic configuration * ********************************************* read (1,*) ilat read (1,*) strfile open (5, file = strfile, status= 'old') * --- if ilat.ne.0, the intial structure is generated using * --- subroutine CRYSTAL if (ilat.ne.0) then write(6, 130) strfile call CRYSTAL(co) ! Generate a perfect lattice block call PERFECT(co,E0) ! and calculate its energy else * --- Otherwise the intial structure is read from a file write(6, 131) strfile read(5,403) dummy, xmn, ymn, zmn read(5,403) dummy, xmx, ymx, zmx read(5,403) dummy, xtmn, ytmn, ztmn read(5,403) dummy, xtmx, ytmx, ztmx read(5,404) dummy, nbas, ntot, nb, nfree if(ntot.gt.maxatom) then write(6,135) maxatom stop endif read(5,405) dummy, search, mdx, mdy, mdz read(5,404) dummy, ibcx, ibcy, ibcz read(5,404) dummy, ipo, ipl read(5,403) dummy, E0 xdisp = xtmx - xtmn ydisp = ytmx - ytmn zdisp = ztmx - ztmn search2 = search**2 if(ibcx.lt.0) then ! decode the boundary conditions bcx = .true. else bcx = .false. endif if(ibcy.lt.0) then bcy = .true. else bcy = .false. endif if(ibcz.lt.0) then bcz = .true. else bcz = .false. endif if(ipo.gt.0) then ipoint = .true. else ipoint = .false. endif if(ipl.gt.0) then iplane = .true. else iplane = .false. endif do i = 1,ntot k1 = i*3 j1 = k1 -1 i1 = j1 -1 read(5,127) idum, co(i1), co(j1), co(k1), lco(i) enddo endif close(5) * Introduce defects if required read (1,*) idef if (idef.ne.0) then read (1,*) deffile open (5, file = deffile, status= 'old') write(6, 132) deffile call DEFECT(co) endif ************************************* * Construct tables of neighbours * ************************************* call TABLES(co,.true.) ************************************************************* * Calculate the energy of the initial (defected) block * ************************************************************* if (ipoint.or.iplane.or.ilat.eq.0) then write(6, 100) ENERGY(co) endif ************************************ * Deform the block if requested * ************************************ ntot3 = 3*ntot do i = 1, ntot3 ! save the initial structure co1(i) = co(i) enddo xdisp1 = xdisp ydisp1 = ydisp zdisp1 = zdisp read(1,*) xs1,ys1,zs1,rbd1 if(xs1.ne.1.D0.or.ys1.ne.1.D0.or.zs1.ne.1.D0.or.rbd1.ne.0.D0) & then dfrm = .true. CALL DEFORM(xs1,ys1,zs1,rbd1,co,co1) write(6,133) xs1,ys1,zs1,rbd1 do i = 1, ntot3 ! save the structure after deformation co1(i) = co(i) enddo xdisp1 = xdisp ydisp1 = ydisp zdisp1 = zdisp else dfrm = .false. endif *========= Begin the minimization process ============ * Read in the minimization mode * * minmode = 0 No energy minimization. Proceed to label 77 * minmode = 1 One energy minimization with fixed sizes of the block * minmode = 2 Energy minimization with respect to both local atomic * displacements and uniform isotropic deformation * (expansion or contraction) of the block * minmode = 3 Energy minimization with respect to both local atomic * displacements and elastic deformation of the block * (expansion or contraction) in Y direction * minmode = 4 Energy minimization with respect to both local atomic * displacements and rigid-body translations of buffer * and fixed atoms in Y direction * minmode = 5 Introduce a mode I crack and minimize the energy for * several values of the stress intensity factor KI * read (1,*) minmode read (1,*) mute read (1,*) iplot if (iplot.eq.1) read (1,*) plotfile * minmode = 0 if (minmode.eq.0) then write(6,134) go to 77 endif tol = 1.D0 1 tol = tol/2 if (1.D0 + tol .gt. 1.D0) go to 1 tol = 2*tol eps = 1.D-5 tol1D = 1.D-6 epsE =1.D-5 mxfun = 500 idev = 6 if (mute.eq.0) then iout = 0 else iout = 1 endif nfree3 = 3*nfree Deform0 = 1.D0 rbd = 0.D0 * minmode = 1 if(minmode.eq.1) then write(6,113) call CONJUG(nfree3,co,E,g,ifun,iter,eps,nflag,mxfun,w, & iout,idev,tol) endif * minmode = 2 or 3 if(minmode.eq.2.or.minmode.eq.3) then if (minmode.eq.2) write(6,119) if (minmode.eq.3) write(6,117) E1 = 0.D0 iternum = 0 2 iternum = iternum + 1 Defor = FMIN(0.95D0,1.05D0,tol1D,co,co1) if (minmode.eq.2) CALL DEFORM(Defor,Defor,Defor,0.D0,co,co1) if (minmode.eq.3) CALL DEFORM(1.D0,Defor,1.D0,0.D0,co,co1) call CONJUG(nfree3,co,E,g,ifun,iter,eps,nflag,mxfun,w, & iout,idev,tol) Deform0 = Deform0*Defor if (DABS(E1-E).gt.epsE) then E1 = E write(6, 120) iter1D, Deform0, E iter1D = 0 do i = 1, ntot3 co1(i) = co(i) enddo ydisp1 = ydisp if(minmode.eq.2) then xdisp1 = xdisp zdisp1 = zdisp endif go to 2 endif write(6, 121) iternum, Deform0 endif * minmode = 4 if(minmode.eq.4) then write(6,114) E1 = 0.D0 iternum = 0 3 iternum = iternum + 1 rigid = FMIN(-0.1D0,0.1D0,tol1D,co,co1) CALL DEFORM(1.D0,1.D0,1.D0,rigid,co,co1) call CONJUG(nfree3,co,E,g,ifun,iter,eps,nflag,mxfun,w, & iout,idev,tol) rbd = rbd + rigid if (DABS(E1-E).gt.epsE) then E1 = E write(6, 116) iter1D, rbd, E iter1D = 0 do i = 1, ntot3 co1(i) = co(i) enddo ydisp1 = ydisp go to 3 endif write(6, 115) iternum, rbd endif ys = ys1*Deform0 xs = xs1 zs = zs1 if(minmode.eq.2) then xs = xs*Deform0 zs = zs*Deform0 endif * minmode = 5 if(minmode.eq.5) then read (1,*) crcfile open (5, file = crcfile, status= 'old') write(6, 136) crcfile read(5,*) x0,y0 ! initial position of crack tip read(5,*) mu ! shear modulus read(5,*) nu ! Poisson's ratio read(5,*) KI, dKI, nk ! initial stress intensity factor, ! its increment and number of steps read(5,*) r1, nmin ! parameters for crack tip location close(5) istep = 1 fact = KI xnew = x0 ynew = y0 write(6,138) 5 CALL MODE1A(co,co1) ! open the crack call CONJUG(nfree3,co,E,g,ifun,iter,eps,nflag,mxfun,w, & iout,idev,tol) ! relax the structure E = ENERGY(co) * Find the new location of the crack tip CALL LOCATE(r1,nmin,co,xnew,ynew) write(6,137) istep, fact, E-nb*E0, xnew, ynew, nflag fact = fact + dKI KI = dKI ! increase the load if(istep.lt.nk) then istep = istep + 1 do i = 1, ntot3 ! save the relaxed structure co1(i) = co(i) ! and repeat the cycle enddo go to 5 endif endif * ========= End of minimization ============= if (nflag.eq.0) write (6, 101) if (nflag.eq.1) write (6, 102) if (nflag.eq.2) write (6, 103) if (nflag.eq.3) write (6, 104) write(6, 105) iter, ifun gmax = 0.D0 dmax = 0.D0 do i = 1, nfree k1 = 3*i j1 = k1 -1 i1 = j1 -1 dr = (co(i1)-co1(i1))**2 +(co(j1)-co1(j1))**2 + & (co(k1)-co1(k1))**2 dg = g(i1)**2 + g(j1)**2 + g(k1)**2 if (dr.gt.dmax) dmax = dr if (dg.gt.gmax) gmax = dg enddo dmax = DSQRT(dmax) gmax = DSQRT(gmax) write(6, 109) dmax, gmax write(6, 106) E if (iplane.and..not.ipoint) then S = xdisp*zdisp gamma = (E-E0*nb)/S write(6, 107) gamma else if (ipoint) then write(6, 108) E-E0*nb endif 77 continue ******************************************************** * Save the final configuration in a file if required * ******************************************************** if (iplot.eq.1) then search = DSQRT(search2) * encode the boundary conditions if(bcx) then ibcx = -1 else ibcx = 1 endif if(bcy) then ibcy = -1 else ibcy = 1 endif if(bcz) then ibcz = -1 else ibcz = 1 endif if(ipoint) then ipo = 1 else ipo = 0 endif if(iplane) then ipl = 1 else ipl = 0 endif * update the block sizes if (dfrm.or.minmode.gt.1) then xmn = xs*xmn ymn = ys*ymn - rbd1*Deform0 - rbd zmn = zs*zmn xmx = xs*xmx ymx = ys*ymx + rbd1*Deform0 + rbd zmx = zs*zmx xtmn = xs*xtmn ytmn = ys*ytmn - rbd1*Deform0 - rbd ztmn = zs*ztmn xtmx = xs*xtmx ytmx = ys*ytmx + rbd1*Deform0 + rbd ztmx = zs*ztmx endif open (7, file = plotfile, status= 'unknown' ) c write(7,400) xmn, ymn, zmn c write(7,400) xmx, ymx, zmx c write(7,400) xtmn, ytmn, ztmn c write(7,400) xtmx, ytmx, ztmx c write(7,401) nbas, ntot, nb, nfree c write(7,402) search, mdx, mdy, mdz c write(7,401) ibcx, ibcy, ibcz c write(7,401) ipo, ipl c write(7,*) ntot c modified by R.Kriz: change ntot to nfreem1 (assumes only one vacancy) nfreem1=nfree-1 write(7,*) nfreem1 write(7,400) E0 c do i =1,ntot c modified by R.Kriz: change ntot to nfreem1 do i =1,nfree k1 = i*3 j1 = k1 -1 i1 = j1 -1 c write(7,127) i, co(i1), co(j1), co(k1), lco(i) c write (7,128) element(lco(i)), co(i1), co(j1), co(k1) c modified by R.Kriz: test if vacancy exists and don't write if(lco(i).ne.0) write (7,128) element(lco(i)), & co(i1), co(j1), co(k1) enddo close(7) write(6, 125) plotfile endif write (6, 126) close(5) close(6) * ----------- FORMAT statements ------------------------------ 100 format(' INITIAL ENERGY OF THE DEFECTED BLOCK:',D15.7) 101 format(/' *** THE MINIMIZATION ALGORITHM SAFELY CONVERGED ***', & /' Congratulations !!! ') 102 format(/' THE MINIMIZATION PROCESS TERMINATED', &/ ' BECAUSE THE MAX NUMBER OF FUNCTION CALLS ATTAINED', &/' Check the maximum residual force below:', &/' If it is roughly < 0.01, the result is to be trusted.', &/' If it is roughly > 0.1, the relaxation is incomplete.') 103 format(' THE MINIMIZATION PROCESS TERMINATED', &/ ' BECAUSE OF COMPUTATIONAL PROBLEMS', &/' The linear search has failed to improve the energy value.', &/' The result can be quite reasonable though. Check the ', &/' maximum residual force below:', &/' If it is roughly < 0.01, the result is to be trusted.', &/' If it is roughly > 0.1, the relaxation is incomplete.') 104 format(' THE MINIMIZATION PROCESS TERMINATED', &/ ' BECAUSE OF COMPUTATIONAL PROBLEMS', &/' The search vector was not a descent direction. This can', &/' only be caused by roundoff, and may suggest that the ' &/' convergence criterion is too strict.', &/' Check the maximum residual force below:', &/' If it is roughly < 0.01, the result is to be trusted.', &/' If it is roughly > 0.1, the relaxation is incomplete.') 105 format(/' TOTAL NUMBER OF INTERATONS MADE:', i4, &/ ' TOTAL NUMBER OF FUNCTION CALLS:', I5/) 106 format(' ----------------------------------------------------', &/ ' MINIMUM ENERGY OF THE BLOCK:',D15.7) 107 format(' SPECIFIC ENERGY OF THE PLANAR DEFECT:',D15.7, &/ ' ----------------------------------------------------'/) 108 format(' DEFECT ENERGY:',D15.7, &/ ' ----------------------------------------------------'/) 109 format(' MAXIMUM ATOMIC DISPLACEMENT DUE TO RELAXATION:',f8.5, &/' MAXIMUM FORCE AFTER RELAXATION:', f9.6/) 110 format(' WRONG MINIMIZATION MODE (minmode =',i2, &' IN FUNCTION F)') 112 format(/' START THE MINIMIZATION PROCESS... '/) 113 format(/' MINIMIZATION IN MODE 1:', &/ ' ATOMIC RELAXATION IN THE FREE BLOCK. NO VARIATION OF', &/ ' THE BLOCK DIMENSIONS. NO RIGID-BODY DISPLACEMENTS '/) 114 format(/' MINIMIZATION IN MODE 4:', &/ ' ATOMIC RELAXATION OF THE FREE BLOCK ALONG WITH', &/ ' RIGID-BODY DISPLACEMENTS IN Y-DIRECTION'/) 115 format(/i15,' ITERATIONS MADE FOR RIGID-BODY DISPLACEMENTS', &/ ' EQUILIBRIUM RIGID-BODY DISPLACEMENT OF EACH GRAIN:', &f9.6,// ' RESULTS FOR THE LAST ATOMIC RELAXATION:'/) 116 format(' NO. OF FUNC. CALLS =', i3,' R-B DISPL. =', f9.6, &' ENERGY =',D15.7) 117 format(/' MINIMIZATION IN MODE 3:', &/ ' ATOMIC RELAXATION OF THE FREE BLOCK ALONG WITH', &/ ' ELASTIC DEFORMATION in Y-DIRECTION'/) 119 format(/' MINIMIZATION IN MODE 2:', &/ ' ATOMIC RELAXATION OF THE FREE BLOCK ALONG WITH', &/ ' UNIFORM ISOTROPIC DEFORMATION'/) 120 format(' NO. OF FUNC. CALLS =', i3,' DEFORM. =', f9.6, &' ENERGY =',D15.7) 121 format(/i6,' ITERATIONS MADE FOR DEFORMATION OF THE BLOCK', & / ' EQUILIBRIUM DEFORMATION FOUND:', &f9.6,// ' RESULTS FOR THE LAST ATOMIC RELAXATION:') 122 format( ' MINIMIZATION IN MODE 2:', &/ ' ATOMIC RELAXATION OF THE FREE BLOCK WITH', &/ ' A FIXED UNIFORM DEFORMATION BY',f7.4) 124 format(i5, 6f10.5, i4) 125 format(/' THE FINAL CONFIGURATION HAS BEEN SAVED IN FILE:', &/18x,a72/) 126 format(40x,'... END OF OUTPUT.'/) 127 format(I5,3E19.8, i3) 128 format (A4,3E14.3) 130 format(/' INITIAL CONFIGURATION IS GENERATED BY', & ' ''CRYSTAL'' USING FILE:',/X,A72) 131 format(/' INITIAL CONFIGURATION IS READ FROM FILE:',/X,A72) 132 format(/' LATTICE DEFECTS ARE GENERATED USING FILE:',/X,A72) 133 format(' DEFORMATION OF THE BLOCK BY FACTORS:', & /' X =', D15.7,' Y =', D15.7,' Z =', D15.7, & /' RBD =',D15.7) 134 format(/3X, ' ...... NO RELAXATION OF THE BLOCK ......') 135 format(//' TOO MANY ATOMS IN THE BLOCK (ntot >', & I8,' ) !!!'//) 136 format(/' CRACK INFORMATION IS READ FROM FILE:',/X,A72,/) 137 format(I3,t5,D11.4,t18,D14.7,' X =',D11.4, & ' Y =',D11.4,t64,i2) 138 format(2x,'N',t10,'K_I',t24,'ENERGY',t38,'CRACK TIP POSITION', & t64,'NF') 300 format(/' =======================================', &/ ' | OUTPUT OF PROGRAM "STATICS" |', &/ ' | --------------------------- |' &/ ' | (Yuri Mishin, Feb.-March 1996) |', &/ ' | [revised in June 1998] |', &/ ' =======================================', &/' For comments, suggestions and bugs contact mishin@vt.edu'/) 301 format(' REMEMBER: in this program the energy is measured', &' in eV,',/5x,' distances in Angstroems, forces in', &' eV/Angstroem,',/11x,' planar defect energy in eV/Angstroem^2') 302 format(4f12.6,i5) 400 format('#',4E18.7) 401 format('#',4I8) 402 format('#',E18.7,3I8) 403 format(A1,4E18.7) 404 format(A1,4I8) 405 format(A1,E18.7,3I8) STOP END SUBROUTINE CRYSTAL(co) *----------------------------------------- * Yuri Mishin, Feb. 1996; June 1998 *=================================================================== * Subroutine CRYSTAL establishes essential crystallographic * characteristics of a given lattice for a given orientation * of a rectangular simulation block, such as * - the lattice period in directions parallel to the block axes, * - d-spacing in directions parallel to the block axes, * - the number of Bravais lattice planes before the lattice * repeats in directions parallel to the block axes, etc. * It embeds this inner block in a bigger outer block according to * to a given search radius. It then locates free atoms (in the * inner block) as well as buffer atoms and fixed atoms (in the * outer block in case of fixed boundary conditions). For periodic * boundary conditions a table of atomic displacements is * constructed. * * * Subroutines used: FACTOR * * Part of this subroutine is based on routines LATTIS and REGION * developed by M. J. Norgett in 1972 (DEVIL code) *================================================================== * Illustration of the block geometry for fixed boundary conditions * X and Y * ^ Y * | * | * ------------------------- * | FIXED | * | ----------------- | * | | BUFFER | | * | | B ---------- B | | * | F | U | | U | F | * | I | F | FREE | F | I | * | X | F | | F | X | ---> X * | E | E | | E | E | * | D | R ---------- R | D | * | | BUFFER | | * | ----------------- | * | FIXED | * ------------------------- * *================================================================ implicit real*8 (a-h,o-z) logical bcx,bcy,bcz,bx,by,bz parameter (nspmax = 4, ! max number of chem. species & nbasmax = 30000, ! max number of basis atoms & maxatom = 100000 ! max number of atoms in block &) character*4 element(nspmax) common /accuracy/ acc common /chemical/ cutoff2,np(nspmax,nspmax),numat,npairs,npr, & element common /basis/ xbas(nbasmax),ybas(nbasmax),zbas(nbasmax), & lbas(nbasmax),nbas common /block/ search2,xmn,ymn,zmn,xmx,ymx,zmx, & xtmn,ytmn,ztmn,xtmx,ytmx,ztmx, & xdisp,ydisp,zdisp,mdx,mdy,mdz, & ntot,nb,nfree common /atoms/ lco(maxatom) common /cyc/ bcx,bcy,bcz dimension xlat(3,3),xrcp(3,3),xmi(3,3),xmap(3,3) dimension co(*), co1(3*maxatom), lco1(maxatom), & ind(maxatom),isort(maxatom) ********************************************************* * Input lattice translation vectors stored by columns * ********************************************************* do j=1,3 read(5,*) xlat(1,j), xlat(2,j), xlat(3,j) enddo write(6, 100) ((xlat(i,j),i=1,3),j=1,3) ***************************************************************** * Calculate the reciprocal lattice vectors. Vc is the unit cell * * volume. * ***************************************************************** xrcp(1,1) = xlat(2,2)*xlat(3,3) - xlat(2,3)*xlat(3,2) xrcp(2,1) = xlat(2,3)*xlat(3,1) - xlat(2,1)*xlat(3,3) xrcp(3,1) = xlat(2,1)*xlat(3,2) - xlat(2,2)*xlat(3,1) Vc = xrcp(1,1)*xlat(1,1) + xrcp(2,1)*xlat(1,2) + & xrcp(3,1)*xlat(1,3) xrcp(1,1) = xrcp(1,1)/Vc xrcp(2,1) = xrcp(2,1)/Vc xrcp(3,1) = xrcp(3,1)/Vc xrcp(1,2) = (xlat(3,2)*xlat(1,3) - xlat(3,3)*xlat(1,2))/Vc xrcp(1,3) = (xlat(1,2)*xlat(2,3) - xlat(1,3)*xlat(2,2))/Vc xrcp(2,2) = (xlat(3,3)*xlat(1,1) - xlat(3,1)*xlat(1,3))/Vc xrcp(2,3) = (xlat(1,3)*xlat(2,1) - xlat(1,1)*xlat(2,3))/Vc xrcp(3,2) = (xlat(3,1)*xlat(1,2) - xlat(3,2)*xlat(1,1))/Vc xrcp(3,3) = (xlat(1,1)*xlat(2,2) - xlat(1,2)*xlat(2,1))/Vc Vc = DABS(Vc) ********************************************** * Input basis vectors and atomic species * ********************************************** read(5,*) nbas write(6,101) nbas if (nbas.gt.nspmax) then write(6, 102) nspmax stop endif if (nbas.lt.0) then write(6, 121) stop endif write(6,104) do nl = 1,nbas read (5,*) xbas(nl),ybas(nl),zbas(nl),lbas(nl) write(6, 105) xbas(nl),ybas(nl),zbas(nl),element(lbas(nl)) enddo ******************************************************** * Input lattice block axes. Vectors stored by columns * ******************************************************** do j=1,3 read(5,*) xmi(1,j), xmi(2,j), xmi(3,j) enddo write(6, 106) ((xmi(i,j),i=1,3),j=1,3) ************************************************************ * Check if the block axes are orthogonal and normalize * * their lengths * ************************************************************ if( & (DABS(xmi(1,1)*xmi(1,2)+xmi(2,1)*xmi(2,2)+xmi(3,1)*xmi(3,2)) & .gt.acc). or. & (DABS(xmi(1,1)*xmi(1,3)+xmi(2,1)*xmi(2,3)+xmi(3,1)*xmi(3,3)) & .gt.acc). or. & (DABS(xmi(1,2)*xmi(1,3)+xmi(2,2)*xmi(2,3)+xmi(3,2)*xmi(3,3)) & .gt.acc)) then write (6, 107) stop endif xm2 = xmi(1,1)*xmi(1,1) + xmi(2,1)*xmi(2,1) + xmi(3,1)*xmi(3,1) ym2 = xmi(1,2)*xmi(1,2) + xmi(2,2)*xmi(2,2) + xmi(3,2)*xmi(3,2) zm2 = xmi(1,3)*xmi(1,3) + xmi(2,3)*xmi(2,3) + xmi(3,3)*xmi(3,3) xm1 = SQRT(xm2) ym1 = SQRT(ym2) zm1 = SQRT(zm2) xmi(1,1) = xmi(1,1)/xm1 xmi(2,1) = xmi(2,1)/xm1 xmi(3,1) = xmi(3,1)/xm1 xmi(1,2) = xmi(1,2)/ym1 xmi(2,2) = xmi(2,2)/ym1 xmi(3,2) = xmi(3,2)/ym1 xmi(1,3) = xmi(1,3)/zm1 xmi(2,3) = xmi(2,3)/zm1 xmi(3,3) = xmi(3,3)/zm1 ***************************************************************** * Switch over to the reference system fixed to the block axes. * * During this transformation we determine projections of the * * lattice translation vectors onto the block axis vectors. * * We also find the largest real factors (xl, yl and zl) for * * those projections. Those factors are equal the d-spacings of * * the Bravais lattice in directions parallel to the block axes * ***************************************************************** do nl = 1,nbas x1 = xbas(nl)*xmi(1,1)+ybas(nl)*xmi(2,1)+zbas(nl)*xmi(3,1) x2 = xbas(nl)*xmi(1,2)+ybas(nl)*xmi(2,2)+zbas(nl)*xmi(3,2) x3 = xbas(nl)*xmi(1,3)+ybas(nl)*xmi(2,3)+zbas(nl)*xmi(3,3) xbas(nl) = x1 ybas(nl) = x2 zbas(nl) = x3 enddo x1 = xlat(1,1)*xmi(1,1)+xlat(2,1)*xmi(2,1)+xlat(3,1)*xmi(3,1) x2 = xlat(1,2)*xmi(1,1)+xlat(2,2)*xmi(2,1)+xlat(3,2)*xmi(3,1) x3 = xlat(1,3)*xmi(1,1)+xlat(2,3)*xmi(2,1)+xlat(3,3)*xmi(3,1) call FACTOR(x1,x2,x12) call FACTOR(x12,x3,xl) x1 = xlat(1,1)*xmi(1,2)+xlat(2,1)*xmi(2,2)+xlat(3,1)*xmi(3,2) x2 = xlat(1,2)*xmi(1,2)+xlat(2,2)*xmi(2,2)+xlat(3,2)*xmi(3,2) x3 = xlat(1,3)*xmi(1,2)+xlat(2,3)*xmi(2,2)+xlat(3,3)*xmi(3,2) call FACTOR(x1,x2,x12) call FACTOR(x12,x3,yl) x1 = xlat(1,1)*xmi(1,3)+xlat(2,1)*xmi(2,3)+xlat(3,1)*xmi(3,3) x2 = xlat(1,2)*xmi(1,3)+xlat(2,2)*xmi(2,3)+xlat(3,2)*xmi(3,3) x3 = xlat(1,3)*xmi(1,3)+xlat(2,3)*xmi(2,3)+xlat(3,3)*xmi(3,3) call FACTOR(x1,x2,x12) call FACTOR(x12,x3,zl) ************************************************************* * Project the reciprocal lattice vectors onto the block * * axes and multiply them by factors xl, yl and zl * ************************************************************* xmap(1,1) = xl * (xrcp(1,1)*xmi(1,1) + xrcp(1,2)*xmi(2,1) + & xrcp(1,3)*xmi(3,1)) xmap(2,1) = xl * (xrcp(2,1)*xmi(1,1) + xrcp(2,2)*xmi(2,1) + & xrcp(2,3)*xmi(3,1)) xmap(3,1) = xl * (xrcp(3,1)*xmi(1,1) + xrcp(3,2)*xmi(2,1) + & xrcp(3,3)*xmi(3,1)) xmap(1,2) = yl * (xrcp(1,1)*xmi(1,2) + xrcp(1,2)*xmi(2,2) + & xrcp(1,3)*xmi(3,2)) xmap(2,2) = yl * (xrcp(2,1)*xmi(1,2) + xrcp(2,2)*xmi(2,2) + & xrcp(2,3)*xmi(3,2)) xmap(3,2) = yl * (xrcp(3,1)*xmi(1,2) + xrcp(3,2)*xmi(2,2) + & xrcp(3,3)*xmi(3,2)) xmap(1,3) = zl * (xrcp(1,1)*xmi(1,3) + xrcp(1,2)*xmi(2,3) + & xrcp(1,3)*xmi(3,3)) xmap(2,3) = zl * (xrcp(2,1)*xmi(1,3) + xrcp(2,2)*xmi(2,3) + & xrcp(2,3)*xmi(3,3)) xmap(3,3) = zl * (xrcp(3,1)*xmi(1,3) + xrcp(3,2)*xmi(2,3) + & xrcp(3,3)*xmi(3,3)) ***************************************************************** * Determine the largest real factors (xr, yr, zr) for the * * projections of the reciprocal lattice vectors onto the block * * axes. These factors equal the Bravais lattice periods in the * * directions parallel to the block axes. In fact, using the * * array xmap we determine the combined factors * * * * x1=xl*xr, y1=yl*yr and z1=zl*zr * * * ***************************************************************** call FACTOR(xmap(1,1),xmap(2,1),x12) call FACTOR(xmap(1,2),xmap(2,2),y12) call FACTOR(xmap(1,3),xmap(2,3),z12) call FACTOR(x12,xmap(3,1),x1) call FACTOR(y12,xmap(3,2),y1) call FACTOR(z12,xmap(3,3),z1) ****************************************************************** * Determine the lattice periodicity in the block axes directions * ****************************************************************** ndx = INT(1.0D0/x1 + 0.5D0) ndy = INT(1.0D0/y1 + 0.5D0) ndz = INT(1.0D0/z1 + 0.5D0) write(6, 109) ndx,ndy,ndz ***************************************************************** * Determine the packing factor npack. In an imaginary lattice * * with a unit cell xl by yl by zl [the sides of this cell are * * parallel to the block axes] every npack-th site coincides * * with a site of the actual Bravais lattice * ***************************************************************** npack = INT(Vc/(xl * yl * zl) + 0.5D0) **************************************************************** * Input dimensions of the inner block and the search radius. * * Negative block dimensions correspond to cyclic boundary * * conditions. If search radius < 0, use the cut-off radius * * of the potential * **************************************************************** read(5,*) nfx, nfy, nfz read(5,*) search if (search.le.0.D0) search=DSQRT(cutoff2) write(6, 110) search search2 = search*search ******************************************** * Determine dimensions of the outer block. * ******************************************** nxadd = INT(search/xl) + 1 nyadd = INT(search/yl) + 1 nzadd = INT(search/zl) + 1 ntx = IABS(nfx) + 2*nxadd nty = IABS(nfy) + 2*nyadd ntz = IABS(nfz) + 2*nzadd ***************************************************************** * Prepare information required for block construction. * * For periodic boundary conditions in x-direction, find * * the number mdx of x-dimensions of the inner block (xdisp) * * required to reach and cover the border of the outer block, * * so that abs(2*nfx*mdx + nfx) > ntx. Variables mdy and mdz * * have a similar meaning * ***************************************************************** *----------- x-direction: if (nfx.gt.0) then mdx = 0 bcx = .false. nbx = ntx ntx = ((ntx + 2*nxadd)/ndx)*ndx + ndx else nfx = IABS(nfx) if (MOD(nfx,ndx).ne.0) then write(6, 201) stop endif mdx = (ntx + nfx - 1)/(nfx + nfx) bcx=.true. endif nx3 = nfx/2 nx2 = nx3 - nfx + 1 nx4 = ntx/2 nx1 = nx4 - ntx + 1 *----------- y-direction: if (nfy.gt.0) then mdy = 0 bcy = .false. nby = nty nty = ((nty + 2*nyadd)/ndy)*ndy + ndy else nfy = IABS(nfy) if (MOD(nfy,ndy).ne.0) then write(6, 202) stop endif mdy = (nty + nfy - 1)/(nfy + nfy) bcy = .true. endif ny3 = nfy/2 ny2 = ny3 - nfy + 1 ny4 = nty/2 ny1 = ny4 - nty + 1 *----------- z-direction: if (nfz.gt.0) then mdz = 0 bcz = .false. nbz = ntz ntz = ((ntz + 2*nzadd)/ndz)*ndz + ndz else nfz = IABS(nfz) if (MOD(nfz,ndz).ne.0) then write(6, 203) stop endif mdz = (ntz + nfz - 1)/(nfz + nfz) bcz = .true. endif nz3 = nfz/2 nz2 = nz3 - nfz + 1 nz4 = ntz/2 nz1 = nz4 - ntz + 1 ************************************************** * Print out the inner and outer block dimensions * ************************************************** xdisp = nfx*xl ydisp = nfy*yl zdisp = nfz*zl if (bcx.and.bcy.and.bcz) go to 1 write(6,112) ntx, nty, ntz write(6,113) nxadd, nyadd, nzadd 1 write(6,111) nfx, nfy, nfz write(6,117) xdisp, ydisp, zdisp if (bcx) write(6,114) if (bcy) write(6,115) if (bcz) write(6,116) *********************************************************** * Prepare the block sizes and other information required * * for the block construction * *********************************************************** c--- For X-direction ----------- mx1 = nfx/2 mx2 = nfx - mx1 xmx = xl*mx1 - acc xmn = -xl*mx2 - acc if(.not.bcx) then ntx1 = ntx/2 ntx2 = ntx - ntx1 xtmx = xl*ntx1 - acc xtmn = -xl*ntx2 - acc nbx1 = nbx/2 nbx2 = nbx - nbx1 xbmx = xl*nbx1 - acc xbmn = -xl*nbx2 - acc else xtmx = xmx xtmn = xmn xbmx = xmx xbmn = xmn endif c--- For Y-direction ----------- my1 = nfy/2 my2 = nfy - my1 ymx = yl*my1 - acc ymn = -yl*my2 - acc if(.not.bcy) then nty1 = nty/2 nty2 = nty - nty1 ytmx = yl*nty1 - acc ytmn = -yl*nty2 - acc nby1 = nby/2 nby2 = nby - nby1 ybmx = yl*nby1 - acc ybmn = -yl*nby2 - acc else ytmx = ymx ytmn = ymn ybmx = ymx ybmn = ymn endif c--- For Z-direction ----------- mz1 = nfz/2 mz2 = nfz - mz1 zmx = zl*mz1 - acc zmn = -zl*mz2 - acc if(.not.bcz) then ntz1 = ntz/2 ntz2 = ntz - ntz1 ztmx = zl*ntz1 - acc ztmn = -zl*ntz2 - acc nbz1 = nbz/2 nbz2 = nbz - nbz1 zbmx = zl*nbz1 - acc zbmn = -zl*nbz2 - acc else ztmx = zmx ztmn = zmn zbmx = zmx zbmn = zmn endif ******************************* * Construction of the block * ******************************* zrpt = ndz*zl i = nx1 - 1 xco = i*xl x1 = i*xmap(1,1) y1 = i*xmap(2,1) z1 = i*xmap(3,1) 11 i = i + 1 if (i.gt.nx4) goto 28 xco = xco + xl x1 = x1 + xmap(1,1) y1 = y1 + xmap(2,1) z1 = z1 + xmap(3,1) j = ny1 - 1 yco = j*yl x2 = x1 + j*xmap(1,2) y2 = y1 + j*xmap(2,2) z2 = z1 + j*xmap(3,2) jplus = 1 i1 = i i2 = 0 bx = (i.lt.nx2).or.(i.gt.nx3) if (bx.and.bcx) goto 12 goto 14 12 continue if (i1.ge.nx2) goto 13 i2 = i2 - 1 i1 = i1 + nfx goto 12 13 continue if (i1.le.nx3) goto 14 i2 = i2 + 1 i1 = i1 - nfx goto 13 14 j = j + jplus if (j.gt.ny4) goto 11 yco = yco + jplus*yl x2 = x2 + jplus*xmap(1,2) y2 = y2 + jplus*xmap(2,2) z2 = z2 + jplus*xmap(3,2) k = nz1 - 1 zco = k*zl x3 = x2 + k*xmap(1,3) y3 = y2 + k*xmap(2,3) z3 = z2 + k*xmap(3,3) j1 = j j2 = 0 by = (j.lt.ny2).or.(j.gt.ny3) if (by.and.bcy) goto 15 goto 17 15 continue if (j1.ge.ny2) goto 16 j2 = j2 - 1 j1 = j1 + nfy goto 15 16 continue if (j1.le.ny3) goto 17 j2 = j2 + 1 j1 = j1 - nfy goto 16 17 k = k + 1 if ((k - nz1).gt.ndz) goto 14 zco = zco + zl x3 = x3 + xmap(1,3) y3 = y3 + xmap(2,3) z3 = z3 + xmap(3,3) i3 = INT(x3 + DSIGN(0.5D0,x3)) j3 = INT(y3 + DSIGN(0.5D0,y3)) k3 = INT(z3 + DSIGN(0.5D0,z3)) if (DABS(x3 - i3).gt.acc) goto 17 if (DABS(y3 - j3).gt.acc) goto 17 if (DABS(z3 - k3).gt.acc) goto 17 jplus = npack/ndz 18 continue if (k.gt.nz4) goto 14 k1 = k k2 = 0 bz = (k.lt.nz2).or.(k.gt.nz3) if (bz.and.bcz) goto 19 goto 21 19 continue if (k1.ge.nz2) goto 20 k2 = k2 - 1 k1 = k1 + nfz goto 19 20 continue if (k1.le.nz3) goto 21 k2 = k2 + 1 k1 = k1 - nfz goto 20 21 continue do 27 nl1 = 1,nbas x4 = xco + xbas(nl1) y4 = yco + ybas(nl1) z4 = zco + zbas(nl1) if (bx.or.by.or.bz) goto 23 nfree = nfree + 1 k4 = 3*nfree j4 = k4 - 1 i4 = j4 - 1 ntot = ntot + 1 if (ntot.gt.maxatom) then write(6, 120) maxatom stop endif if (nfree.eq.ntot) goto 22 k3 = 3*ntot j3 = k3 - 1 i3 = j3 - 1 co1(i3) = co1(i4) co1(j3) = co1(j4) co1(k3) = co1(k4) lco1(ntot) = lco1(nfree) 22 continue co1(i4) = x4 co1(j4) = y4 co1(k4) = z4 lco1(nfree) = lbas(nl1) go to 27 23 continue if (bx.and.bcx) goto 27 if (by.and.bcy) goto 27 if (bz.and.bcz) goto 27 ntot = ntot + 1 if (ntot.gt.maxatom) then write(6, 120) maxatom stop endif k4 = 3*ntot j4=k4-1 i4 = j4 - 1 co1(i4) = x4 co1(j4) = y4 co1(k4) = z4 lco1(ntot) = lbas(nl1) 27 continue k = k + ndz zco = zco + zrpt goto 18 28 continue ********************************************* * Smoothen out the block borders * * and center the block at x=0, y=0 and z=0 * ********************************************* ntot3 = 3*ntot if(.not.bcx) xd = xtmx - xtmn if(.not.bcy) yd = ytmx - ytmn if(.not.bcz) zd = ztmx - ztmn c--- For X-direction ----------- do n=3,ntot3,3 ncx = n - 2 if(co1(ncx).gt.xtmx.or.co1(ncx).lt.xtmn) then if (bcx) then do i = -mdx,mdx dsp = xdisp*i if(DABS(dsp).gt.acc) then xc = co1(ncx) + dsp if((xc.ge.xtmn).and.(xc.le.xtmx)) then co1(ncx) = xc go to 31 endif endif enddo write(6,210) write(6,211) stop else do i=-1,1,2 xc = co1(ncx) + xd*i if((xc.ge.xtmn).and.(xc.le.xtmx)) then co1(ncx) = xc go to 31 endif enddo write(6,210) write(6,211) stop endif 31 continue endif enddo c--- For Y-direction ----------- do n=3,ntot3,3 ncy = n - 1 if(co1(ncy).gt.ytmx.or.co1(ncy).lt.ytmn) then if (bcy) then do i=-mdy,mdy dsp = ydisp*i if(DABS(dsp).gt.acc) then yc = co1(ncy) + dsp if((yc.ge.ytmn).and.(yc.le.ytmx)) then co1(ncy) = yc go to 32 endif endif enddo write(6,210) write(6,212) stop else do i=-1,1,2 yc = co1(ncy) + yd*i if((yc.ge.ytmn).and.(yc.le.ytmx)) then co1(ncy) = yc go to 32 endif enddo write(6,210) write(6,212) stop endif 32 continue endif enddo c--- For Z-direction ----------- do n=3,ntot3,3 ncz = n if(co1(ncz).gt.ztmx.or.co1(ncz).lt.ztmn) then if (bcz) then do i=-mdz,mdz dsp = zdisp*i if(DABS(dsp).gt.acc) then zc = co1(ncz) + dsp if((zc.ge.ztmn).and.(zc.le.ztmx)) then co1(ncz) = zc go to 33 endif endif enddo write(6,210) write(6,213) stop else do i=-1,1,2 zc = co1(ncz) + zd*i if((zc.ge.ztmn).and.(zc.le.ztmx)) then co1(ncz) = zc go to 33 endif enddo write(6,210) write(6,213) stop endif 33 continue endif enddo if(bcx.and.bcy.and.bcz) go to 41 ***************************************************** * For fixed boundary conditions smoothen out the * * free/buffer and buffer/fixed interfaces and * * center the free region at x=0, y=0 and z=0 * ***************************************************** * re-sort the atoms into 3 classes: free (ind=0), * buffer (ind=1) and fixed (ind=2) do iat = 1,ntot k = 3*iat z = co1(k) y = co1(k-1) x = co1(k-2) if(x.le.xmx.and.x.ge.xmn.and. & y.le.ymx.and.y.ge.ymn.and. & z.le.zmx.and.z.ge.zmn) then ind(iat) = 0 go to 40 endif if(x.le.xbmx.and.x.ge.xbmn.and. & y.le.ybmx.and.y.ge.ybmn.and. & z.le.zbmx.and.z.ge.zbmn) then ind(iat) = 1 go to 40 endif ind(iat) = 2 40 enddo * Re-order arrays co1 and lco1 according to the new classes. * Array isort establishes the correspondence between the * two orders nfr = 0 nb = 0 nt = 0 * find all free atoms do iat = 1,ntot if(ind(iat).eq.0) then nt = nt + 1 nb = nb + 1 nfr = nfr + 1 isort(nt) = iat endif enddo * find all buffer atoms do iat = 1,ntot if(ind(iat).eq.1) then nt = nt + 1 nb = nb + 1 isort(nt) = iat endif enddo * find all fixed atoms do iat = 1,ntot if(ind(iat).eq.2) then nt = nt + 1 isort(nt) = iat endif enddo * re-write arrays co and lco according to the new order do i = 1, ntot j = isort(i) i3 = 3*i j3 = 3*j co(i3) = co1(j3) co(i3-1) = co1(j3-1) co(i3-2) = co1(j3-2) lco(i) = lco1(j) enddo go to 42 41 do i=1,ntot3 co(i) = co1(i) enddo do i=1,ntot lco(i) = lco1(i) enddo 42 write(6,119) ntot, nfree if(bcx.and.bcy.and.bcz) then nb = nfree go to 43 endif write(6,204) nb - nfree write(6,205) ntot-nb 43 write(6,206) *--------------- FORMAT statements ---------------------------- 100 format(/' PRIMITIVE TRANSLATION VECTORS:',3(/ 3f12.5)) 101 format(/' NUMBER OF BASIS ATOMS = ',i2) 102 format(//,' ***** TOO MANY BASIS ATOMS ( >', I3,' ) *****'//) 104 format(6x,' COORDINATES',24x,' ELEMENT') 105 format(6x, 3(f10.6,2x), 4x,a4) 106 format(/' EDGE VECTORS OF THE SIMULATION BLOCK:',3(/ 3f12.5)) 107 format(//' EDGE VECTORS OF THE BLOCK ARE NOT ORTHOGONAL !!!'//) 109 format(/' NUMBER OF PLANES BEFORE LATTICE REPEATS IN BLOCK', &' DIRECTIONS:',/6x,3(i3,7x)) 110 format(/5x,'RADIUS OF SEARCH = ',f7.4/) 111 format(5x,'INNER BLOCK IS ',i3,' BY ',i3,' BY ',i3) 112 format(5x,'OUTER BLOCK IS ',i3,' BY ',i3,' BY ',i3) 113 format(5x,'BUFFER THICKNESS ',i3,' BY ',i3,' BY ',i3) 114 format(5x,'CYCLIC BOUNDARY IN X DIRECTION') 115 format(5x,'CYCLIC BOUNDARY IN Y DIRECTION') 116 format(5x,'CYCLIC BOUNDARY IN Z DIRECTION') 117 format(5x,'DIMENSIONS OF THE INNER BLOCK:',3f10.5,/) 118 format(//x,'DIMENSION OF ARRAY DISP (',I5, & ') IS TOO SMALL !!!'//) 119 format(/,' THE PERFECT LATTICE BLOCK CONTAINS:',/, & ' TOTAL NUMBER OF ATOMS: NTOT =',i7, &/' NUMBER OF FREE ATOMS: NFREE =',i7) 120 format(//' TOO MANY ATOMS IN THE BLOCK (ntot >', & I8,' ) !!!'//) 121 format(/' WRONG NUMBER OF BASIS ATOMS!!!'/) 201 format(//' CYCLIC BOUNDARY CONDITIONS IN X-DIRECTION'/ & ' ARE INCOMPATIBLE WITH THE SPECIFIED BLOCK SIZE !!!'//) 202 format(//' CYCLIC BOUNDARY CONDITIONS IN Y-DIRECTION'/ & ' ARE INCOMPATIBLE WITH THE SPECIFIED BLOCK SIZE !!!'//) 203 format(//' CYCLIC BOUNDARY CONDITIONS IN Z-DIRECTION'/ & ' ARE INCOMPATIBLE WITH THE SPECIFIED BLOCK SIZE !!!'//) 204 format(' NUMBER OF BUFFER ATOMS:',i16) 205 format(' NUMBER OF FIXED ATOMS:',i16) 206 format(/, & ' ---------------------------------------------------'/) 210 format(//' THE ROUTINE FAILS TO CORRECT THE BLOCK SHAPE') 211 format(' THERE IS A PROBLEM IN X-DIRECTION !!!'//) 212 format(' THERE IS A PROBLEM IN Y-DIRECTION !!!'//) 213 format(' THERE IS A PROBLEM IN Z-DIRECTION !!!'//) RETURN END SUBROUTINE DEFECT (co) *----------------------------------------- * Yuri Mishin, Feb. 1996, June 1998 *===================================================================== * Subroutine DEFECT introduces point defects (vacancies, antisites * or interstitials) or a planar defect (e.g. grain boundary, SF or * APB) in X-Z plane. Vacancies and antisites are referred to by * numbers of lattice sites; interstitial atoms are given by their * coordinates. A symmetrical tilt grain boundary is produced by * a 180 deg rotation of the upper grain (y>0) around axis Y. * A planar fault is produced by a given rigid body translation of * the upper grain relative to the lower one. If a planar defect * was introduced, the rectangular shape of the block is restored * by applying translations in X and Z. The defects are introduced * in the following order: * 1. grain boundary, 2. planar fault, 3. vacancies, 4. antisites, * 5. interstitials *===================================================================== implicit real*8 (a-h,o-z) logical ipoint,iplane parameter (nspmax = 4, ! max number of chem. species & maxatom = 100000,! max number of atoms in block & maxint = 100 ! max number of interstitials &) character*4 element(nspmax) common /accuracy/ acc common /chemical/ cutoff2,np(nspmax,nspmax),numat,npairs,npr, & element common /block/ search2,xmn,ymn,zmn,xmx,ymx,zmx, & xtmn,ytmn,ztmn,xtmx,ytmx,ztmx, & xdisp,ydisp,zdisp,mdx,mdy,mdz, & ntot,nb,nfree common /atoms/ lco(maxatom) common /defct/ nvac,ilat,ipoint,iplane dimension co(*), na(nspmax), xint(3,maxint),lcoint(maxint) *=== 180 degree rotation to create a grain boundary if requested === not = 3*ntot read(5,*) igb if (igb.gt.0) then iplane=.true. do n=3,not,3 if(co(n-1).gt.acc) then co(n) = -co(n) co(n-2)= -co(n-2) endif enddo write(6, 115) endif *=== Rigid-body displacement to create a planar fault read(5,*) ipl if (ipl.ne.0) then iplane=.true. read(5,*) rigidx, rigidy, rigidz do n = 3,not,3 if(co(n-1).gt.acc) then co(n-2) = co(n-2) + rigidx co(n-1) = co(n-1) + rigidy co(n) = co(n) + rigidz endif enddo write(6, 116) rigidx, rigidy, rigidz endif *************************************************************** * Restore the rectangular shape of the simulation block after * * the introduction of a planar defect * *************************************************************** if (iplane) then c--- For X-direction ----------- do n=3,not,3 ncx = n - 2 if(co(ncx).gt.xtmx.or.co(ncx).lt.xtmn) then do i = -mdx,mdx dsp = xdisp*i if(DABS(dsp).gt.acc) then xc = co(ncx) + dsp if((xc.ge.xtmn).and.(xc.le.xtmx)) then co(ncx) = xc go to 31 endif endif enddo write(6,110) write(6,111) stop 31 continue endif enddo c--- For Z-direction ----------- do n=3,not,3 ncz = n if(co(ncz).gt.ztmx.or.co(ncz).lt.ztmn) then do i = -mdz,mdz dsp = zdisp*i if(DABS(dsp).gt.acc) then zc = co(ncz) + dsp if((zc.ge.ztmn).and.(zc.le.ztmx)) then co(ncz) = zc go to 32 endif endif enddo write(6,110) write(6,112) stop 32 continue endif enddo endif *====== Introduce vacancies if requested ============= read(5,*) nvac if (nvac.lt.0.or.nvac.gt.nfree) then write(6,103) nvac stop endif if (nvac.ne.0) then ipoint= .true. if (nvac.eq.1) then write(6,100) else write(6,101) nvac endif do i=1,nvac read(5,*) ivac if (ivac.lt.1.or.ivac.gt.nfree) then write(6,102) ivac stop endif write(6, 121) ivac lco(ivac)=0 enddo endif *===== Introduce antisites if requested =========== read(5,*) nanti if (nanti.lt.0.or.nanti.gt.nfree) then write(6,104) nanti stop endif if (nanti.ne.0) then ipoint=.true. if (nanti.eq.1) then write(6,105) else write(6,106) nanti endif do i=1,nanti read(5,*) ianti, new if (ianti.lt.1.or.ianti.gt.nfree) then write(6,107) ianti stop endif if (new.lt.1.or.new.gt.numat) then write(6,108) new, ianti stop endif write(6,109) ianti, lco(ianti), & element(lco(ianti)), new, element(new) lco(ianti) = new enddo endif *====== Introduce interstitials if requested ============= read(5,*) nint if (nint.gt.maxint) then write(6,123) nint,maxint stop endif if (nint.lt.0) then write(6,128) nint stop endif if (nint.gt.0) then ipoint= .true. if (nint.eq.1) then write(6,124) else write(6,125) nint endif do i=1,nint read(5,*) (xint(j,i), j=1,3), lcoint(i) if (xint(1,i).gt.xmx.or.xint(1,i).lt.xmn.or. & xint(2,i).gt.ymx.or.xint(2,i).lt.ymn.or. & xint(3,i).gt.zmx.or.xint(3,i).lt.zmn) then write(6,126) i stop endif if (lcoint(i).lt.1.or.lcoint(i).gt.numat) then write(6,129) lcoint(i) stop endif write(6, 127) i,(xint(j,i), j=1,3), lcoint(i), & element(lcoint(i)) enddo endif * if there are new atoms, update arrays co and lco if(nint.gt.0) then ntot = ntot + nint if(ntot.gt.maxatom) then write(6,130) maxatom stop endif nb = nb + nint nfree = nfree + nint nint3 = 3*nint do i = ntot, nfree + 1, -1 lco(i) = lco(i - nint) j = 3*i co(j) = co(j- nint3) co(j-1) = co(j-1-nint3) co(j-2) = co(j-2-nint3) enddo do i = 1, nint j = nfree - nint + i k = 3*j co(k) = xint(3,i) co(k-1) = xint(2,i) co(k-2) = xint(1,i) lco(j) = lcoint(i) enddo endif *************************************** * Recalculate the block composition * *************************************** if(ipoint.or.ilat.eq.0) then do i=1,numat na(i)=0 enddo numvac=0 do i=1, nfree ilco=lco(i) if (ilco.ne.0) then j=ilco na(j) = na(j) + 1 else numvac = numvac + 1 endif enddo nn=0 do i=1,numat nn = nn + na(i) enddo if ((nn+numvac.ne.nfree.or.numvac.ne.nvac).and. & ilat.ne.0) then write(6,117) stop else write(6,118) nfree-numvac, numvac do i=1,numat write(6,119) na(i), element(i) enddo write(6,120) endif endif *------------- FORMAT statements ------------------------------ 100 format(/' THE BLOCK CONTAINS A VACANCY') 101 format(/' THE BLOCK CONTAINS',i4,' VACANCIES:') 102 format(//' SORRY, CANNOT CREATE A VACANCY.', &/' SITE NO.',i6,' LIES OUTSIDE THE BLOCK !!!'//) 103 format(//' WRONG NUMBER OF VACANCIES:', i5, ' !!!'//) 104 format(//' WRONG NUMBER OF ANTISITES:', i5, ' !!!'//) 105 format(/' THE BLOCK CONTAINS AN ANTISITE') 106 format(/' THE BLOCK CONTAINS',i4,' ANTISITES:') 107 format(//' SORRY, CANNOT CREATE AN ANTISITE.', &/' SITE NO.',i6,' LIES OUTSIDE THE BLOCK !!!'//) 108 format(//' SORRY, CANNOT CREATE AN ANTISITE', &/' WITH WRONG SORT =',i3, ' AT SITE NO.',i5, ' !!!'//) 109 format(' AT SITE',i5, ' WITH OLD SORT =', i2,' ( ',a4,')', & ' AND NEW SORT =', i2,' ( ',a4,')') 110 format(//' THE ROUTINE FAILS TO RESTORE THE BLOCK SHAPE', & /' AFTER INTRODUCTION OF A PLANAR DEFECT') 111 format(' THERE IS A PROBLEM IN X-DIRECTION !!!'//) 112 format(' THERE IS A PROBLEM IN Z-DIRECTION !!!'//) 115 format(/' A GRAIN BOUNDARY IN X-Z PLANE HAS BEEN CREATED', &/ ' BY A 180 DEGREE ROTATION OF THE UPPPER GRAIN') 116 format(/' A PLANAR DEFECT (in X-Z plane) HAS BEEN CREATED', &/' BY RIGID-BODY DISPLACEMENT OF THE UPPPER GRAIN BY VECTOR: ', &/ 3f12.5) 117 format(//' ATOMIC DISBALANCE IN THE BLOCK !!!'//) 118 format(/' ======= COMPOSITION OF THE DEFECTED BLOCK: =======', &//' THE BLOCK CONTAINS',i6,' ATOMS AND', &i4,' VACANCIES') 119 format(i6, ' ATOMS ARE ',a4) 120 format(2x,51('=')) 121 format(' AT SITE', i5) 123 format(//' TOO MANY INTERSTITIALS:', i4, ' >',I4,' !!!'//) 124 format(/' THE BLOCK CONTAINS AN INTERSTITIAL') 125 format(/' THE BLOCK CONTAINS',i4,' INTERSTITIALS:') 126 format(//' INTERSTITIAL', I4,' LIES OUTSIDE THE FREE ', & 'REGION !!!'//) 127 format(I4,X,' X =',F10.5,' Y =',F10.5,' Z =',F10.5, & ' SORT =',I2,' ( ',a4,')') 128 format(//' WRONG NUMBER OF INTERSTITIALS:', i5, ' !!!'//) 129 format(//' SORRY, CANNOT CREATE AN INTERSTITIAL', &/' WITH WRONG SORT =',i3) 130 format(//' THE INTERSTITIALS OVERFILL ARRAY co (ntot >', & I8,' ) !!!'//) RETURN END SUBROUTINE DEFORM (x,y,z,rbd,co,co1) *----------------------------------------- * Yuri Mishin, June 1998 *==================================================================== * Deformation of the simulation cell. * x, y, z are scaling factors in the respective directions * rbd is rigid-body displacement in Y direction (for planar defects) *==================================================================== implicit real*8 (a-h,o-z) common /block/ search2,xmn,ymn,zmn,xmx,ymx,zmx, & xtmn,ytmn,ztmn,xtmx,ytmx,ztmx, & xdisp,ydisp,zdisp,mdx,mdy,mdz, & ntot,nb,nfree common/minimize/ xdisp1,ydisp1,zdisp1,minmode,iter1D dimension co(*), co1(*) do i = 1,ntot k1 = 3*i j1 = k1 -1 i1 = j1 -1 co(i1) = x*co1(i1) co(j1) = y*co1(j1) co(k1) = z*co1(k1) if(i.gt.nfree) then if(co1(j1).gt.0.D0) co(j1) = co(j1) + rbd if(co1(j1).lt.0.D0) co(j1) = co(j1) - rbd endif enddo xdisp = x*xdisp1 ydisp = y*ydisp1 + 2.D0*rbd zdisp = z*zdisp1 RETURN END SUBROUTINE MODE1A(co,co1) *=================================================================== * This subroutine changes the atomic coordinates according to the * elasticity theory expressions for the displacement field around * a sharp mode I crack. * The original version of this subroutine was written by * Vijay Shastry, 3 Jan 1996 *=================================================================== implicit real*8(a-h,o-z) real*8 KI, mu, nu common /block/ search2,xmn,ymn,zmn,xmx,ymx,zmx, & xtmn,ytmn,ztmn,xtmx,ytmx,ztmx, & xdisp,ydisp,zdisp,mdx,mdy,mdz, & ntot,nb,nfree common /crack/ x0,y0,KI,mu,nu dimension co(*),co1(*) parameter (pi=3.141592654D0) Xk = 3.D0 - 4.D0*nu do i= 1,3*ntot,3 x = co1(i) - x0 y = co1(i+1) - y0 r = DSQRT(x**2 + y**2) ! distance from crack tip * Theta is the angle made by vector (x,y) with the crack plane theta = DATAN(y/x) if (x.lt.0.D0.and.y.gt.0.D0) then theta = theta + pi elseIf (x.lt.0.D0.and.y.le.0.D0) then theta = theta - pi endif * Elastic displacements around the crack ux = 0.5D0*KI/mu*DSQRT(0.5D0*r/Pi)*DCOS(0.5D0*theta) & *(Xk - 1.D0 + 2.D0*(DSIN(0.5D0*theta))**2) uy = 0.5D0*KI/mu*DSQRT(0.5D0*r/Pi)*DSIN(0.5D0*theta) & *(Xk + 1.D0 - 2.D0*(DCOS(0.5D0*theta))**2) * Update the atomic coordinates co(i) = co1(i) + ux co(i+1) = co1(i+1) + uy enddo RETURN END SUBROUTINE LOCATE(r1,nmin,co,xnew,ynew) *----------------------------------------- * Yuri Mishin, June 1998 *==================================================================== * This subroutine finds the current location of a crack tip. * It first finds all atoms located at the crack surface. * Such atoms are identified by their reduced coordination number. * Namely, an atom is recognized to be a surface atom if its first * coordination number (number of atoms within a search sphere of * radius r1, which is slightly larger than the first coodination * raduis in the perfect lattice) is smaller than a chosen number * nmin (e.g. 9-10 for fcc structure). The subroutine then finds * the crack tip from the position of the surface atom with the * largest X-coordinate *==================================================================== implicit real*8(a-h,o-z) parameter (maxatom = 100000,! max number of atoms in block & maxnb = 150 ! max number of neighbours & , nspmax=4) logical realn character*4 element(nspmax) common /table1/ neighbrs(maxatom,maxnb),idisp(maxatom,maxnb,3) common /table2/ numnbr(maxatom),realn(maxatom,maxnb) common /block/ search2,xmn,ymn,zmn,xmx,ymx,zmx, & xtmn,ytmn,ztmn,xtmx,ytmx,ztmx, & xdisp,ydisp,zdisp,mdx,mdy,mdz, & ntot,nb,nfree common/chemical/cutoff2,np(nspmax,nspmax),numat,npairs,npr, & element dimension co(*),nnn(maxatom) common/atoms/ lco(maxatom) r2 = r1*r1 do i=1,nb nnn(i) = 0 enddo do n1=1,nb k3 = 3*n1 j3 = k3-1 i3 = j3 -1 z1 = co(k3) y1 = co(j3) x1 = co(i3) do nbr = 1,numnbr(n1) n2 = neighbrs(n1,nbr) k4 = 3*n2 j4 = k4 -1 i4 = j4 -1 z2 = co(k4) y2 = co(j4) x2 = co(i4) if (.not.realn(n1,nbr)) then z2 = co(k4) + zdisp*idisp(n1,nbr,3) y2 = co(j4) + ydisp*idisp(n1,nbr,2) x2 = co(i4) + xdisp*idisp(n1,nbr,1) endif x12 = x1 - x2 y12 = y1 - y2 z12 = z1 - z2 u = x12*x12 + y12*y12 + z12*z12 if(u.gt.r2) go to 2 nnn(n1) = nnn(n1) + 1 if(realn(n1,nbr).and.n2.le.nb) nnn(n2) = nnn(n2) + 1 2 enddo enddo 1 format (i4,3f14.6) open(89, file='surface.dat',status='unknown') xnew = -1.D10 nscount =0 do n1=1,nb k3 = 3*n1 j3 = k3-1 i3 = j3 -1 x1 = co(i3) if (nnn(n1).lt.nmin) then write(89,1) & lco(n1), co(i3),co(j3),co(k3) nscount= nscount+1 endif if (nnn(n1).lt.nmin.and.x1.gt.xnew) then xnew = x1 ynew = co(j3) endif enddo rewind (89) open(99, file='surface.xyz',status='unknown') write(99,*) nscount write(99,*) nscount do i=1,nscount read(89,*) lc, a,b,c write(99,*) element(lc), a,b,c enddo close(99) close(89) RETURN END FUNCTION FMIN(ax,bx,tol,co,co1) c================================================================= c A slightly modified version of the 1D-minimization code from: c c "Computer Methods for Mathematical Computations" c by G.E. Forsythe, M.A. Malcolm and C.B. Moler c (Englewood Cliffs, NJ: Prentice-Hall, 1977) c c Retrieves the best estimate (with accuracy tol) of the point c within the interval (ax,bx) where function F(x) attains a c minimum c================================================================= implicit real*8 (a-h,o-z) dimension co(*), co1(*) c c is the squared inverse of the golden ratio c = 0.5D0*(3.D0 - DSQRT(5.0D0)) c eps is approximately the square root of the relative machine c precision. eps = 1.0D00 10 eps = eps/2.0D00 tol1 = 1.0D0 + eps if (tol1 .gt. 1.0D00) go to 10 eps = DSQRT(eps) c initialization a = ax b = bx v = a + c*(b - a) w = v x = v e = 0.0D0 fx = F(x,co,co1) fv = fx fw = fx c main loop starts here 20 xm = 0.5D0*(a + b) tol1 = eps*DABS(x) + tol/3.0D0 tol2 = 2.0D0*tol1 c check stopping criterion if (DABS(x - xm) .le. (tol2 - 0.5D0*(b - a))) go to 90 c is golden-section necessary if (DABS(e) .le. tol1) go to 40 c fit parabola r = (x - w)*(fx - fv) q = (x - v)*(fx - fw) p = (x - v)*q - (x - w)*r q = 2.0D00*(q - r) if (q .gt. 0.0D0) p = -p q = DABS(q) r = e e = d c is parabola acceptable 30 if (DABS(p) .ge. DABS(0.5D0*q*r)) go to 40 if (p .le. q*(a - x)) go to 40 if (p .ge. q*(b - x)) go to 40 c a parabolic interpolation step d = p/q u = x + d c f must not be evaluated too close to ax or bx if ((u - a) .lt. tol2) d = DSIGN(tol1, xm - x) if ((b - u) .lt. tol2) d = DSIGN(tol1, xm - x) go to 50 c a golden-section step 40 if (x .ge. xm) e = a - x if (x .lt. xm) e = b - x d = c*e c f must not be evaluated too close to x 50 if (DABS(d) .ge. tol1) u = x + d if (DABS(d) .lt. tol1) u = x + DSIGN(tol1, d) fu = F(u,co,co1) c update a, b, v, w, and x if (fu .gt. fx) go to 60 if (u .ge. x) a = x if (u .lt. x) b = x v = w fv = fw w = x fw = fx x = u fx = fu go to 20 60 if (u .lt. x) a = u if (u .ge. x) b = u if (fu .le. fw) go to 70 if (w .eq. x) go to 70 if (fu .le. fv) go to 80 if (v .eq. x) go to 80 if (v .eq. w) go to 80 go to 20 70 v = w fv = fw w = u fw = fu go to 20 80 v = u fv = fu go to 20 c end of main loop 90 fmin = x return end FUNCTION F(x,co,co1) c----------------------------------------- c Yuri Mishin, June 1998 c==================================================== c Calculation of the objective function f(x) for the c minimization routine FMIN. Depending on the c minimization mode minmode, the argument of this c function is either rigid-body displacement or c the uniform deformation factor c==================================================== implicit real*8 (a-h,o-z) common/minimize/ xdisp1,ydisp1,zdisp1,minmode,iter1D dimension co(*),co1(*) iter1D = iter1D +1 if (minmode.eq.2) then CALL DEFORM(x,x,x,0.D0,co,co1) F = ENERGY(co) return else if (minmode.eq.3) then CALL DEFORM(1.D0,x,1.D0,0.D0,co,co1) F = ENERGY(co) return else if (minmode.eq.4) then CALL DEFORM(1.D0,1.D0,1.D0,x,co,co1) F = ENERGY(co) return else write(6, 100) minmode stop endif 100 format(' *** WRONG MINIMIZATION MODE (minmode =',i2, &' IN FUNCTION F) ***') end SUBROUTINE TABLES (co,newlc) *----------------------------------------- * Yuri Mishin, June 1998 *===================================================================== * Develops tables of neighbours according to a given search radius: * - numnbr(i) is the number of neighbours of atom i. * - neighbrs(i,j) is the atomic # of j-th neighbour of atom # i. * - idisp(i,j,k) is displacement of j-th neighbour of atom # i. * - realn(i,j) is true if j-th neighbour is a real atom and false * if it is a periodic image of a real atom * The construction is based on a link cell map with lcx, lcy and * lcz cells in the respective directions. The cell sizes (dx, dy and * dz) are slightly larger than the search radius. The minimum * admissible number of cells in each direction is 3 for fixed * boundary conditions and 4 for periodic boundary conditions; * otherwise only one cell is used. c===================================================================== implicit real*8 (a-h,o-z) logical bcx,bcy,bcz,sx,sy,sz,disp,newlc,lcready, & real,realc,realn parameter (maxatom = 100000,! max number of atoms in block & maxnb = 150, ! max number of neighbours & maxlc = 1000, ! max number of link cells & maxnlc = 27 ! max number of neighboring cells &) common /accuracy/ acc common /block/ search2,xmn,ymn,zmn,xmx,ymx,zmx, & xtmn,ytmn,ztmn,xtmx,ytmx,ztmx, & xdisp,ydisp,zdisp,mdx,mdy,mdz, & ntot,nb,nfree common /cyc/ bcx,bcy,bcz common /table1/ neighbrs(maxatom,maxnb),idisp(maxatom,maxnb,3) common /table2/ numnbr(maxatom),realn(maxatom,maxnb) common /lc/ ntop(maxlc),link(maxatom),lcx,lcy,lcz, & neilc(maxlc,maxnlc),lcd(maxlc,maxnlc,3), & numnlc(maxlc),realc(maxlc,maxnlc) dimension co(*) data sx/.false./, sy/.false./, sz/.false./, lcready/.false./ dmin = DSQRT(search2) if (.not.newlc) go to 10 ***************************** * Prepare a link cell map * ***************************** *--- Divide the block into cells ---- xlc = xtmx - xtmn lcx = INT(xlc/dmin) if (bcx.and.lcx.le.3) then lcx = 1 sx = .true. endif if (.not.bcx.and.lcx.le.2) lcx = 1 ylc = ytmx - ytmn lcy = INT(ylc/dmin) if (bcy.and.lcy.le.3) then lcy = 1 sy = .true. endif if (.not.bcy.and.lcy.le.2) lcy = 1 zlc = ztmx - ztmn lcz = INT(zlc/dmin) if (bcz.and.lcz.le.3) then lcz = 1 sz = .true. endif if (.not.bcz.and.lcz.le.2) lcz = 1 dx = xlc/lcx ! cell size in x direction dy = ylc/lcy ! cell size in y direction dz = zlc/lcz ! cell size in z direction lcxy = lcx*lcy lctot = lcz*lcxy if (lctot.gt.maxlc) then write(6,103) lctot, maxlc stop endif write(6, 105) lcx, lcy, lcz *--- Construct the table of neighbours of link cells ----- do ix = 1, lcx do iy = 1, lcy do iz = 1, lcz mlc = lcxy*(iz-1) + lcx*(iy-1) + ix numnlc(mlc) = 0 do jx = -1,1 do jy = -1,1 do jz = -1,1 kx = ix + jx ky = iy + jy kz = iz + jz if((kx.lt.1.or.kx.gt.lcx).and. & (.not.bcx.or.lcx.eq.1)) go to 11 if((ky.lt.1.or.ky.gt.lcy).and. & (.not.bcy.or.lcy.eq.1)) go to 11 if((kz.lt.1.or.kz.gt.lcz).and. & (.not.bcz.or.lcz.eq.1)) go to 11 i1 = 0 i2 = 0 i3 = 0 if(lcx.gt.1.and.bcx) then if(kx.lt.1) then kx = lcx i1 = -1 endif if(kx.gt.lcx) then kx = 1 i1 = 1 endif endif if(lcy.gt.1.and.bcy) then if(ky.lt.1) then ky = lcy i2 = -1 endif if(ky.gt.lcy) then ky = 1 i2 = 1 endif endif if(lcz.gt.1.and.bcz) then if(kz.lt.1) then kz = lcz i3 = -1 endif if(kz.gt.lcz) then kz = 1 i3 = 1 endif endif klc = lcxy*(kz-1) + lcx*(ky-1) + kx numnlc(mlc) = numnlc(mlc) + 1 neilc(mlc,numnlc(mlc)) = klc lcd(mlc,numnlc(mlc),1) = i1 lcd(mlc,numnlc(mlc),2) = i2 lcd(mlc,numnlc(mlc),3) = i3 if (i1.ne.0.or.i2.ne.0.or.i3.ne.0) then realc(mlc,numnlc(mlc)) = .false. else realc(mlc,numnlc(mlc)) = .true. endif 11 enddo enddo enddo enddo enddo enddo *--- Distribute atoms among the link cells ------ do i=1,maxlc ntop(i) = 0 enddo do ia = 1, ntot i = 3*ia x1 = co(i-2) y1 = co(i-1) z1 = co(i) ix = INT((x1 - xtmn)/dx) + 1 iy = INT((y1 - ytmn)/dy) + 1 iz = INT((z1 - ztmn)/dz) + 1 mlc = lcxy*(iz-1) + lcx*(iy-1) + ix j = ntop(mlc) ntop(mlc) = ia link(ia) = j enddo lcready = .true. *==== The link cell map is ready ======= 10 continue if (.not.lcready) then write(6, 104) stop endif ********************************************************* * Construct the atomic table of neighbours using the * * link cell map * ********************************************************* write(6, 106) do n1 = 1, nb numnbr(n1) = 0 enddo do i = 1,nb do j=1, maxnb realn(i,j) = .false. enddo enddo maxnbr = 0 do ic = 1, lctot ! choose a cell (ic) do i = 1, numnlc(ic) ! look at all neighbouring cells (j) j = neilc(ic,i) * if necessary, apply additional translations of cell j * in periodic directions do ii = -mdx,mdx if (.not.sx.and.ii.ne.0) go to 16 do jj = -mdy,mdy if (.not.sy.and.jj.ne.0) go to 15 do kk = -mdz,mdz if (.not.sz.and.kk.ne.0) go to 14 if (ii.eq.0.and.jj.eq.0.and.kk.eq.0) then disp = .false. else disp = .true. endif real = .not.disp.and.realc(ic,i) * if (real), then cell j consists of real atoms, * otherwise it's a periodic image of a real cell n1 = ntop(ic) ! choose and atom in cell ic if(n1.eq.0) go to 21 ! skip cell if it's empty 12 continue if (n1.gt.nb) go to 18 k1 = 3*n1 j1 = k1-1 i1 = j1-1 x1 = co(i1) y1 = co(j1) z1 = co(k1) xsmin = x1 - dmin ! draw a box around atom n1 xsmax = x1 + dmin ysmin = y1 - dmin ysmax = y1 + dmin zsmin = z1 - dmin zsmax = z1 + dmin n2 = ntop(j) ! choose and atom in cell j if(n2.eq.0) go to 20 ! skip cell if it's empty 13 continue if (real.and.n2.le.n1) go to 19 k2 =3*n2 j2 = k2 - 1 i2 = j2 - 1 x2 = co(i2) y2 = co(j2) z2 = co(k2) * apply all necessary translations to atom n2 if (.not.realc(ic,i)) then x2 = x2 + xdisp*lcd(ic,i,1) y2 = y2 + ydisp*lcd(ic,i,2) z2 = z2 + zdisp*lcd(ic,i,3) endif if(disp) then if(ii.ne.0) x2 = x2 + xdisp*ii if(jj.ne.0) y2 = y2 + ydisp*jj if(kk.ne.0) z2 = z2 + zdisp*kk endif if ((x2.le.xsmax). and. ! is atom n2 in the box? & (x2.ge.xsmin). and. & (y2.le.ysmax). and. & (y2.ge.ysmin). and. & (z2.le.zsmax). and. & (z2.ge.zsmin)) then x12 = x2 - x1 y12 = y2 - y1 z12 = z2 - z1 r12 = x12**2 + y12**2 + z12**2 * is atom n2 within the search sphere? if (r12.le.search2.and.r12.gt.acc) then * add atom n2 to the neighbour list numnbr(n1) = numnbr(n1) + 1 nbr = numnbr(n1) if(nbr.gt.maxnb) then write(6, 101) maxnb stop endif neighbrs(n1,nbr)= n2 idisp(n1,nbr,1) = lcd(ic,i,1) idisp(n1,nbr,2) = lcd(ic,i,2) idisp(n1,nbr,3) = lcd(ic,i,3) realn(n1,nbr) = real if(.not.disp) go to 17 if(ii.ne.0) & idisp(n1,nbr,1) = idisp(n1,nbr,1) + ii if(jj.ne.0) & idisp(n1,nbr,2) = idisp(n1,nbr,2) + jj if(kk.ne.0) & idisp(n1,nbr,3) = idisp(n1,nbr,3) + kk 17 continue endif endif 19 n2 = link(n2) if(n2.gt.0) go to 13 if (nbr.gt.maxnbr) maxnbr = nbr 18 n1 = link(n1) if (n1.gt.0) go to 12 14 enddo 15 enddo 16 enddo 20 enddo 21 enddo write(6, 102) maxnbr 101 format(/' NUMBER OF NEIGHBOURS >',I6,' !!!', &/' INCREASE 2ND DIMENSION OF ARRAYS neighbrs and idisp.'/) 102 format(' MAXIMUM NUMBER OF NEIGHBOURS IS ',i3/) 103 format(//' TOO MANY LINK CELLS:',I4,' >',I4,' !!!'//) 104 format(//' THE LINK CELL MAP IS NOT READY !!!'//) 105 format(/' THE LINK CELL MAP IS', i3,' by',I3,' by',i3) 106 format(/' CONSTRUCTING TABLE OF NEIGHBOURS') RETURN END SUBROUTINE PERFECT (co,E) *----------------------------------------- * Yuri Mishin, June 1998 *================================================================ * This simple subroutine calculates the cohesive energy and the * electron densities on atoms in the perfect lattice *================================================================ implicit real*8 (a-h,o-z) parameter (nspmax = 4, ! max number of chem. species & nbasmax = 30000, ! max number of basis atoms & maxatom = 100000 ! max number of atoms in block &) logical bcx,bcy,bcz character*4 element(nspmax) common /accuracy/ acc common /chemical/ cutoff2,np(nspmax,nspmax),numat,npairs,npr, & element common /basis/ xbas(nbasmax),ybas(nbasmax),zbas(nbasmax), & lbas(nbasmax),nbas common /block/ search2,xmn,ymn,zmn,xmx,ymx,zmx, & xtmn,ytmn,ztmn,xtmx,ytmx,ztmx, & xdisp,ydisp,zdisp,mdx,mdy,mdz, & ntot,nb,nfree common /atoms/ lco(maxatom) common /cyc/ bcx,bcy,bcz dimension co(*), rho(nbasmax) E=0.D0 do i=1,nbas rho(i) = 0.D0 enddo *--- Choose a basis atom ------- do n1=1,nbas k1 = 3*n1 j1 = k1-1 i1 = j1 -1 z1 = co(k1) y1 = co(j1) x1 = co(i1) *--- Choose a neighbour ------- do n2=1,ntot k2 = 3*n2 j2 = k2-1 i2 = j2-1 x2 = co(i2) y2 = co(j2) z2 = co(k2) kp = np(lco(n1),lco(n2)) kf2 = npairs + lco(n2) * Try possible displacement of atom n2 (for cyclic boundaries) do ii = -mdx,mdx if(.not.bcx.and.ii.ne.0) go to 3 do jj = -mdy,mdy if(.not.bcy.and.jj.ne.0) go to 2 do kk = -mdz,mdz if(.not.bcz.and.kk.ne.0) go to 1 *--- Determine interatomic distance between atoms n1 and n2 --- u = (x2 + xdisp*ii - x1)**2 + & (y2 + ydisp*jj - y1)**2 + & (z2 + zdisp*kk - z1)**2 if (u.lt.acc) go to 1 u = DSQRT(u) *--- Get the pair interaction energy and electron density ---- E = E + 0.5D0*SPL(kp,u,0) rho(n1) = rho(n1) + SPL(kf2,u,0) 1 enddo 2 enddo 3 enddo enddo enddo * The pair interaction energy and electron densities * on the basis atoms are known. Now add the embedding energy * and get the total energy. do n1=1, nbas k = npr + lco(n1) E = E + SPL(k,rho(n1),0) enddo E = E/nbas write(6, 100) E write(6, 101) do i=1,nbas write(6, 102) i, element(lbas(i)),rho(i) enddo write(6, 103) 100 format(' COHESIVE ENERGY OF THE PERFECT LATTICE:',f10.6/) 101 format(' ELECTRON DENSITY IN THE PERFECT LATTICE:', & /' LCO ELEMENT EL. DENSITY') 102 format(i9,6x,a4,4x,f10.6) 103 format(/' --------------------------------------------------') RETURN END FUNCTION ENERGY (co) *----------------------------------------- * Yuri Mishin, June 1998 *==================================================================== * Calculation of the total energy of an extended simulation block. * An extended block includes both free and buffer atoms. * Because this calculation uses tables of neighbours, subroutine c TABLES must be called before calling ENERGY. *==================================================================== implicit real*8 (a-h,o-z) logical realn parameter (nspmax = 4, ! max number of chem. species & maxatom = 100000,! max number of atoms in block & maxnb = 150 ! max number of neighbours &) character*4 element(nspmax) common /block/ search2,xmn,ymn,zmn,xmx,ymx,zmx, & xtmn,ytmn,ztmn,xtmx,ytmx,ztmx, & xdisp,ydisp,zdisp,mdx,mdy,mdz, & ntot,nb,nfree common /table1/ neighbrs(maxatom,maxnb),idisp(maxatom,maxnb,3) common /table2/ numnbr(maxatom),realn(maxatom,maxnb) common /chemical/ cutoff2,np(nspmax,nspmax),numat,npairs,npr, & element common /atoms/ lco(maxatom) dimension co(*), rho(maxatom) * Initialization E=0.D0 do i=1,nb rho(i) = 0.D0 enddo ********************************** * Choose a free of buffer atom * ********************************** do n1=1,nb if (lco(n1).eq.0) go to 1 ! skip it if it is a vacancy k3 = 3*n1 j3 = k3-1 i3 = j3 -1 z1 = co(k3) y1 = co(j3) x1 = co(i3) kf1 = npairs + lco(n1) ************************ * Choose a neighbour * ************************ do nbr = 1,numnbr(n1) n2 = neighbrs(n1,nbr) if (lco(n2).eq.0) go to 2 ! skip it if it is a vacancy kp = np(lco(n1),lco(n2)) kf2 = npairs + lco(n2) k4 = 3*n2 j4 = k4 -1 i4 = j4 -1 z2 = co(k4) y2 = co(j4) x2 = co(i4) if (.not.realn(n1,nbr)) then z2 = co(k4) + zdisp*idisp(n1,nbr,3) y2 = co(j4) + ydisp*idisp(n1,nbr,2) x2 = co(i4) + xdisp*idisp(n1,nbr,1) endif * --- Calculate the interatomic distance x12 = x1 - x2 y12 = y1 - y2 z12 = z1 - z2 u = x12*x12 + y12*y12 + z12*z12 * --- Compare u with the cut-off radius if(u.gt.cutoff2) go to 2 u = DSQRT(u) * --- Get the pair interaction energy V V = SPL(kp,u,0) if (.not.realn(n1,nbr).or.n2.gt.nb) V = 0.5D0*V E = E + V * --- Get the electron density on atom n1 rho(n1) = rho(n1) + SPL(kf2,u,0) * --- if n2 is a free or a buffer atom, then n1 makes a contribution * to the electron density on it if(realn(n1,nbr).and.n2.le.nb) & rho(n2) = rho(n2) + SPL(kf1,u,0) 2 enddo 1 enddo * The pair interaction energy and electron densities * on all free and buffer atoms are known. Now add the embedding * energy and get the total energy. do n1=1, nb if (lco(n1).eq.0) go to 3 ! skip n1 if it is a vacancy k = npr + lco(n1) E = E + SPL(k,rho(n1),0) 3 enddo ENERGY = E RETURN END SUBROUTINE ENFORCE (co,E,g) *----------------------------------------- * Yuri Mishin, June 1998 *==================================================================== * Calculation of the 'extended' energy of the simulation block and * the forces acting on the atoms. The 'extended' energy equals the * energy of the free region plus the buffer zone. The forces * on the free atoms are derivatives of the 'extended' energy with * respect to the atomic coordinates. Except when all boundaries are * cyclic, the value of E returned by this routine does NOT coincide * with the total energy of the free region. * * Because this calculation uses tables of neighbours, subroutine c TABLES must be called before calling ENFORCE. *==================================================================== implicit real*8 (a-h,o-z) logical realn parameter (nspmax = 4, ! max number of chem. species & maxatom = 100000,! max number of atoms in block & maxnb = 150 ! max number of neighbours &) character*4 element(nspmax) common /block/ search2,xmn,ymn,zmn,xmx,ymx,zmx, & xtmn,ytmn,ztmn,xtmx,ytmx,ztmx, & xdisp,ydisp,zdisp,mdx,mdy,mdz, & ntot,nb,nfree common /table1/ neighbrs(maxatom,maxnb),idisp(maxatom,maxnb,3) common /table2/ numnbr(maxatom),realn(maxatom,maxnb) common /chemical/ cutoff2,np(nspmax,nspmax),numat,npairs,npr, & element common /atoms/ lco(maxatom) dimension co(*), g(*), rho(maxatom), Fi(maxatom) equivalence (rho,Fi) * Initialization E=0.D0 do i=1,nb rho(i) = 0.D0 if(i.le.nfree) then i3 = 3*i G(i3-2) = 0.D0 G(i3-1) = 0.D0 G(i3) = 0.D0 endif enddo *********************************** * Choose a free or buffer atom * *********************************** do n1=1,nb if (lco(n1).eq.0) go to 1 ! skip n1 if it is a vacancy k3 = 3*n1 j3 = k3-1 i3 = j3 -1 z1 = co(k3) y1 = co(j3) x1 = co(i3) kf1 = npairs + lco(n1) ************************ * Choose a neighbour * ************************ do nbr = 1,numnbr(n1) n2 = neighbrs(n1,nbr) if (lco(n2).eq.0) go to 2 ! skip n2 if it is a vacancy kp = np(lco(n1),lco(n2)) kf2 = npairs + lco(n2) k4 = 3*n2 j4 = k4 -1 i4 = j4 -1 z2 = co(k4) y2 = co(j4) x2 = co(i4) if (.not.realn(n1,nbr)) then z2 = co(k4) + zdisp*idisp(n1,nbr,3) y2 = co(j4) + ydisp*idisp(n1,nbr,2) x2 = co(i4) + xdisp*idisp(n1,nbr,1) endif * --- Calculate the interatomic distance x12 = x1 - x2 y12 = y1 - y2 z12 = z1 - z2 u = x12*x12 + y12*y12 + z12*z12 * --- Compare u with the cut-off radius if(u.gt.cutoff2) go to 2 u = DSQRT(u) * --- If n1 is a free or buffer atom, get the pair interaction energy * V and the pair interaction forces F V = SPL(kp,u,0) if (.not.realn(n1,nbr).or.n2.gt.nb) V = 0.5D0*V E = E + V if(n1.eq.n2) go to 6 ! no interaction with * self-images if(n1.le.nfree) then ff = SPL(kp,u,1)/u Fx = ff*x12 Fy = ff*y12 Fz = ff*z12 G(k3) = G(k3) + Fz G(j3) = G(j3) + Fy G(i3) = G(i3) + Fx * --- if n2 is a free atom, then n1 makes a contribution * to the force acting on it if (realn(n1,nbr).and.n2.le.nfree) then G(k4) = G(k4) - Fz G(j4) = G(j4) - Fy G(i4) = G(i4) - Fx endif 6 continue endif * --- Get the electron density on atom n1: rho(n1) = rho(n1) + SPL(kf2,u,0) * --- if n2 is a free or buffer atom, then n1 makes a contribution * to the electron density on it if(realn(n1,nbr).and.n2.le.nb) & rho(n2) = rho(n2) + SPL(kf1,u,0) 2 enddo 1 enddo * The pair interaction energy and electron densities on all free * and buffer atoms are known. Now add the embedding energy * and get the total energy. In the same loop we determine the * derivatives of the embedding function on all free and buffer atoms do n1=1, nb if (lco(n1).eq.0) go to 3 ! skip n1 if it is a vacancy k = npr + lco(n1) E = E + SPL(k,rho(n1),0) Fi(n1) = SPL(k,rho(n1),1) 3 enddo * ===== Start a new loop to calculate the embedding forces ===== ************************* * Choose a free atom * ************************* do n1=1,nfree if (lco(n1).eq.0) go to 4 ! skip n1 if it is a vacancy k3 = 3*n1 j3 = k3-1 i3 = j3 -1 z1 = co(k3) y1 = co(j3) x1 = co(i3) kf1 = npairs + lco(n1) ************************ * Choose a neighbour * ************************ do nbr = 1,numnbr(n1) n2 = neighbrs(n1,nbr) if (n2.eq.n1) go to 5 ! no interaction with * self-images if (lco(n2).eq.0) go to 5 ! skip n2 if it is a vacancy if (n2.gt.nb) go to 5 ! skip n2 if it is fixed kp = np(lco(n1),lco(n2)) kf2 = npairs + lco(n2) k4 = 3*n2 j4 = k4 -1 i4 = j4 -1 z2 = co(k4) y2 = co(j4) x2 = co(i4) if (.not.realn(n1,nbr)) then z2 = co(k4) + zdisp*idisp(n1,nbr,3) y2 = co(j4) + ydisp*idisp(n1,nbr,2) x2 = co(i4) + xdisp*idisp(n1,nbr,1) endif * --- Calculate the interatomic distance x12 = x1 - x2 y12 = y1 - y2 z12 = z1 - z2 u = x12*x12 + y12*y12 + z12*z12 * --- Compare u with the cut-off radius if(u.gt.cutoff2) go to 5 u = DSQRT(u) * --- Get the embedding force acting on atom n1 ff = Fi(n1)*SPL(kf2,u,1)/u & + Fi(n2)*SPL(kf1,u,1)/u Fx = ff*x12 Fy = ff*y12 Fz = ff*z12 G(k3) = G(k3) + Fz G(j3) = G(j3) + Fy G(i3) = G(i3) + Fx * --- if n2 is a free atom, then n1 makes a contribution * to the force acting on it if (realn(n1,nbr).and.n2.le.nfree) then G(k4) = G(k4) - Fz G(j4) = G(j4) - Fy G(i4) = G(i4) - Fx endif 5 enddo 4 enddo RETURN END SUBROUTINE CONJUG (n,x,f,g,ifun,iter,eps,nflag,mxfun,w, &iout,idev,acc) c=================================================================== c This subroutine minimizes a given function f by the conjugate c gradient method. c n - The number of variables in the function to be minimized c x - The vector containing the current estimate to the minimizer c On entry, x contains an initial estimate supplied by the c calling program. On exiting, x holds the best estimate to c the minimizer obtained c f - On exiting, f contains the lowest value of the object c function obtained c g - On exiting, g contains the elements of the gradient of f c evaluated at point x c ifun - On exiting,ifun contains the number of times the function c and gradient have been evaluated c iter - On exiting, contains the total number of search directions c calculated to obtain the current estimate to the minizer c eps - Convergence parameter. Convergence occurs when the norm of c the gradient is less than or equal to eps times the maximum c of one and the norm of the vector x c nflag- On exiting, indicates which condition caused the exit: c * if nflag=0, the algorithm has converged c * if nflag=1, the maximum number of function evaluations c have been used c * if nflag=2, the linear search has failed to improve the c function value. This is the usual exit if c either the function or the gradient is c incorrectly coded c * if nflag=3, the search vector was not a descent direction c This can only be caused by roundoff,and may c suggest that the convergence criterion is c too strict c mxfun- The maximum allowed number of function and gradient calls c w - A vector of working storage. Must be dimensioned 5*n+2 c iout - An output parameter. If iout = 0, there is no printed c output from this subroutine. If iout > 0, the value of f, c the norm of the gradient squared, iter and ifun are written c every iout iterations c idev - Number of the output device on which output is to be c written when iout>0 c acc - An estimate of machine accuracy. A linear search is c unsuccessfully terminated when the norm of the step size c becomes smaller than acc. In practice, acc=10.d-20 has c proved satisfactory. c c The function to be minimized is calculated in a subroutine c ENFORCE which calculates the function and gradient at c x and places them in f and g(1),...,g(n) respectively. c The subroutine must have the form: c c SUBROUTINE ENFORCE(x,f,g) c c The vector dimension n must be delivered to ENFORCE via a common c block. c c This is a modified version of the algorithm developed by c D.F. Shanno and K.H. Phua, c ACM Transactions on Mathematical Software 6 (December 1980), c 618-622. implicit real*8 (a-h,o-z) dimension x(*),g(*),w(*) logical rsw iter=0 ifun=0 nflag=0 c Set parameters to extract vectors from w. c w(i) holds the search vector,w(nx+i) holds the best current c estimate to the minimizer,and w(ng+i) holds the gradient c at the best current estimate. nx=n ng=nx+n c w(nry+i) holds the restart y vector and w(nrd+i) holds the c restart search vector. nry=ng+n nrd=nry+n ncons=5*n ncons1=ncons+1 ncons2=ncons+2 c Calculate the function and gradient at the initial c point and initialize nrst,which is used to determine c whether a beale restart is being done. nrst=n means that this c iteration is a restart iteration. initialize rsw,which indicates c that the current search direction is a gradient direction. 20 call ENFORCE(x,f,g) ifun=ifun+1 nrst=n rsw=.true. c Calculate the initial search direction , the norm of x squared, c and the norm of g squared. dg1 is the current directional c derivative,while xsq and gsq are the squared norms. dg1=0.D0 xsq=0.D0 do i=1,n w(i)=-g(i) xsq=xsq+x(i)*x(i) dg1=dg1-g(i)*g(i) enddo gsq=-dg1 c Test if the initial point is the minimizer. if(gsq.le.eps*eps*DMAX1(1.0D0,xsq)) return c Begin the major iteration loop. ncalls is used to guarantee that c at least two points have been tried. fmin is the c current function value. 40 fmin=f ncalls=ifun c If output is desired,test if this is the correct iteration c and if so, write output. if (iout.ne.0) then if(MOD(iter,iout).eq.0) then write(idev,500) iter,ifun,fmin,gsq endif endif c Begin linear search. alpha is the steplength. c set alpha to the nonrestart conjugate gradient alpha. alpha=alpha*dg/dg1 c If a restart has been performed, set alpha=1.0. if(nrst.eq.1) alpha=1.D0 c If a gradient direction is used, set alpha=1.0/DSQRT(gsq), c which scales the initial search vector to unity. if(rsw) alpha=1.D0/DSQRT(gsq) c The linear search fits a cubic to f and dal, the function and its c derivative at alpha, and to fp and dp,the function c and derivative at the previous trial point ap. c Initialize ap ,fp,and dp. ap=0.D0 fp=fmin dp=dg1 c Save the current derivative to scale the next search vector. dg=dg1 c Update the iteration. iter=iter+1 c Calculate the current steplength and store the current x and g. step=0.D0 do i=1,n step=step+w(i)*w(i) nxpi=nx+i ngpi=ng+i w(nxpi)=x(i) w(ngpi)=g(i) enddo step=DSQRT(step) c Begin the linear search iteration. c Test for failure of the linear search. c Test if direction is a gradient direction. 80 if(alpha*step.le.acc) then if(.not.rsw)go to 20 nflag=2 return endif c Calculate the trial point. do i=1,n nxpi=nx+i x(i)=w(nxpi)+alpha*w(i) enddo c Evaluate the function at the trial point. call ENFORCE(x,f,g) c Test if the maximum number of function calls have been used. ifun=ifun+1 if(ifun.gt.mxfun) then nflag=1 return endif c Compute the derivative of f at alpha. dal=0.D0 do i=1,n dal=dal+g(i)*w(i) enddo c Test whether the new point has a negative slope but a higher c function value than alpha=0. If this is the case,the search c has passed through a local max and is heading for a distant local c minimum. if(f.gt.fmin.and.dal.lt.0.) go to 160 c If not, test whether the steplength criteria have been met. if(f.gt.(fmin+.0001D0*alpha*dg).or.DABS(dal/dg) &.gt.(0.9)) go to 130 c If they have been met, test if two points have been tried c and if the true line minimum has not been found. if((ifun-ncalls).le.1.and.DABS(dal/dg).gt.eps) go to 130 go to 170 c A new point must be tried. Use cubic interpolation to find c the trial point at. 130 u1=dp+dal-3.D0*(fp-f)/(ap-alpha) u2=u1*u1-dp*dal if(u2.lt.0.D0) u2=0.D0 u2=DSQRT(u2) at=alpha-(alpha-ap)*(dal+u2-u1)/(dal-dp+2.D0*u2) c Test whether the line minimum has been bracketed. if((dal/dp).gt.0.D0) go to 140 c The minimum has been bracketed. Test whether the trial point lies c sufficiently within the bracketed interval. c If it does not, choose at as the midpoint of the interval. if(at.lt.(1.01D0*DMIN1(alpha,ap)).or.at.gt.(.99D0*DMAX1 &(alpha,ap))) at=(alpha+ap)/2.D0 go to 150 c The minimum has not been bracketed. Test if both points are c greater than the minimum and the trial point is sufficiently c smaller than either. 140 if(dal .gt.0.D0.and.0.D0.lt.at.and.at.lt. & (.99D0*DMIN1(ap,alpha))) go to 150 c Test if both points are less than the minimum and the trial point c is sufficiently large. if(dal.le.0.D0.and.at.gt.(1.01D0*DMAX1(ap,alpha))) go to 150 c If the trial point is too small,double the largest prior point. if(dal.le.0.D0) at=2.D0*DMAX1(ap,alpha) c If the trial point is too large, halve the smallest prior point. if(dal.gt.0.D0) at=DMIN1(ap,alpha)/2.D0 c Set ap=alpha, alpha=at,and continue search. 150 ap=alpha fp=f dp=dal alpha=at go to 80 c A relative max has been passed. c Reduce alpha and restart the search. 160 alpha=alpha/3.D0 ap=0.D0 fp=fmin dp=dg go to 80 c The line search has converged. c Test for convergence of the algorithm. 170 gsq=0.D0 xsq=0.D0 do i=1,n gsq=gsq+g(i)*g(i) xsq=xsq+x(i)*x(i) enddo if(gsq.le.eps*eps*DMAX1(1.D0,xsq)) return c Search continues. set w(i)=alpha*w(i),the full step vector. do i=1,n w(i)=alpha*w(i) enddo c Compute the new search vector. first test whether a c conjugate gradient or a variable metric vector is used. c Conjugate gradient update section. c Test if a powell restart is indicated. rtst=0.D0 do i=1,n ngpi=ng+i rtst=rtst+g(i)*w(ngpi) enddo if(DABS(rtst/gsq).gt.0.2D0) nrst=n c If a restart is indicated, save the current d and y c as the beale restart vectors and save d'y and y'y c in w(ncons+1) and w(ncons+2). if(nrst.ne.n) go to 220 w(ncons+1)=0.D0 w(ncons+2)=0.D0 do i=1,n nrdpi=nrd+i nrypi=nry+i ngpi=ng+i w(nrypi)=g(i)-w(ngpi) w(nrdpi)=w(i) w(ncons1)=w(ncons1)+w(nrypi)*w(nrypi) w(ncons2)=w(ncons2)+w(i)*w(nrypi) enddo c Calculate the restart hessian times the current gradient. 220 u1=0.D0 u2=0.D0 do i=1,n nrdpi=nrd+i nrypi=nry+i u1=u1-w(nrdpi)*g(i)/w(ncons1) u2=u2+w(nrdpi)*g(i)*2.D0/w(ncons2)-w(nrypi)*g(i)/w(ncons1) enddo u3=w(ncons2)/w(ncons1) do i=1,n nxpi=nx+i nrdpi=nrd+i nrypi=nry+i w(nxpi)=-u3*g(i)-u1*w(nrypi)-u2*w(nrdpi) enddo c If this is a restart iteration,w(nx+i) contains the new search c vector. if(nrst.eq.n) go to 300 c Not a restart iteration. Calculate the restart hessian c times the current y. 250 u1=0.D0 u2=0.D0 u3=0.D0 u4=0.D0 do i=1,n ngpi=ng+i nrdpi=nrd+i nrypi=nry+i u1=u1-(g(i)-w(ngpi))*w(nrdpi)/w(ncons1) u2=u2-(g(i)-w(ngpi))*w(nrypi)/w(ncons1) & +2.D0*w(nrdpi)*(g(i)-w(ngpi))/w(ncons2) u3=u3+w(i)*(g(i)-w(ngpi)) enddo step=0.D0 do i=1,n ngpi=ng+i nrdpi=nrd+i nrypi=nry+i step=(w(ncons2)/w(ncons1))*(g(i)-w(ngpi)) & +u1*w(nrypi)+u2*w(nrdpi) u4=u4+step*(g(i)-w(ngpi)) w(ngpi)=step enddo c Calculate the doubly updated hessian times the current c gradient to obtain the search vector. u1=0.D0 u2=0.D0 do i=1,n u1=u1-w(i)*g(i)/u3 ngpi=ng+i u2=u2+(1.D0+u4/u3)*w(i)*g(i)/u3-w(ngpi)*g(i)/u3 enddo do i=1,n ngpi=ng+i nxpi=nx+i w(nxpi)=w(nxpi)-u1*w(ngpi)-u2*w(i) enddo c Calculate the derivative along the new search vector. 300 dg1=0.D0 do i=1,n nxpi=nx+i w(i)=w(nxpi) dg1=dg1+w(i)*g(i) enddo c If the new direction is not a descent direction, stop. if (dg1.gt.0.D0) go to 320 c Update nrst to assure at least one restart every n iterations. if(nrst.eq.n) nrst=0 nrst = nrst+1 rsw=.false. go to 40 c Roundoff has produced a bad direction. 320 nflag=3 return 500 format(' ITERATION',i4,' NUMBER OF FUNCTION CALLS:',i5, &/' E =', d15.8,' G-SQUARED =',d15.8/) end SUBROUTINE FACTOR(xx,yy,z) *============================================================== * This subroutine determines the largest common real factor z * for two given real numbers xx and yy, * i.e. it represents xx and yy as * * xx=nz and yy=mz * * where n and m are two smallest integers. * * This is a slightly modified version of routine FACTOR * developed by M. J. Norgett in 1972 *=============================================================== implicit real*8 (a-h,o-z) common/accuracy/ acc x = DABS(xx) y = DABS(yy) if ((x*y).lt.acc) then z = DMAX1(x,y) return else if (x.lt.y) then 1 i = INT(y/x) if (i.eq.0) then write(6,100) return else y = y - i * x if (y.lt.acc) then z=x return else 2 i = INT(x/y) if (i.eq.0) then write(6,100) return else x = x - i * y if (x.lt.acc) then z=y return else go to 1 endif endif endif endif else if (x.gt.y) then go to 2 else z=x return endif endif endif 100 format(//6x,'*** Error in FACTOR routine !!! ***') end SUBROUTINE INTERPOL *----------------------------------------- * Yuri Mishin, Feb. 1996; revised in June 1998 *=================================================================== * This subroutine reads in data from potential files and develops * cubic spline coefficients that are subsequently used by function * SPL to retrieve interpolated values of interatomic potentials and * their derivatives if requested. The input data and the spline * coefficients are stored in arrays a, b, c, and d in common * block /pot/. * * Subroutines used: SPLINE * *=================================================================== implicit real*8 (a-h,o-z) parameter (nspmax = 4, ! max number of chem. species & maxpnt = 3000, ! max number of tabulated values & maxf = nspmax*(nspmax+1)/2 + 2*nspmax & ! max number of potential files &) character*4 element(nspmax) character*72 filename(maxf) common /chemical/ cutoff2,np(nspmax,nspmax),numat,npairs,npr, & element common/pot/a(maxf,maxpnt),b(maxf,maxpnt),c(maxf,maxpnt), & d(maxf,maxpnt),x0(maxf),xn(maxf), & step(maxf),n(maxf) & dimension x(maxpnt),y(maxpnt),w1(maxpnt),w2(maxpnt),w3(maxpnt) * ==== Input the number of chemical species ========== read(1,*) numat if ((numat.lt.1).or.(numat.gt.4)) then write(6,100) stop endif *********************************************************** * Input names of potential files. * * The first numat*(numat+1)/2 files represent pair * * interaction (e.g. for numat = 3, p11, p12, p13, p22, * * p23, p33). These are followed by numat electron-density * * files and numat embedding-energy files. * *********************************************************** npairs = numat*(numat+1)/2 npr=npairs + numat nfiles = npr + numat write(6,103) do i=1,numat read(1,*) element(i) write(6, 104) i, element(i) enddo write(6,101) do i=1,nfiles read (1,*) filename(i) write(6,*) filename(i) enddo *--- Array np(i,j) stores numberss of pair interaction files *--- for species i and j nf = 0 do i=1,numat do j=i,numat nf = nf + 1 np(i,j) = nf np(j,i) = nf enddo enddo ******************************************************** * Read in data from the potential files and calculate * * the spline coefficients * ******************************************************** do i=1, nfiles nunit=10+i open (nunit,file=filename(i),status='old') read(nunit,*) n(i), x0(i), xn(i), step(i), idumm1, idumm2 read(nunit,*) (y(j), j=1,n(i)) step(i)=1.0D0/step(i) do j=1,n(i) x(j)=x0(i)+step(i)*j enddo call SPLINE(n(i),x,y,w1,w2,w3) do j=1,n(i) a(i,j)=y(j) b(i,j)=w1(j) c(i,j)=w2(j) d(i,j)=w3(j) enddo close(nunit) enddo ***************************************************** *Calculate the cut-off radius of atomic interaction * ***************************************************** cutoff=0.D0 do i=1, npr if (xn(i).gt.cutoff) cutoff=xn(i) enddo cutoff2=cutoff**2 write(6, 102) cutoff 100 format(//' WRONG NUMBER OF CHEMICAL SPECIES !!!'//) 101 format(/' THE FOLLOWING POTENTIAL FILES ARE USED:'/) 102 format(/' CUT-OFF RADIUS OF ATOMIC INTERACTION: ',f8.4) 103 format(/' CHEMICAL ELEMENTS:',' TYPE ELEMENT'/) 104 format(20x,i5,6x,a4) RETURN END SUBROUTINE SPLINE (n, x, y, b, c, d) **************************************************************** * Spline interpolation subroutine from * "Computer Methods for Mathematical Computations" * by G.E. Forsythe, M.A. Malcolm and C.B. Moler * (Englewood Cliffs, NJ: Prentice-Hall, 1977) **************************************************************** implicit real*8 (a-h,o-z) dimension x(*), y(*), b(*), c(*), d(*) * The coefficients b(i), c(i), and d(i), i=1,2,...,n are computed * for a cubic interpolating spline * s(x) = y(i) + b(i)*(x-x(i)) + c(i)*(x-x(i))**2 + d(i)*(x-x(i))**3 * for x(i) .le. x .le. x(i+1) * input.. * n = the number of data points or knots (n.ge.2) * x = the abscissas of the knots in strictly increasing order * y = the ordinates of the knots * output.. * b, c, d = arrays of spline coefficients as defined above. * using p to denote differentiation, * y(i) = s(x(i)) * b(i) = sp(x(i)) * c(i) = spp(x(i))/2 * d(i) = sppp(x(i))/6 (derivative from the right) * the accompanying function subprogram seval can be used * to evaluate the spline. nm1 = n-1 if ( n .lt. 2 ) return if ( n .lt. 3 ) go to 50 * set up tridiagonal system * b = diagonal, d = offdiagonal, c = right hand side. d(1) = x(2) - x(1) c(2) = (y(2) - y(1))/d(1) do 10 i = 2, nm1 d(i) = x(i+1) - x(i) b(i) = 2*(d(i-1) + d(i)) c(i+1) = (y(i+1) - y(i))/d(i) c(i) = c(i+1) - c(i) 10 continue * end conditions. third derivatives at x(1) and x(n) * obtained from divided differences b(1) = -d(1) b(n) = -d(n-1) c(1) = 0.D0 c(n) = 0.D0 if ( n .eq. 3 ) go to 15 c(1) = c(3)/(x(4)-x(2)) - c(2)/(x(3)-x(1)) c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3)) c(1) = c(1)*d(1)**2/(x(4)-x(1)) c(n) = -c(n)*d(n-1)**2/(x(n)-x(n-3)) * forward elimination 15 do 20 i = 2, n t = d(i-1)/b(i-1) b(i) = b(i) - t*d(i-1) c(i) = c(i) - t*c(i-1) 20 continue * back substitution c(n) = c(n)/b(n) do 30 ib = 1, nm1 i = n-ib c(i) = (c(i) - d(i)*c(i+1))/b(i) 30 continue * c(i) is now the sigma(i) of the text * compute polynomial coefficients b(n) = (y(n) - y(nm1))/d(nm1) + d(nm1)*(c(nm1) + 2*c(n)) do 40 i = 1, nm1 b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2*c(i)) d(i) = (c(i+1) - c(i))/d(i) c(i) = 3*c(i) 40 continue c(n) = 3*c(n) d(n) = d(n-1) return 50 b(1) = (y(2)-y(1))/(x(2)-x(1)) c(1) = 0.D0 d(1) = 0.D0 b(2) = b(1) c(2) = 0.D0 d(2) = 0.D0 RETURN END FUNCTION SPL(k,u,m) *----------------------------------------- * Yuri Mishin, Feb. 1996; revised in June 1998 *=================================================================== * Calculates the cubic spline function * spl = y(i) + b(i)*(u-x(i)) + c(i)*(u-x(i))**2 + d(i)*(u-x(i))**3 * and its 1st, 2nd and 3rd derivatives for an arbitrary argument u, * where x(i) .lt. u .lt. x(i+1). * If u .lt. x(1) then i = 1 is used. * If u .ge. x(n) then i = n is used for the embedding function * while spl = 0 is set for pair interaction * and electronic density * x and y are arrays of abscissas and ordinates, b, c and d are * arrays of spline coefficients computed by subroutine SPLINE. * Input parameter m indicates that the m-th derivative (1<= m <=3) * of spl(u) or only spl(u) (m = 0) must be returned. * Up to maxf data sets, enumerated by index k, can be used in this * function. * ================================================================= implicit real*8 (a-h,o-z) parameter (nspmax = 4, ! max number of chem. species & maxpnt = 3000, ! max number of tabulated values & maxf = nspmax*(nspmax+1)/2 + 2*nspmax & ! max number of potential files &) character*4 element(nspmax) common /chemical/ cutoff2,np(nspmax,nspmax),numat,npairs,npr, & element common/pot/a(maxf,maxpnt),b(maxf,maxpnt),c(maxf,maxpnt), & d(maxf,maxpnt),x0(maxf),xn(maxf), & step(maxf),n(maxf) if(u.gt.xn(k)) then if(k.le.npr) then spl = 0.D0 return endif endif *------- Find i such that x(i) <= u < x(i+1) ux=u-x0(k) s = step(k) if (ux.lt.s) then i=1 else i=INT(ux/s) if (i.gt.n(k)) i=n(k) endif *------------ calculate spl(u) or its derivatives dx = ux - i*s if (m.eq.0) then spl = a(k,i) + dx*(b(k,i) + dx*(c(k,i) + dx*d(k,i))) else if(m.eq.1) then spl = b(k,i) + dx*(2*c(k,i) + 3*dx*d(k,i)) else if (m.eq.2) then spl = 2*c(k,i) + 6*dx*d(k,i) else if (m.eq.3) then spl = 6*d(k,i) else write(6,100) m stop endif endif endif endif 100 format(//' WRONG INPUT INDEX m = ', i3, &' IN FUNCTION SPL !!!'//) RETURN END