program NMD C C PERFORMS MOLECULAR DYNAMICS ON CHARGED PARTICLES C INTERACTING WITH EMBEDDED ATOM LIKE POTENTIALS. CHARGE C TRANSFER IS SOLVED EXPLICITELY USING EWALD SUMMATION IN C EITHER 2 OR 3 DIMENSIONS. C implicit double precision (a-h,o-z) include "PARAMETERS.h" include "units.inc" parameter (SQRAMU=4.074960123D01) double precision kcell(3,3),Omega dimension rnuc(3,MAXATM),vnuc(3,MAXATM) dimension rold(3,MAXATM) dimension h0(3,3),h1(3,3) dimension papp(3,3),amass(MAXATM) dimension force(3,MAXATM) dimension rrvec(0:3,MAXVEC) dimension rkvec(0:3,MAXVEC) dimension rxvec(0:3,MAXVEC) dimension indx(MAXIND) dimension pden(0:3,MAXATM) dimension ia(MAXATM) dimension zchg(MAXATM) dimension fvold(3,MAXATM) dimension vewald(MXHERM) dimension vwork(MXHERM) dimension catom(MAXATM),work(MAXATM,7) dimension xvec(6,MAXATM),pvec(6,MAXATM) dimension dqdeps(6,MAXATM) dimension prs(6) dimension strain(6) dimension nnuc(MAXTYP) dimension hard(MAXTYP) dimension chi(MAXTYP) dimension zeff(MAXTYP) dimension zeta(MAXTYP) dimension potparm(10,3) dimension nbrlist(0:2,MAXNBR) dimension xnbr(0:3,MAXNBR) dimension bot(3),top(3) dimension rmove(3,MAXATM) dimension hv(3) character*2 atsym(MAXTYP) character*9 file character*3 suf logical EQUIL,FINISHCONF,FINISHSTRAIN,QUP time = 0.d0 call GetInpV(NDIM,MAXATM,MAXTYP,natom,h0,kcell,Omega,rnuc, & zchg,vnuc,amass,nttot,ia,nnuc,hard,chi,zeff,zeta,potparm,time) do k=1,nttot call GtAtSy(atsym(k),nnuc(k)) enddo call GetPar(alphew,rcut,rmax,ecut,NDIM,Omega,h0,kcell) call MkLatt(NDIM,nrtot,rrvec,h0,kcell,rcut,rxvec,indx,MAXIND, & MAXVEC) print *,nrtot,' direct lattice vectors found' call MkLatt(NDIM,nktot,rkvec,kcell,h0,ecut,rxvec,indx,MAXIND, & MAXVEC) print *,nktot,' reciprocal lattice vectors found' c c READ STUFF FROM CONFIG FILE... c jconf=1 call GetConfig(jconf,nstep,tprint,itemp,jtest,istyle, & pmass,nconf,ilist,iconf,dt,teq,gamma,nbot,ntop, & dstrain,nstrain,pq,papp,hv,hf) c c INITIALIZE NEIGHBOR LIST ... c call makenbr(MAXNBR,NDIM,natom,h0,rnuc,nrtot,rrvec,rmax,nbrtot, & nbrlist,work) call update(natom,nrtot,nbrtot,nbrlist,rnuc,rrvec,xnbr) c c ... FORCES AND ENERGIES call aclear(3*natom,force) call aclear(6,prs) c call espotl(NDIM,natom,nrtot,nktot,nbrtot,nbrlist,xnbr,Omega, & rrvec,rkvec,rnuc,zchg,ia,hard,chi,zeff,zeta,potparm,alphew,ees, & eew,eembd,epair,etot,vewald,catom,force,pden,vwork,work,prs) call kinet(natom,nttot,vnuc,ia,prs,amass,rk,boxk) temp = EtoK*2.0*rk/float(3*natom-3) energy = (etot+rk)*EtoeV potential = etot*EtoeV rk = rk*EtoeV call writcn(NDIM,nttot,atsym,rnuc,vnuc,h0,natom,ia,zchg,time, & 'indata') call writxyz(natom,nttot,rnuc,vnuc,nnuc,ia,zchg,time,'xyz.ini') write (1,10000) time,energy,potential,rk,pressure,temp,zchg(1) if(istyle.eq.2.OR.istyle.eq.4) then call zeroforce(natom,nbot,ntop,force,bot,top) open (unit=17,file='zpull',status='unknown') write (17,11000) 0.d0,energy,bot,top close (17) endif close (1) open (unit=20,file='zhmat',status='unknown') open (unit=21,file='zpress',status='unknown') n = 0 if (jtest.eq.0) then isold = 0 isold2 = 0 else isold = 100000000 isold2 = 100000000 endif kntstrain = 0 kntequil = 0 if(istyle.eq.2) then jconf = 0 else jconf = 1 endif dist = 0.d0 length = 6 rtol = 1.0d-9 c c START OF AN MD STEP.... c do istep=1,nstep time = time + dt c c CHECK FOR EQUILIBRATION c if (jtest.ne.0) then call tcheck(1,natom,energy,teq,LENGTH,RTOL,EQUIL) endif c c IF EQUILIBRATED, OR OUT OF TIME, GET NEXT PRESSURE,TEMP, ETC. c FINISHCONF = (istep.eq.isold+iconf) FINISHSTRAIN = (istep.eq.isold2+nstrain) if (EQUIL.OR.FINISHCONF.OR.FINISHSTRAIN) then call writcn(NDIM,nttot,atsym,rnuc,vnuc,h0,natom,ia, & zchg,time,'restart') kntequil = kntequil + 1 write (suf,12000) kntequil file = 'xyz'//suf call writxyz(natom,nttot,rnuc,vnuc,nnuc,ia,zchg,time,file) open (unit=1,file='zoutput',status='old',access='append') write (1,10000) time,energy,potential,rk,pressure,temp, & zchg(1) close (unit=1) write (20,10000) time,h0 write (21,13000) time,prs if (EQUIL) then print *,'equilibrated at: ',istep if (istyle.EQ.1) then goto 100 else if (istyle.EQ.2) then FINISHSTRAIN = .TRUE. if (kntstrain.eq.nstrain) FINISHCONF = .TRUE. else if (istyle.EQ.3.OR.istyle.EQ.4) then FINISHCONF = .TRUE. endif else print *,'moving on at: ',istep endif print *,istep,FINISHSTRAIN,FINISHCONF,isold,isold2, & kntstrain,jconf c c GET NEXT CONFIG c if (FINISHCONF) then jconf = jconf + 1 if (jconf.gt.nconf) goto 100 if (jtest.eq.0) then isold2 = istep isold = istep else kntstrain = 0 endif call GetConfig(jconf,nstep,tprint,itemp,jtest,istyle, & pmass,nconf,ilist,iconf,dt,teq,gamma,nbot,ntop, & dstrain,nstrain,pq,papp,hv,hf) endif c c APPLY NEXT STRAIN c if (FINISHSTRAIN) then kntstrain = kntstrain + 1 call strainedge(natom,nbot,ntop,rnuc,dstrain) call makenbr(MAXNBR,NDIM,natom,h0,rnuc,nrtot,rrvec, & rmax,nbrtot,nbrlist,work) call update(natom,nrtot,nbrtot,nbrlist,rnuc,rrvec, & xnbr) if (jconf.eq.0) then call zeroforce(natom,nbot,ntop,force,bot,top) jconf = 1 endif open (unit=17,file='zpull',status='unknown', & access='append') write (17,11000) dist,energy,bot,top close (17) dist = dist + dstrain if (jtest.eq.0) isold2 = istep endif endif c c START VERLET STEP ... c call ACopy(3*natom,rnuc,1,rold) oldpot = potential call ACopy(3*natom,force,1,fvold) call verlet(natom,rnuc,vnuc,force,ia,amass,dt,1,rmove,qup) c c REMAKE NEIGHBOR LIST IF NECESSARY... c if(QUP) then call AClear(3*natom,rmove) call makenbr(MAXNBR,NDIM,natom,h0,rnuc,nrtot,rrvec,rmax, & nbrtot,nbrlist,work) endif c c CALCULATE NEW FORCES, ENERGIES c call update(natom,nrtot,nbrtot,nbrlist,rnuc,rrvec,xnbr) call espotl(NDIM,natom,nrtot,nktot,nbrtot,nbrlist,xnbr,Omega, & rrvec,rkvec,rnuc,zchg,ia,hard,chi,zeff,zeta,potparm,alphew, & ees,eew,eembd,epair,etot,vewald,catom,force,pden,vwork, & work) potential = etot*EtoeV if (istyle.eq.2.OR.istyle.eq.4.and.jconf.ne.0) then call zeroforce(natom,nbot,ntop,force,bot,top) endif if (istyle.eq.4) then call assignforce(natom,nbot,ntop,force,halfzforce,sum1,sum2) endif c c SOLVE FOR ELASTIC STRESSES AND MODULI, c REMAKE NEIGHBOR LIST c CURRENTLY DISABLED !! c if (istep.eq.0) then nptot = nttot*(nttot+1)/2 if (NDIM.eq.3) then call esnewton(NDIM,natom,nherm,strain,nbrtot,nbrlist, & xnbr,nrtot,rrvec,nktot,rkvec,h0,kcell,Omega,rnuc, & zchg,nttot,nptot,potparm,ia,zeff,zeta,alphew,pvec, & xvec,dqdeps,vewald,catom,pden,work) open (unit=19,file='zstrain',access='append', & status='unknown') write (19,14000) time,strain close (19) endif call makenbr(MAXNBR,NDIM,natom,h0,rnuc,nrtot,rrvec,rmax, & nbrtot,nbrlist,work) call update(natom,nrtot,nbrtot,nbrlist,rnuc,rrvec,xnbr) endif c c IF ILIST < 0 THEN WRITE XYZ FILE EVERY (-ILIST) STEPS c if ((ilist.lt.0).AND.(mod(istep,-ilist).eq.0)) then kntequil = kntequil + 1 write (suf,12000) kntequil file = 'xyz'//suf call writxyz(natom,nttot,rnuc,vnuc,nnuc,ia,zchg,time,file) endif c c c THERMALIZE SYSTEM, CALCULATE KINETIC ENERGIES c call kinet(natom,nttot,vnuc,ia,prs,amass,rk,boxk) temp = EtoK*2.0*rk/float(3*natom-NDIM) call thermal(natom,itemp,vnuc,h1,force,rk,teq,temp,gamma) rk = rk*EtoeV pressure = 0.d0 do k=1,3 pressure = pressure + abs(prs(k)) enddo pressure = pressure*EtoeV/3.d0 if(istyle.EQ.4) then call shearpull(natom,nbot,ntop,vnuc,halfvel) endif c c ... FINISH VERLET STEP c call verlet(natom,rnuc,vnuc,force,ia,amass,dt,2,rmove,QUP) c c PERFORM SUNDRY SUMMATIONS AND WHATNOT... c du = potential - oldpot c call summations(natom,rnuc,rold,force,fvold,du,dusum) energy = potential + rk c c OPTIONALLY PRINT INFORMATION EVERY TPRINT FS c ASSUME NO ONE WILL SPECIFY LESS THAN HUNDRETHS OF A FS, SO C WE CAN DO INTEGER MOD'S ON 100*TIME AND 100*TPRINT C it1 = 100*(time + 1.e-4) it2 = 100*(tprint + 1.e-4) if ((mod(it1,it2).eq.0).AND.(istep.ne.isold)) then open (unit=1,file='zoutput',status='old',access='append') open (unit=21,file='zstrain',status='unknown', & access='append') call fflush(6) write (1,10000) time,energy,potential,rk,pressure,temp, & zchg(1) write (20,10000) time,((h0(j,i),j=1,3),i=1,3) write (21,13000) time,(prs(j),j=1,6) close (1) call writcn(NDIM,nttot,atsym,rnuc,vnuc,h0,natom,ia, & zchg,time,'restart') call writxyz(natom,nttot,rnuc,vnuc,nnuc,ia,zchg,time, & 'xyz.dat') endif c c .... FINISH OF THE MD STEP c enddo 100 continue c c WRITE FINAL CONFIGURATION c open (unit=1,file='zoutput',status='old',access='append') open (unit=21,file='zstrain',status='old',access='append') if(istyle.eq.2.OR.istyle.eq.4) then open (unit=17,file='zpull',status='old',access='append') write (17,11000) dist,energy,bot,top close (17) endif write (1,10000) time,energy,potential,rk,pressure,temp,zchg(1) write (20,10000) time,((h0(j,i),j=1,3),i=1,3) write (21,13000) time,(prs(j),j=1,6) close (1) close (20) close (21) call writcn(ndim,nttot,atsym,rnuc,vnuc,h0,natom,ia,zchg,time, & 'restart') call writxyz(natom,nttot,rnuc,vnuc,nnuc,ia,zchg,time,'xyz.fin') stop 10000 format(f8.1,9(1x,1pe15.8)) 11000 format(f9.6,1x,7f12.6) 12000 format(I3.3) 13000 format(f8.1,9(1x,1pe15.8)) 14000 format(f6.2,1p6e16.8) end subroutine GetInp(NDIM,MAXATM,MAXTYP,natom,ucell,kcell,Omega,rnuc, & amass,nttot,ntype,nnuc,hard,chi,zeff,zeta,potparm) implicit double precision (a-h,o-z) c c use rationalized atomic scale units c c time: femtoseconds [10**(-15) sec] c c distance: Angstrom [10**(-10) m] c c mass: 10**(-28) kg c c energy: 10**(-18) J [62.41808 eV] c c velocity A/fs [10**5 m/sec] c c pressure: 10**(12) Pa (1 TPa, 10 Mbar) c parameter (EtoK=72435.42d0) parameter (EtoeV=6.241808d0) parameter (AMUtoM=16.6053d0) parameter (EtoH=0.229387274d0) parameter (AUtoA=0.52917d0) double precision ucell(3,3),kcell(3,3),Omega dimension rnuc(3,MAXATM) dimension amass(MAXATM) dimension ntype(MAXATM) dimension nnuc(MAXTYP) dimension hard(MAXTYP) dimension chi(MAXTYP) dimension zeff(MAXTYP) dimension zeta(MAXTYP) dimension potparm(10,3) dimension snuc(3) character*240 line parameter (HARTREE=27.21078584d0) parameter (AU=0.52917d0) c if 1st line begins with #, use RdInpt() call getline(5,nline,line) if (line(1:1).eq.'#') then call RdInp(NDIM,MAXATM,MAXTYP,natom,ucell,kcell,Omega,rnuc, & amass,nttot,ntype,nnuc,hard,chi,zeff,zeta,potparm) return endif c read in number of dimensions read (line,*) NDIM c read in periodic boundary conditions read (5,*) ucell c read in atomic positions natom = 0 nttot = 0 do i=1,MAXATM read (5,*) nn,snuc if (nn.le.0) goto 200 natom = natom + 1 do k=1,3 rnuc(k,i) = 0.d0 do j=1,3 rnuc(k,i) = rnuc(k,i) + snuc(j)*ucell(k,j) enddo enddo do n=1,nttot if (nn.eq.nnuc(n)) then ntype(i) = n goto 100 endif enddo nttot = nttot + 1 if (nttot.gt.MAXTYP) then print *,'Too many atom types used: ',MAXTYP stop endif nnuc(nttot) = nn ntype(i) = nttot 100 continue enddo 200 continue c read in values of hard(*),chi(*),zeff(*),zeta(*) do n=1,nttot nnuc(n) = -nnuc(n) enddo do n=1,nttot read (5,*,END=300,ERR=300) nn,cc,hh,xx,zz do i=1,nttot if (nn.eq.-nnuc(i)) then nnuc(i) = nn hard(i) = hh/EtoeV chi(i) = cc/EtoeV zeff(i) = zz c zeta(i) = 1.5d0/xx zeta(i) = xx goto 400 endif enddo 300 print *,'Error in reading Coulomb parameters' stop 400 continue enddo c set atomic masses do i=1,natom amass(i) = 0.d0 enddo c read in values of potential parameters c c potparm(1,*): A c potparm(2,*): r* c potparm(3,*): alpha c potparm(4,*): beta c potparm(5,*): B c potparm(6,*): eta c potparm(7,*): C c nptot = nttot*(nttot+1)/2 do n=1,nptot read (5,*) Ea,rstar,al,beta,Eb,eta,Ec potparm(1,n) = Ea/EtoeV potparm(2,n) = rstar potparm(3,n) = al potparm(4,n) = beta potparm(5,n) = Eb/EtoeV potparm(6,n) = eta potparm(7,n) = Ec/EtoeV enddo c set up reciprocal lattice vectors call RecLatt(NDIM,Omega,ucell,kcell) c get atomic masses do n=1,natom amass(n) = GetMass(nnuc(ntype(n))) enddo return end SUBROUTINE BSORT(ITEM,N,NDIM,KEY,INDEX,MAXIND) C....................................................................... C C USE A BINARY SORT TO ORDER ITEMS IN ASCENDING ORDER. THE ARGUMENTS C ARE C C ITEM(NDIM,N)....ITEMS TO BE SORTED. C N...............NUMBER OF ITEMS TO SORT. C NDIM............LOW INDEX DIMENSION OF ITEMS TO BE SORTED. C KEY.............POSITION OF ITEM TO BE SORTED C INDEX(*)........INDEX ARRAY WHICH ON RETURN WILL CONTAIN C IN INDEX(K) THE INDEX OF THE K-TH ITEM IN C THE ORDERED LIST. IT IS ASSUMED TO BE C INITIALIZED PRIOR TO ENTRY. C MAXIND..........DIMENSION OF INDEX(*). (.GE. 2*N). C C NC..............NUMBER OF COMPARISONS OF ITEMS DONE. C NF..............NUMBER OF INDEX FETCHES MADE. C NS..............NUMBER OF INDEX STORES MADE. C C AUTHOR: JOHN W. MINTMIRE, QUANTUM THEORY PROJECT, C UNIVERSITY OF FLORIDA, GAINESVILLE, FLORIDA C FLORIDA, GAINESVILLE, FLORIDA C REVISED: (12-APR-78) C REVISED: (18-FEB-92) C C....................................................................... INTEGER N,INDEX,MAXIND DOUBLE PRECISION ITEM(NDIM,N) DIMENSION INDEX(MAXIND) C C THE FOLLOWING SECTION DOES THE FIRST ITERATION OF BSORT C K1 = 1 L1 = 2 DO WHILE ((L1.LE.N).AND.(K1.LE.N)) KINDX = INDEX(K1) LINDX = INDEX(L1) IF (ITEM(KEY,KINDX).GT.ITEM(KEY,LINDX)) THEN INDEX(K1) = LINDX INDEX(L1) = KINDX ENDIF K1 = K1 + 2 L1 = L1 + 2 ENDDO C C END OF FIRST ITERATION C IF (N.LE.2) RETURN I1 = 0 I2 = N KOUNT = 0 ITER = 0 100 ITER = ITER + 1 I2T = 2**ITER JSAVE = I2 + 1 K1 = 1 DO WHILE (K1.LE.N) K2 = MIN(N,K1+I2T-1) L1 = K2 + 1 L2 = MIN(N,L1+I2T-1) C MERGE FIRST AND SECOND LIST TOGETHER DO WHILE ((L1.LE.L2).AND.(K1.LE.K2)) IF (ITEM(KEY,INDEX(I1+K1)).LE.ITEM(KEY,INDEX(I1+L1))) THEN INDEX(JSAVE) = INDEX(I1+K1) K1 = K1 + 1 ELSE INDEX(JSAVE) = INDEX(I1+L1) L1 = L1 + 1 ENDIF JSAVE = JSAVE + 1 ENDDO DO WHILE (K1.LE.K2) INDEX(JSAVE) = INDEX(I1+K1) K1 = K1 + 1 JSAVE = JSAVE + 1 ENDDO DO WHILE (L1.LE.L2) INDEX(JSAVE) = INDEX(I1+L1) L1 = L1 + 1 JSAVE = JSAVE + 1 ENDDO K1 = L2 + 1 ENDDO I1 = I2 I2 = N - I1 IF (2*I2T.LT.N) GOTO 100 IF (I1.NE.0) THEN DO I=1,N INDEX(I) = INDEX(I1+I) enddo ENDIF RETURN END subroutine MkLatt(ndim,nrtot,rvec,ucell,kcell,CutOff,uvec,indx, & MAXIND,MAXVEC) implicit double precision (a-h,o-z) parameter (twopi=6.2831853071795865d+00) parameter (TOLER=1.d-10) dimension rvec(0:3,MAXVEC),uvec(0:3,MAXVEC) dimension indx(MAXIND) double precision ucell(3,3),kcell(3,3) dimension ncell(3),r(3) ncell(1)=0 ncell(2)=0 ncell(3)=0 do n=1,NDIM sum = 0.d0 do k=1,3 sum = sum + kcell(k,n)**2 enddo ncell(n) = 1 + INT(Cutoff*SQRT(sum)/twopi) enddo TOL2 = CutOff**2 nrtot = 0 do n1=1,ncell(1) dist = 0.d0 do k=1,3 r(k) = n1*ucell(k,1) dist = dist + r(k)**2 enddo if (dist.le.TOL2) then nrtot = nrtot + 1 if (nrtot.gt.MAXVEC) then print *,'MkLatt-too many lattice vectors found' print *,MAXVEC stop endif uvec(0,nrtot) = SQRT(dist) do k=1,3 uvec(k,nrtot) = r(k) enddo endif enddo do n2=1,ncell(2) do n1=-ncell(1),ncell(1) dist = 0.d0 do k=1,3 r(k) = n2*ucell(k,2) + n1*ucell(k,1) dist = dist + r(k)**2 enddo if (dist.le.TOL2) then nrtot = nrtot + 1 if (nrtot.gt.MAXVEC) then print *,'MkLatt-too many lattice vectors found' print *,MAXVEC endif uvec(0,nrtot) = SQRT(dist) do k=1,3 uvec(k,nrtot) = r(k) enddo endif enddo enddo if (NDIM.ge.3) then do n3=1,ncell(3) do n2=-ncell(2),ncell(2) do n1=-ncell(1),ncell(1) dist = 0.d0 do k=1,3 r(k) = n1*ucell(k,1) + n2*ucell(k,2) & + n3*ucell(k,3) dist = dist + r(k)**2 enddo if (dist.le.TOL2) then nrtot = nrtot + 1 if (nrtot.gt.MAXVEC) then print *,'MkLatt-too many lattice vectors found' print *,MAXVEC stop endif uvec(0,nrtot) = SQRT(dist) do k=1,3 uvec(k,nrtot) = r(k) enddo endif enddo enddo enddo endif do i=1,nrtot indx(i) = i enddo call bsort(uvec,nrtot,4,1,indx,MAXIND) if (nrtot.gt.1) then do k=0,3 rvec(k,1) = uvec(k,indx(1)) enddo do n=2,nrtot do k=0,3 rvec(k,n) = uvec(k,indx(n)) enddo enddo endif return end subroutine RecLatt(NDIM,Omega,ucell,kcell) implicit double precision (a-h,o-z) parameter (twopi=6.2831853071795865d+00) double precision ucell(3,3),kcell(3,3),Omega if (NDIM.eq.3) then call Nabla(kcell,ucell(1,2),ucell(1,3)) call Nabla(kcell(1,2),ucell(1,3),ucell(1,1)) call Nabla(kcell(1,3),ucell(1,1),ucell(1,2)) Omega = 0.d0 do k=1,3 Omega = Omega + kcell(k,1)*ucell(k,1) enddo call aScale(9,twopi/Omega,kcell) Omega = ABS(Omega) else if (NDIM.eq.2) then call Nabla(kcell(1,3),ucell(1,1),ucell(1,2)) Omega = 0.d0 do k=1,3 Omega = Omega + kcell(k,3)**2 enddo call aScale(3,twopi/Omega,kcell(1,3)) call Nabla(kcell(1,1),ucell(1,2),kcell(1,3)) call Nabla(kcell(1,2),kcell(1,3),ucell(1,1)) Omega = SQRT(Omega) else if( NDIM.ne.0) then print *,'RecLatt-incorrect value of NDIM: ',NDIM stop endif return end subroutine aClear(n,a) integer n double precision a(n) do i=1,n a(i) = 0.d0 enddo return end subroutine iCopy(n,ia,ib) integer n integer ia(n),ib(n) do i=1,n ib(i) = ia(i) enddo return end subroutine VNorm(vec) implicit double precision (a-h,o-z) dimension vec(3) sum = 0.d0 do i=1,3 sum = sum + vec(i)**2 enddo sum = SQRT(sum) do i=1,3 vec(i) = vec(i)/sum enddo return end subroutine Nabla(a,b,c) implicit double precision (a-h,o-z) dimension a(3),b(3),c(3) c vector cross product a = b X c a(1) = b(2)*c(3) - b(3)*c(2) a(2) = b(3)*c(1) - b(1)*c(3) a(3) = b(1)*c(2) - b(2)*c(1) return end subroutine amxv(n,a,x,y) implicit double precision (a-h,o-z) dimension a(n,n) dimension x(n) dimension y(n) do i=1,n y(i) = 0.d0 enddo do i=1,n do j=1,n y(j) = y(j) + a(j,i)*x(i) enddo enddo return end subroutine ewald(NDIM,natom,nherm,rnuc,nbrtot,nbrlist,xnbr, & nrtot,rrvec,nktot,rkvec,Omega,alpha,vewald) implicit double precision (a-h,o-z) c calculate Ewald contributions to electrostatic energy c c NDIM: number of PBC's c natom: number of atoms in unit cell c nherm: packed dimension of vewald(*) c rnuc(3,natom): nuclear coordinates of atoms c nrtot: total number of real lattice points c rrvec(0:3,nrtot): coordinates of lattice sites c nktot: number of reciprocal lattice points c rkvec(0:3,nktot): coordinates of reciprocal lattice sites c Omega: area of unit cell c alpha: Gaussian exponent for Ewald procedure c vewald(nherm): packed array of quadratic coefficients of ese dimension rnuc(3,natom) dimension nbrlist(0:2,nbrtot) dimension xnbr(0:3,nbrtot) dimension rrvec(0:3,nrtot) dimension rkvec(0:3,nktot) dimension vewald(nherm) call aClear(nherm,vewald) c call timex(t1) call ewd(natom,nherm,rnuc,nrtot,rrvec,alpha,vewald) c call timex(t2) if (NDIM.eq.2) then call ew2(natom,nherm,rnuc,nktot,rkvec,Omega,alpha,vewald) else if (NDIM.eq.3) then call ew3(natom,nherm,rnuc,nktot,rkvec,Omega,alpha,vewald) else if (NDIM.ne.0) then print *,'*ewald3-wrong number of PBCs: ',NDIM stop endif c call timex(t3) c print *, 'ewd:',t2-t1,' ew2/3: ',t3-t2 return end subroutine ewd(natom,nherm,rnuc,nrtot,rrvec,alpha,vewald) implicit double precision (a-h,o-z) c calculate Ewald contributions to electrostatic energy c c natom: number of atoms in unit cell c nherm: packed dimension of vewald(*) c rnuc(3,natom): nuclear coordinates of atoms c nrtot: total number of real lattice points c Omega: area of unit cell c alpha: Gaussian exponent for Ewald procedure c vewald(nherm): packed array of quadratic coefficients of ese double precision derfc dimension rnuc(3,natom) dimension rrvec(0:3,nrtot) dimension rr(3) dimension vewald(nherm) parameter (TOL=1.d-16) TAU = -DLOG(TOL) PI = DACOS(-1.d0) TWOPI = 2*PI sqalph = SQRT(ALPHA) if(alpha.gt.0.0) alphk = 0.25d0/alpha x0 = -SQRT(alpha/PI) sum1 = 0.d0 do n=nrtot,1,-1 sum1 = sum1 + derfc(sqalph*rrvec(0,n))/rrvec(0,n) enddo zsum = 2*(x0+sum1) ij = 0 do i=1,natom do j=1,i-1 dist = 0.d0 do k=1,3 rr(k) = rnuc(k,i) - rnuc(k,j) dist = dist + rr(k)**2 enddo dist = SQRT(dist) sum1 = 0.d0 do n=nrtot,1,-1 dist1 = 0.d0 dist2 = 0.d0 do k=1,3 dist1 = dist1 + (rr(k)+rrvec(k,n))**2 dist2 = dist2 + (rr(k)-rrvec(k,n))**2 enddo if (alpha*dist1.lt.TAU) then dist1 = SQRT(dist1) sum1 = sum1 + derfc(sqalph*dist1)/dist1 endif if (alpha*dist2.lt.TAU) then dist2 = SQRT(dist2) sum1 = sum1 + derfc(sqalph*dist2)/dist2 endif enddo sum1 = sum1 + derfc(sqalph*dist)/dist ij = ij + 1 vewald(ij) = vewald(ij) + sum1 enddo ij = ij + 1 vewald(ij) = vewald(ij) + zsum enddo return end subroutine ew3(natom,nherm,rnuc,nktot,rkvec,Omega,alpha,vewald) implicit double precision (a-h,o-z) c calculate Ewald contributions to electrostatic energy c c natom: number of atoms in unit cell c nherm: packed dimension of vewald(*) c rnuc(3,natom): nuclear coordinates of atoms c nktot: number of reciprocal lattice points c rkvec(0:3,nktot): coordinates of reciprocal lattice sites c Omega: area of unit cell c alpha: Gaussian exponent for Ewald procedure c vewald(nherm): packed array of quadratic coefficients of ese dimension rnuc(3,natom) dimension rkvec(0:3,nktot) dimension vewald(nherm) PI = DACOS(-1.d0) TWOPI = 2*PI sqalph = SQRT(ALPHA) alphk = 0.25d0/alpha piom = PI/Omega zsum = 0.d0 do n=nktot,1,-1 dd = rkvec(0,n)*rkvec(0,n) factor = (8*piom)*EXP(-alphk*dd)/dd zsum = zsum + factor ij = 0 do i=1,natom rphase = 0.d0 do k=1,3 rphase = rphase+rkvec(k,n)*rnuc(k,i) enddo do j=1,i-1 phase = rphase do k=1,3 phase = phase-rkvec(k,n)*rnuc(k,j) enddo vewald(ij+j) = vewald(ij+j) + factor*COS(phase) enddo ij = ij+i enddo enddo ij = 0 do i=1,natom ij = ij + i vewald(ij) = vewald(ij) + zsum enddo return end subroutine ew2(natom,nherm,rnuc,nktot,rkvec,Omega,alpha,vewald) implicit double precision (a-h,o-z) c calculate Ewald contributions to electrostatic energy c c natom: number of atoms in unit cell c nherm: packed dimension of vewald(*) c rnuc(3,natom): nuclear coordinates of atoms c nktot: number of reciprocal lattice points c rkvec(0:3,nktot): coordinates of reciprocal lattice sites c Omega: area of unit cell c alpha: Gaussian exponent for Ewald procedure c vewald(nherm): packed array of quadratic coefficients of ese dimension rnuc(3,natom) dimension rkvec(0:3,nktot) dimension rr(3) dimension vewald(nherm) parameter (TOL = 1.d-12) PI = DACOS(-1.d0) TWOPI = 2*PI sqalph = SQRT(ALPHA) TAU = -DLOG(TOL) y0 = -SQRT(PI/alpha)/Omega piom = pi/Omega alphK = 0.5d0/sqalph sum2 = 0.d0 do n=nktot,1,-1 aa = alphK*rkvec(0,n) term = derfc(aa)/rkvec(0,n) sum2 = sum2 + term enddo zsum = 2*(y0+2*piom*sum2) ij = 0 c call timex(t0) do i=1,natom c call timex(t1) do j=1,i-1 do k=1,3 rr(k) = rnuc(k,i) - rnuc(k,j) enddo zmn = rr(3) azmn = sqalph*zmn azmn2 = azmn*azmn xp = EXP(-azmn2) term1 = 2*piom*zmn*(1.d0-derfc(azmn)) term2 = -2*y0*xp sum0 = term1 + term2 sum2 = 0.d0 if (xp.gt.TOL) then do n=1,nktot aa = alphK*rkvec(0,n) zk = zmn*rkvec(0,n) if (MAX(aa,ABS(zk)).lt.TAU) then pre = EXP(zk) term = derfc(aa-azmn)/pre term = term + derfc(aa+azmn)*pre phase = rkvec(1,n)*rr(1) + rkvec(2,n)*rr(2) xt = DCOS(phase)*term/rkvec(0,n) sum2 = sum2 + xt else goto 100 endif enddo 100 sum2 = 2*piom*sum2 endif vewald(ij+j) = vewald(ij+j) + sum2 - sum0 enddo ij = ij + i c call timex(t2) c print *, ij, t2-t1,t2-t0 vewald(ij) = vewald(ij) + zsum enddo return end function ewsum(natom,znuc,vewald) implicit double precision (a-h,o-z) dimension znuc(natom),vewald(*) v2 = 0.d0 ij = 0 do i=1,natom do j=1,i-1 ij = ij + 1 v2 = v2 + 2*znuc(i)*znuc(j)*vewald(ij) enddo ij = ij + 1 v2 = v2 + vewald(ij)*znuc(i)*znuc(i) enddo ewsum = v2/2.d0 return end subroutine GetMem(pntr,mem) pointer (pntr,v) integer mem pntr = malloc(mem) if (pntr.le.0) then print *,'GetMem is unable to obtain memory ',mem stop else print *,'Getting memory: ',mem endif return end subroutine FreMem(pntr,mem) pointer (pntr,v) call free(pntr) return end subroutine abort() print *,'aborting . . .' stop end subroutine timex(T) implicit double precision (a-h,o-z) real tarray(2) c UNIX OS version call etime(tarray) T = tarray(1) + tarray(2) return end subroutine GetLib(dirsymb,filename,libname) implicit integer(a-z) c c GetLib takes an input character string "dirsymb" which is c the symbolic reference to where a library directory is c located, and returns in "libname" the ASCII string for c the directory c character*(*) libname character*(*) dirsymb character*(*) filename call getenv(dirsymb,libname) ntmp = len(libname) do while ((ntmp.gt.0).and.(libname(ntmp:ntmp).eq.' ')) ntmp = ntmp - 1 enddo if (ntmp.gt.0) then ntmp = ntmp + 1 libname(ntmp:ntmp) = '/' endif libname = libname(1:ntmp)//filename return end subroutine esnewton(NDIM,natom,nherm,strain,nbrtot,nbrlist,xnbr, & nrtot,rrvec,nktot,rkvec,ucell,kcell,Omega,rnuc,zchg,nttot, & nptot,potparm,ntype,zeff,zeta,alphew,pvec,xvec,dqdeps,vewald, & catom,pden,work) implicit double precision (a-h,o-z) dimension strain(6) dimension nbrlist(0:2,nbrtot) dimension xnbr(0:3,nbrtot) dimension rrvec(0:3,nrtot) dimension rkvec(0:3,nktot) double precision ucell(3,3),kcell(3,3) dimension rnuc(3,natom) dimension zchg(natom) dimension ntype(natom) dimension potparm(10,nptot) dimension zeff(nttot) dimension zeta(nttot) dimension pvec(natom,6) dimension xvec(natom,6) dimension dqdeps(natom,6) dimension vewald(nherm) dimension catom(natom) dimension pden(natom) dimension work(natom,7) dimension prs(6) dimension esprs(6) dimension esijkl(21) dimension vijkl(21) dimension cijkl(21) dimension cpair(21),ppair(6) dimension istrain(6) call timex(t1) call esel(NDIM,natom,nbrtot,nbrlist,xnbr,nrtot,rrvec,nktot,rkvec, & Omega,rnuc,zchg,ntype,zeff,zeta,alphew,pvec,xvec,esijkl,vijkl) call timex(t2) call elpair(natom,nptot,ntype,nbrtot,nbrlist,xnbr,potparm,ppair, & cpair,pden,work) call timex(t3) call SvStrain(natom,vewald,catom,zchg,dqdeps,pvec,xvec,esprs, & esijkl,work) call aCopy(6,esprs,1,prs) call aCopy(21,esijkl,1,cijkl) call timex(t4) do m=1,6 prs(m) = prs(m) + ppair(m) strain(m) = -prs(m) enddo do m=1,21 cijkl(m) = cijkl(m) + cpair(m) vijkl(m) = cijkl(m) enddo c calculate strains call DSPSV('U',6,1,vijkl,istrain,strain,6,info) call MkStrain(NDIM,natom,strain,ucell,kcell,Omega,rnuc,nrtot, & rrvec,nktot,rkvec) call timex(t5) print 10000,strain print * print 11000,t5-t1,t2-t1,t3-t2,t4-t3,t5-t4 call flush(6) return 10000 format(6f12.6) 11000 format('esnewton timings: ',5f10.2) end subroutine SvStrain(natom,vewald,catom,zchg,dqdeps,pvec,xvec,prs, & cijkl,work) implicit double precision (a-h,o-z) dimension vewald(*) dimension catom(natom) dimension zchg(natom) dimension dqdeps(6,natom) dimension pvec(6,natom) dimension xvec(6,natom) dimension prs(6) dimension cijkl(21) dimension work(natom,0:6) dimension dlambda(6) call aclear(6,prs) do i=1,natom work(i,0) = 1.d0 do m=1,6 work(i,m) = xvec(m,i) + pvec(m,i) enddo enddo call DSPSV('U',natom,7,vewald,dqdeps,work,natom,info) sum = 0.d0 do i=1,natom sum = sum + work(i,0) enddo do m=1,6 ssum = 0.d0 do i=1,natom ssum = ssum + work(i,m) enddo dlambda(m) = ssum/sum do i=1,natom dqdeps(m,i) = -work(i,m) + dlambda(m)*work(i,0) enddo enddo ij = 0 do i=1,6 do m=1,natom prs(i) = prs(i) + zchg(m)*(pvec(i,m)/2.d0+xvec(i,m)) enddo do j=1,i ij = ij + 1 do m=1,natom cijkl(ij) = cijkl(ij) & + 0.5d0*(dqdeps(i,m)*(pvec(j,m)+xvec(j,m)) & +dqdeps(j,m)*(pvec(i,m)+xvec(i,m))) enddo enddo enddo return end subroutine makenbr(MAXNBR,NDIM,natom,ucell,rnuc,nrtot,rrvec,rcut, & nbrtot,nbrlist,work) implicit double precision (a-h,o-z) c makes neighbor list using cutoff value c c (input parameters) c MAXNBR: maximum number of neighbors dimensioned c NDIM: number of Cartesion directions for periodicity c (can be either 2 or 3) c natom: number of dimer pairs c ucell: periodic repeat distance array c rnuc: nuclear coordinates c nrtot: total number of real periodic lattice vectors c rrvec: array of lattice vectors c rcut: cutoff distance for interdimer potential c c (output parameters) c nbrtot: number of neighbor list elements found c nbrlist: neighbor list stored as indices of neighbor pairs c and cell of second neighbor c c (work arrays) c work: stores intermediate separation vectors c c this current version DOES NOT use the minimum image convention dimension rnuc(3,natom) dimension ucell(3,NDIM) dimension rrvec(0:3,nrtot) dimension nbrlist(0:2,MAXNBR) dimension work(3,natom) nbrtot = 0 rcut2 = rcut*rcut call timex(tstart) do i=2,natom do j=1,i-1 work(1,j) = rnuc(1,j) - rnuc(1,i) work(2,j) = rnuc(2,j) - rnuc(2,i) work(3,j) = rnuc(3,j) - rnuc(3,i) enddo do j=1,i-1 dist = 0.d0 do k=1,3 dist = dist + work(k,j)**2 enddo if (dist.lt.rcut2.or.NDIM.eq.0) then nbrtot = nbrtot + 1 if (nbrtot.le.MAXNBR) then nbrlist(0,nbrtot) = 0 nbrlist(1,nbrtot) = i nbrlist(2,nbrtot) = j endif endif enddo enddo do n=1,nrtot do i=1,natom work(1,i) = rnuc(1,i) - rrvec(1,n) work(2,i) = rnuc(2,i) - rrvec(2,n) work(3,i) = rnuc(3,i) - rrvec(3,n) enddo do i=1,natom do j=1,natom dist = 0.d0 do k=1,3 dist = dist + (work(k,i)-rnuc(k,j))**2 enddo if (dist.lt.rcut2) then nbrtot = nbrtot + 1 if (nbrtot.le.MAXNBR) then nbrlist(0,nbrtot) = n nbrlist(1,nbrtot) = i nbrlist(2,nbrtot) = j endif endif enddo enddo enddo if (nbrtot.gt.MAXNBR) then print *,'MakeNbr-MAXNBR exceeded, need ',nbrtot,' neighbors' stop endif call timex(tend) tend = tend - tstart print 10000,nbrtot,tend return 10000 format('MakeNbr found ',i8,' neighbors (',f8.1,' secs)') end subroutine RdInp(NDIM,MAXATM,MAXTYP,natom,ucell,kcell,Omega,rnuc, & amass,nttot,ntype,nnuc,hard,chi,zeff,zeta,potparm) implicit double precision (a-h,o-z) c c use rationalized atomic scale units c c time: femtoseconds [10**(-15) sec] c c distance: Angstrom [10**(-10) m] c c mass: 10**(-28) kg c c energy: 10**(-18) J [62.41808 eV] c c velocity A/fs [10**5 m/sec] c c pressure: 10**(12) Pa (1 TPa, 10 Mbar) c parameter (EtoK=72435.42d0) parameter (EtoeV=6.241808d0) parameter (AMUtoM=16.6053d0) parameter (EtoH=0.229387274d0) parameter (AUtoA=0.52917d0) double precision ucell(3,3),kcell(3,3),Omega dimension rnuc(3,MAXATM) dimension amass(MAXATM) dimension ntype(MAXATM) dimension nnuc(MAXTYP) dimension hard(MAXTYP) dimension chi(MAXTYP) dimension zeff(MAXTYP) dimension zeta(MAXTYP) dimension potparm(10,3) integer GtAtNo character*2 symbol character*240 line parameter (HARTREE=27.21078584d0) parameter (AU=0.52917d0) c read in number of dimensions read (5,*) NDIM c read in periodic boundary conditions read (5,*) ucell c read in atomic positions natom = 0 nttot = 0 call getline(5,nline,line) nn = 0 do while ((nline.gt.0).and.(natom.lt.MAXATM)) natom = natom + 1 call GtLbl(nn,nline,line,symbol) do k=1,3 call skipwt(nn,nline,line) call GtReal(nn,nline,line,rnuc(k,natom)) enddo nz = GtAtNo(symbol) do n=1,nttot if (nz.eq.nnuc(n)) then ntype(natom) = n goto 100 endif enddo nttot = nttot + 1 if (nttot.gt.MAXTYP) then print *,'Too many atom types used: ',MAXTYP stop endif nnuc(nttot) = nz ntype(natom) = nttot 100 continue call getline(5,nline,line) nn = 0 if ((natom.ge.MAXATM).and.(nline.gt.0)) then print *,'Number of atoms limit exceeded' stop endif enddo c read in potential parameters call RdPot(nttot,nnuc,chi,hard,zeta,zeff,potparm) c set up reciprocal lattice vectors call RecLatt(NDIM,Omega,ucell,kcell) c get atomic masses do n=1,natom amass(n) = GetMass(nnuc(ntype(n))) enddo return end subroutine getline(nunit,n,line) character*(*) line integer n logical isspace n = 0 read (nunit,10000,end=100) line n = len(line) do while ((n.gt.0).and.isspace(line(n:n))) n = n - 1 enddo return 100 n = -1 return 10000 format(a) end subroutine skipnwt(n,nline,line) integer n,nline character*240 line logical isspace do while ((n.lt.nline).and..not.isspace(line(n+1:n+1))) n = n + 1 enddo return end subroutine skipwt(n,nline,line) integer n,nline character*240 line logical isspace do while ((n.lt.nline).and.isspace(line(n+1:n+1))) n = n + 1 enddo return end logical function strcmp(n,a,b) character*(*) a,b logical islower character*1 xy,x,y,UPPER UPPER(xy) = CHAR(ICHAR(xy)+ICHAR('A')-ICHAR('a')) strcmp = .FALSE. nlimit = MIN(n,LEN(a),LEN(b)) do i=1,n if (islower(a(i:i))) then x = UPPER(a(i:i)) else x = a(i:i) endif if (islower(b(i:i))) then y = UPPER(b(i:i)) else y = b(i:i) endif if (x.ne.y) return enddo strcmp = .TRUE. return end logical function isalnum(c) character*1 c logical isalpha,isdigit isalnum = isalpha(c).or.isdigit(c) return end logical function isalpha(c) character*1 c logical isupper,islower isalpha = isupper(c).or.islower(c) return end logical function isdigit(c) character*1 c isdigit = (index('1234567890',c).gt.0) return end logical function islower(c) character*1 c integer ic ic = ICHAR(c) islower = ((ic.ge.ICHAR('a')).and.(ic.le.ICHAR('z'))) return end logical function isspace(c) character*1 c character*1 BLANK parameter (BLANK=' ') if ((c.eq.BLANK).or.(ICHAR(c).lt.32)) then isspace = .TRUE. else isspace = .FALSE. endif return end logical function isupper(c) character*1 c integer ic ic = ICHAR(c) isupper = ((ic.ge.ICHAR('A')).and.(ic.le.ICHAR('Z'))) return end integer function GtAtNo(string) implicit integer(a-z) c*************************************************************************** c * c GTATNO.FOR - returns atomic number if string begins with valid * c atomic symbol, 0 if string begins with X (dummy label)* c or * (supporting obsolete usage) , and -1 if * c it can't figure out what the string means * c * c Author: J. W. Mintmire * c Naval Research Laboratory * c Washington, DC 20375 * c * c Date: January 24, 1986 * c Rvsd: April 19, 1990 * c * c*************************************************************************** character*(*) string character*1 UPPER,x character*2 Nstrng,Name logical isalpha,strcmp parameter (MXNAME=101) character*2 RefNam(0:MXNAME) data RefNam/'X ','H ','He','Li','Be','B ','C ','N','O ','F ','Ne', & 'Na','Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti', & 'V ','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se', & 'Br','Kr','Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru','Rh','Pd', & 'Ag','Cd','In','Sn','Sb','Te','I ','Xe','Cs','Ba','La','Ce', & 'Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb', & 'Lu','Hf','Ta','W ','Re','Os','Ir','Pt','Au','Hg','Tl','Pb', & 'Bi','Po','At','Fr','Rn','Ra','Ac','Th','Pa','U ','Np','Pu', & 'Am','Cm','Bk','Cf','Es','Fm','Md'/ UPPER(x) = CHAR(ICHAR(x)+ICHAR('A')-ICHAR('a')) if (string(1:1).eq.'*') then GtAtNo = 0 return else if (isalpha(string(1:1))) then Name = string(1:1) if (isalpha(string(2:2))) then Name(2:2) = string(2:2) endif else print *,'GtAtNo-invalid string encountered' stop endif do i=0,MXNAME if (strcmp(2,Name,RefNam(i))) then GtAtNo = i return endif enddo GtAtNo = -1 return Entry GtAtSy(Nstrng,IatNo) Nstrng = RefNam(IatNo) return end subroutine GtIntg(n,nline,line,intg) implicit double precision (a-h,o-z) character*240 line character*20 string logical isspace i = 0 do while ((n+i.lt.nline).and..not.isspace(line(n+i+1:n+i+1))) i = i + 1 enddo if (i.le.0) then print *,'GtIntg-no integer found' print *,line(1:nline) stop else if (i.gt.20) then print *,'GtIntg-invalid integer string found' print *,line(1:nline) stop endif string = ' ' string(21-i:20) = line(n+1:n+i) n = n + i call skipwt(n,nline,line) read (string,*,err=100) intg return 100 print *,'GtIntg-error reading integer string, input line follows' print *,line(1:nline) call ABORT 10000 format(i20) end subroutine GtLbl(n,nline,line,zsymb) implicit double precision (a-h,o-z) character*240 line character*(*) zsymb logical isalpha,isalnum,isspace i = 0 do while ((i.lt.MIN(len(zsymb),nline-n)).and. & isalnum(line(n+i+1:n+i+1))) c check first character to see if alphabetic i = i + 1 if (i.eq.1) then if (.NOT.isalpha(line(n+i:n+i))) then print *,'GtLbl-invalid first character of label' print *,line(1:nline) stop endif endif enddo zsymb = ' ' zsymb = line(n+1:n+i) n = n + i if (.not.isspace(line(n+1:n+1))) then print *,'GtLbl-invalid character at end of nuclear label' print *,line(1:nline) stop endif call skipwt(n,nline,line) return end subroutine GtReal(n,nline,line,value) implicit double precision (a-h,o-z) character*240 line character*30 string logical isspace i = 0 do while ((n+i.lt.nline).and..not.isspace(line(n+i+1:n+i+1))) i = i + 1 enddo if (i.le.0) then print *,'GtReal-no real value found' print *,line(1:nline) stop else if (i.gt.30) then print *,'GtReal-invalid real string found' print *,line(1:nline) stop endif string = ' ' string(31-i:30) = line(n+1:n+i) n = n + i call skipwt(n,nline,line) read (string,*,err=100) value return 100 print *,'GtReal-error reading float, input line follows' print *,line(1:nline) call ABORT 10000 format(f30.4) end function GetMass(n) implicit double precision (a-h,o-z) parameter (AMUtoM=16.6053d0) if (n.eq.1) then GetMass = 1.0079d0 else if (n.eq.8) then GetMass = 15.9994d0 else if (n.eq.13) then GetMass = 26.98154d0 else if (n.eq.22) then GetMass = 47.90 else if (n.eq.29) then GetMass = 63.546d0 else print *,'GetMass: mass not entered for atomic number ',n stop endif GetMass = AMUtoM*GetMass return end subroutine ewf(NDIM,natom,rnuc,znuc,nbrtot,nbrlist,xnbr,nrtot, & rrvec,nktot,rkvec,Omega,alphew,force) implicit double precision (a-h,o-z) c calculate Ewald contributions to electrostatic energy c c natom: number of atoms in unit cell c c rcut: Ewald cutoff distance c c rnuc(3,natom): nuclear coordinates of atoms c c znuc(natom): nuclear charges c c nrtot: total number of real lattice points c c rrvec(0:3,nrtot): coordinates of lattice sites c c nktot: number of reciprocal lattice points c c rkvec(0:3,nktot): coordinates of reciprocal lattice sites c c Omega: area of unit cell c c alphew: Gaussian exponent for Ewald procedure c c force(3,natom) calculated forces from Ewald summation c dimension rnuc(3,natom) dimension znuc(natom) dimension nbrlist(0:2,nbrtot) dimension xnbr(0:3,nbrtot) dimension rrvec(0:3,nrtot) dimension rkvec(0:3,nktot) dimension force(3,natom) do i=1,natom force(1,i) = 0.d0 force(2,i) = 0.d0 force(3,i) = 0.d0 enddo call ewfd(natom,rnuc,znuc,nrtot,rrvec,alphew,force) if (NDIM.eq.2) then call ewf2(natom,rnuc,znuc,nktot,rkvec,Omega,alphew,force) else if (NDIM.eq.3) then call ewf3(natom,rnuc,znuc,nktot,rkvec,Omega,alphew,force) else if (NDIM.ne.0) then print *,'ewf: invalid number of dimensions: ',NDIM stop endif do i=1,natom do k=1,3 force(k,i) = znuc(i)*force(k,i) enddo enddo return end subroutine ewf2(natom,rnuc,znuc,nktot,rkvec,Omega,alpha,force) implicit double precision (a-h,o-z) c calculate Ewald contributions to electrostatic energy c c natom: number of atoms in unit cell c rnuc(3,natom): nuclear coordinates of atoms c znuc(natom): nuclear charges c nktot: number of reciprocal lattice points c rkvec(0:3,nktot): coordinates of reciprocal lattice sites c Omega: area of unit cell c alpha: Gaussian exponent for Ewald procedure c force(3,natom): forces from direct summation dimension rnuc(3,natom),znuc(natom) dimension rkvec(0:3,nktot) dimension rr(3),vv(3) dimension force(3,natom) parameter (TOL=1.d-12) PI = 3.14159265358979323846 HPI = PI/2.0 sqalph = SQRT(alpha) alphk = 0.5d0/sqalph zfact = sqalph/SQRT(PI) TAU = -DLOG(TOL) scale = 2*PI/Omega c call timex(t1) do i=1,natom do j=1,i-1 do k=1,3 rr(k) = rnuc(k,i) - rnuc(k,j) enddo zmn = rr(3) azmn = sqalph*zmn c now add K=0 components vv(1) = 0.d0 vv(2) = 0.d0 vv(3) = -scale*(1.d0-derfc(azmn)) c add K .ne. 0 components do n=1,nktot zk = zmn*rkvec(0,n) if (ABS(zk).lt.TAU) then pre = EXP(zk) aa = alphK*rkvec(0,n) aa1 = aa + azmn aa2 = aa - azmn term1 = pre*derfc(aa1) term2 = derfc(aa2)/pre dt1 = pre*EXP(-aa1*aa1) dt2 = EXP(-aa2*aa2)/pre phase = rkvec(1,n)*rr(1) + rkvec(2,n)*rr(2) xyterm = -(term1+term2)*COS(phase-HPI)/rkvec(0,n) zterm1 = term1 - term2 zterm2 = zfact*(dt2-dt1)/rkvec(0,n) zterm = (zterm1+zterm2)*COS(phase) xyterm = xyterm*scale zterm = zterm*scale vv(1) = vv(1) + xyterm*rkvec(1,n) vv(2) = vv(2) + xyterm*rkvec(2,n) vv(3) = vv(3) + zterm else goto 100 endif enddo 100 continue do k=1,3 force(k,i) = force(k,i) + vv(k)*znuc(j) force(k,j) = force(k,j) - vv(k)*znuc(i) enddo enddo enddo c call timex(t2) c print *, 'ewf2 total: ',t2-t1 return end subroutine ewf3(natom,rnuc,znuc,nktot,rkvec,Omega,alpha,force) implicit double precision (a-h,o-z) c calculate Ewald contributions to electrostatic energy c c natom: number of atoms in unit cell c rnuc(3,natom): nuclear coordinates of atoms c znuc(natom): nuclear charges c nktot: number of reciprocal lattice points c rkvec(0:3,nktot): coordinates of reciprocal lattice sites c Omega: area of unit cell c alpha: Gaussian exponent for Ewald procedure c force(3,natom): forces from direct summation dimension rnuc(3,natom),znuc(natom) dimension rkvec(0:3,nktot) dimension force(3,natom) PI = 3.14159265358979323846 HPI=PI/2.0 alphk = 0.25d0/alpha pre = 8*PI/Omega do n=nktot,1,-1 d2 = rkvec(0,n)**2 term = pre*EXP(-alphk*d2)/d2 do i=1,natom do j=1,i-1 phase = 0.d0 do k=1,3 phase = phase + rkvec(k,n)*(rnuc(k,i)-rnuc(k,j)) enddo x = term*COS(phase-HPI) do k=1,3 force(k,i) = force(k,i) - x*znuc(j)*rkvec(k,n) force(k,j) = force(k,j) + x*znuc(i)*rkvec(k,n) enddo enddo enddo enddo return end subroutine ewfd(natom,rnuc,znuc,nrtot,rrvec,alpha,force) implicit double precision (a-h,o-z) c calculate Ewald contributions to electrostatic energy c c natom: number of atoms in unit cell c rnuc(3,natom): nuclear coordinates of atoms c znuc(natom): nuclear charges c nrtot: total number of real lattice points c rrvec(0:3,nrtot): coordinates of lattice sites c alpha: Gaussian exponent for Ewald procedure c force(3,natom): forces from direct summation dimension rnuc(3,natom),znuc(natom) dimension rrvec(0:3,nrtot) dimension rr(3),rr1(3),rr2(3),v(3) dimension force(3,natom) parameter (TOL=1.d-12) PI = 3.14159265358979323846 sqalph = SQRT(alpha) sqa2pi = SQRT(alpha/PI) TAU = -DLOG(TOL) do i=1,natom do j=1,i-1 dd = 0.d0 do k=1,3 rr(k) = rnuc(k,i) - rnuc(k,j) dd = dd + rr(k)**2 v(k) = 0.d0 enddo do n=nrtot,1,-1 dd1 = 0.d0 dd2 = 0.d0 do k=1,3 rr1(k) = rr(k) - rrvec(k,n) rr2(k) = rr(k) + rrvec(k,n) dd1 = dd1 + rr1(k)**2 dd2 = dd2 + rr2(k)**2 enddo add1 = alpha*dd1 add2 = alpha*dd2 if (add1.lt.TAU) then dist1 = SQRT(dd1) t1 = (2*sqa2pi*dist1*EXP(-alpha*dd1) & +derfc(sqalph*dist1))/dd1 t = t1/dist1 v(1) = v(1) + t*rr1(1) v(2) = v(2) + t*rr1(2) v(3) = v(3) + t*rr1(3) endif if (add2.lt.TAU) then dist2 = SQRT(dd2) t2 = (2*sqa2pi*dist2*EXP(-alpha*dd2) & +derfc(sqalph*dist2))/dd2 t = t2/dist2 v(1) = v(1) + t*rr2(1) v(2) = v(2) + t*rr2(2) v(3) = v(3) + t*rr2(3) endif enddo dist = SQRT(dd) t = (2*sqa2pi*dist*EXP(-alpha*dd)+derfc(sqalph*dist))/dd do k=1,3 v(k) = v(k) + t*rr(k)/dist force(k,i) = force(k,i) - v(k)*znuc(j) force(k,j) = force(k,j) + v(k)*znuc(i) enddo enddo enddo return end subroutine GetPar(alphew,rcut,rmax,ecut,NDIM,Omega,ucell,kcell) implicit double precision (a-h,o-z) double precision ucell(3,3),kcell(3,3) c RATIO is the ratio of reciprocal lattice points c relative to the direct lattice points used in c the Ewald summation parameter (TOL=1.d-25) parameter (CUTMAX=30.d0) parameter (CUTMIN = 20.d0) PI = 3.14159265358979323846 if (NDIM.eq.3) then RATIO = 3.d0 else RATIO = 1.d0 endif tau = -DLOG(TOL) rmax = 2.0 rcut = CUTMIN ecut = 2.0 alphew = 0.0 if (NDIM.eq.2) then rmax = SQRT(Omega) alphew = SQRT(RATIO) else if (NDIM.eq.3) then rmax = Omega**(1.d0/3.d0) alphew = RATIO**(1.d0/3.d0) else if (NDIM.eq.0) then return else print *,'GetPar: invalid number of dimensions ',NDIM stop endif alphew = PI*alphew/rmax**2 c estimate rcut based on current optimized alphew rcut = SQRT(tau/alphew) c make sure rcut is at least greater than all primitive c lattice vectors and at least equal to CUTMAX cmax = 0.d0 do m=1,NDIM sum = 0.d0 do k=1,3 sum = sum + ucell(k,m)**2 enddo sum = SQRT(sum) cmax = MAX(sum,cmax) enddo rcut = MAX(rcut,CUTMIN) c set alphew and ecut based on current values of rcut alphew = tau/rcut**2 ecut = 2*SQRT(tau*alphew) rcut = rcut + cmax rmax = CUTMIN print *,'Ewald alpha value: ',alphew print *,'Ewald real lattice cutoff: ',rcut print *,'Ewald reciprocal lattice cutoff: ',ecut print *,'Real space cutoff: ',rmax return end subroutine esforce(NDIM,natom,nbrtot,nbrlist,xnbr,nrtot,rrvec, & nktot,rkvec,Omega,rnuc,zchg,ntype,zeff,zeta,alphew,force) implicit double precision (a-h,o-z) c c use rationalized atomic scale units c c time: femtoseconds [10**(-15) sec] c c distance: Angstrom [10**(-10) m] c c mass: 10**(-28) kg c c energy: 10**(-18) J [6.241808 eV] c c velocity A/fs [10**5 m/sec] c c pressure: 10**(13) Pa (10 TPa, 100 Mbar) c parameter (EtoK=72435.42d0) parameter (EtoeV=6.241808d0) parameter (AMUtoM=16.6053d0) parameter (EtoH=0.229387274d0) parameter (AUtoA=0.52917d0) dimension nbrlist(0:2,nbrtot) dimension xnbr(0:3,nbrtot) dimension rrvec(0:3,nrtot) dimension rkvec(0:3,nktot) dimension zchg(natom) dimension rnuc(3,natom),ntype(natom) dimension zeff(*),zeta(*) dimension force(3,natom) dimension cpu(0:2) call timex(t1) call ewf(NDIM,natom,rnuc,zchg,nbrtot,nbrlist,xnbr,nrtot,rrvec, & nktot,rkvec,Omega,alphew,force) call timex(t2) do m=1,nbrtot i = nbrlist(1,m) j = nbrlist(2,m) ni = ntype(i) nj = ntype(j) dist = xnbr(0,m) call CalVij(dist,zeta(ni),zeta(nj),vij,dvij,d2vij) call CalVn(dist,zeta(ni),vi,dvi,d2vi) call CalVn(dist,zeta(nj),vj,dvj,d2vj) factor = zchg(i)*zchg(j)*dvij + zeff(ni)*zchg(j)*(dvj-dvij) & + zeff(nj)*zchg(i)*(dvi-dvij) factor = factor/dist do k=1,3 force(k,i) = force(k,i) + factor*xnbr(k,m) force(k,j) = force(k,j) - factor*xnbr(k,m) enddo enddo call timex(t3) call aScale(3*natom,AUtoA/EtoH,force) cpu(0) = t3-t1 cpu(1) = t2-t1 cpu(2) = t3-t2 c print 111, cpu 111 format('total:',f8.3,' ewf:',f8.3,' nbrtot:',f8.3) return end subroutine CalVij(R,a,b,vab,vabd,vabd2) implicit double precision (a-h,o-z) double precision kappa parameter (TOL=1.d-8) parameter (RMAX=18.d0) tau = (a-b) if (abs(tau).lt.TOL) then rho = a*R xprho = EXP(-2*rho) vab = -xprho*(1+rho*(11.d0/8.d0+rho*(3.d0/4.d0+rho/6.d0)))/R vabd = -vab/R - 2*a*vab - & a*xprho*(11.d0/8.d0+rho*1.5d0+rho*rho/2.d0)/R vabd2 = -(2*vabd+2*a*vab)/R - 2*a*vabd & + a*a*xprho*(11.d0/4.d0+3*rho+rho*rho)/R & - a*a*xprho*(1.5d0+rho)/R else tau = tau/(a+b) rhoa = a*R rhob = b*R xpa = EXP(-2*rhoa) xpb = EXP(-2*rhob) kappa = (tau+1/tau)/2 vab = -(xpa*(2+kappa+rhoa)*(1-kappa)**2+xpb*(2-kappa+rhob)* & (1+kappa)**2)/R/4.d0 vabd = (-vab+0.5d0* & (a*xpa*(1.5d0+kappa+rhoa)*(1-kappa)**2+b*xpb* & (1.5d0-kappa+rhob)*(1+kappa)**2))/R vabd2 = -(2*vabd+a*xpa*(a*(1+kappa+rhoa)*(1-kappa)**2)+b*xpb* & (b*(1-kappa+rhob)*(1+kappa)**2))/R endif return end subroutine CalVn(R,a,v,vd,vd2) implicit double precision (a-h,o-z) parameter (RMAX=18.d0) c this routine needs testing rho = a*R exprho = 0.d0 exprho = EXP(-2*rho) v = -(1+rho)*exprho/R vd = -v/R + a*(1+2*rho)*exprho/R vd2 = -2*vd/R - 4*a*a*rho*exprho/R return end subroutine getvij(R,a,b,va,vb,vab) implicit double precision (a-h,o-z) double precision kappa parameter (TOL=1.d-8) parameter (RMAX=18.d0) tau = (a-b) if (abs(tau).lt.TOL) then rho = a*R xprho = 0.d0 xprho = EXP(-2*rho) va = -(1+rho)*xprho/R vb = va vab = -xprho*(1+rho*(11.d0/8.d0+rho*(3.d0/4.d0+rho/6.d0)))/R else tau = tau/(a+b) rhoa = a*R rhob = b*R xpa = EXP(-2*rhoa) xpb = EXP(-2*rhob) va = -(1+rhoa)*xpa/R vb = -(1+rhob)*xpb/R kappa = (tau+1/tau)/2 vab = -(xpa*(2+kappa+rhoa)*(1-kappa)**2/4.d0+xpb* & (2-kappa+rhob)*(1+kappa)**2/4.d0)/R endif return end subroutine epforce(natom,nbrtot,ntype,nbrlist,xnbr,potparm,epair, & embed,force,pden) implicit double precision (a-h,o-z) dimension ntype(natom) dimension nbrlist(0:2,nbrtot) dimension xnbr(0:3,nbrtot) dimension potparm(10,3) dimension force(3,natom) dimension pden(0:3,natom) dimension npot(2,2) data npot/1,3,3,2/ epair = 0.d0 call GetDen(natom,nbrtot,nbrlist,xnbr,ntype,potparm,pden) do m=1,nbrtot i = nbrlist(1,m) j = nbrlist(2,m) ni = ntype(i) etai = potparm(6,ni) betai = potparm(4,ni) ri = potparm(2,ni) nj = ntype(j) etaj = potparm(6,nj) betaj = potparm(4,nj) rj = potparm(2,nj) nij = npot(nj,ni) dist = xnbr(0,m) call eamppot(dist,potparm(1,nij),phi,dphi,d2phi) rhoi = -etai*betai*EXP(-betai*(dist-ri)) rhoj = -etaj*betaj*EXP(-betaj*(dist-rj)) f = (dphi+pden(2,i)*rhoj+pden(2,j)*rhoi)/dist do k=1,3 force(k,i) = force(k,i) + f*xnbr(k,m) force(k,j) = force(k,j) - f*xnbr(k,m) enddo epair = epair + phi enddo embed = 0.d0 do i=1,natom embed = embed + pden(1,i) enddo epair = epair + embed return end subroutine espotl(NDIM,natom,nrtot,nktot,nbrtot,nbrlist,xnbr, & Omega,rrvec,rkvec,rnuc,zchg,ntype,hard,chi,zeff,zeta,potparm, & alphew,ees,eew,eembd,epair,etot,vewald,catom,force,pden,vwork, & work,prs) implicit double precision (a-h,o-z) double precision derfc dimension nbrlist(0:2,nbrtot) dimension xnbr(0:3,nbrtot) dimension rrvec(0:3,nrtot) dimension rkvec(0:3,nktot) dimension rnuc(3,natom) dimension force(3,natom) dimension pden(0:3,natom) dimension zchg(natom) dimension ntype(natom) dimension hard(*) dimension chi(*) dimension zeff(*) dimension zeta(*) dimension potparm(10,*) dimension vewald(*),catom(natom) dimension vwork(*) dimension work(*) dimension cpu(0:3) dimension prs(6) call timex(t1) nherm = natom*(natom+1)/2 call escharge(NDIM,natom,nbrtot,nbrlist,xnbr,nrtot,rrvec,nktot, & rkvec,Omega,rnuc,zchg,ntype,hard,chi,zeff,zeta,alphew,ees,eew, & vewald,catom,vwork,work) call timex(t2) cpu(1) = t2 - t1 call AClear(3*natom,force) call AClear(6,prs) call esforce(NDIM,natom,nbrtot,nbrlist,xnbr,nrtot,rrvec,nktot, & rkvec,Omega,rnuc,zchg,ntype,zeff,zeta,alphew,force) call timex(t3) cpu(2) = t3 - t2 call epforce(natom,nbrtot,ntype,nbrlist,xnbr,potparm,epair,eembd, & force,pden) etot = ees + epair call ascale(3*natom,-1.d0,force) call timex(t4) cpu(3) = t4 - t3 cpu(0) = t4 - t1 c print 111, cpu 111 format ('total:',f8.3,' escharge:',f8.3,' esforce:', & f8.3,' epforce:',f8.3) return end subroutine GetDen(natom,nbrtot,nbrlist,xnbr,ntype,potparm,pden) implicit double precision (a-h,o-z) c GetDen: this routine calculates the local atomic densities P_i c within the EAM formalism. dimension nbrlist(0:2,nbrtot) dimension xnbr(0:3,nbrtot) dimension ntype(natom),potparm(10,3) dimension pden(0:3,natom) parameter (TOL=1.d-12) call aClear(4*natom,pden) do m=1,nbrtot i = nbrlist(1,m) j = nbrlist(2,m) n = nbrlist(0,m) ni = ntype(i) etai = potparm(6,ni) betai = potparm(4,ni) ri = potparm(2,ni) nj = ntype(j) etaj = potparm(6,nj) betaj = potparm(4,nj) rj = potparm(2,nj) dist = xnbr(0,m) rhoi = EXP(-betai*(dist-ri)) rhoj = EXP(-betaj*(dist-rj)) pden(0,i) = pden(0,i) + etaj*rhoj pden(0,j) = pden(0,j) + etai*rhoi enddo do i=1,natom ni = ntype(i) call embedpot(potparm(1,ni),pden(0,i)) enddo return end subroutine GetDen2(natom,nrtot,rnuc,rrvec,ntype,potparm,pden) implicit double precision (a-h,o-z) c GetDen: this routine calculates the local atomic densities P_i c within the EAM formalism. dimension rnuc(3,natom) dimension rrvec(0:3,nrtot) dimension ntype(natom),potparm(10,3) dimension pden(0:3,natom) dimension rr(3) parameter (TOL=1.d-12) call aClear(4*natom,pden) do i=1,natom ni = ntype(i) etai = potparm(6,ni) betai = potparm(4,ni) ri = potparm(2,ni) do j=1,i-1 nj = ntype(j) etaj = potparm(6,nj) betaj = potparm(4,nj) rj = potparm(2,nj) dist = 0.d0 do k=1,3 rr(k) = rnuc(k,i) - rnuc(k,j) dist = dist + rr(k)**2 enddo dist = SQRT(dist) rhoi = EXP(-betai*(dist-ri)) rhoj = EXP(-betaj*(dist-rj)) do n=1,nrtot dist1 = 0.d0 dist2 = 0.d0 do k=1,3 dist1 = dist1 + (rr(k)+rrvec(k,n))**2 dist2 = dist2 + (rr(k)-rrvec(k,n))**2 enddo dist1 = SQRT(dist1) dist2 = SQRT(dist2) rhoi = rhoi + EXP(-betai*(dist1-ri)) rhoi = rhoi + EXP(-betai*(dist2-ri)) rhoj = rhoj + EXP(-betaj*(dist1-rj)) rhoj = rhoj + EXP(-betaj*(dist2-rj)) enddo pden(0,i) = pden(0,i) + etaj*rhoj pden(0,j) = pden(0,j) + etai*rhoi enddo rhoi = 0.d0 do n=1,nrtot rhoi = rhoi + EXP(-betai*(rrvec(0,n)-ri)) enddo pden(0,i) = pden(0,i) + 2*etai*rhoi enddo do i=1,natom ni = ntype(i) call embedpot(potparm(1,ni),pden(0,i)) enddo return end subroutine eamppot(dist,potparm,phi,dphi,d2phi) implicit double precision (a-h,o-z) dimension potparm(10) B = potparm(5) C = potparm(7) rstar = potparm(2) alpha = potparm(3) beta = 0.5d0*potparm(4) sig = alpha*(dist-rstar) xp1 = C*EXP(-sig) xp2 = 2*B*EXP(-beta*(dist-rstar)) phi = (xp2-(1.d0+sig)*xp1) dphi = (-beta*xp2+alpha*sig*xp1) d2phi = (beta*beta*xp2+alpha*alpha*(1.d0-sig)*xp1) return end subroutine embedpot(potparm,pden) implicit double precision (a-h,o-z) parameter (TOL=1.d-12) dimension potparm(*) dimension pden(0:3) eta = potparm(6) A = potparm(1)/sqrt(eta) rstar = potparm(2) rho = pden(0) if (rho.gt.TOL) then F1 = -A*SQRT(rho) dF1 = 0.5d0*F1/rho d2F1 = -0.25d0*F1/rho**2 pden(1) = F1 pden(2) = dF1 pden(3) = d2F1 else pden(1) = 0.d0 pden(2) = 0.d0 pden(3) = 0.d0 endif return end subroutine RdPot(nttot,nnuc,chi,hard,zeta,zeff,potparm) implicit double precision (a-h,o-z) dimension nnuc(nttot) dimension chi(nttot) dimension hard(nttot) dimension zeta(nttot) dimension zeff(nttot) dimension potparm(10,3) parameter (MAXSTRING=240) character*(MAXSTRING) potlib character*(MAXSTRING) line character*(MAXSTRING) POTLIBNAME character*(MAXSTRING) label1,label2 integer GtAtNo logical UFLAG,TEST1,TEST2 dimension nz(2) character*1 DELIMIT parameter (DELIMIT=':') parameter (POTLIBNAME='POTLIB') c format of potential.lib file c c atomic potentials: c c atomic-symbol chi(n) hard(n) zeta(n) zeff(n) aembed(n) eta(n) c c Al 0.000000 10.313880 0.934976 0.593190 c c pair potentials c c atom1:atom2 rstar alpha beta B C c Al:O 2.417130 3.425080 3.001020 0.130460 0.130460 c values of potential parameters c c potparm(1,*): A c potparm(2,*): r* c potparm(3,*): alpha c potparm(4,*): beta c potparm(5,*): B c potparm(6,*): eta c potparm(7,*): C c c c use rationalized atomic scale units c c time: femtoseconds [10**(-15) sec] c c distance: Angstrom [10**(-10) m] c c mass: 10**(-28) kg c c energy: 10**(-18) J [6.241808 eV] c c velocity A/fs [10**5 m/sec] c c pressure: 10**(12) Pa (1 TPa, 10 Mbar) c parameter (EtoK=72435.42d0) parameter (EtoeV=6.241808d0) parameter (AMUtoM=16.6053d0) parameter (EtoH=0.229385d0) parameter (AUtoA=0.52917d0) call GetLib(POTLIBNAME,'potential.lib',potlib) nptot = nttot*(nttot+1)/2 do n=1,nttot nnuc(n) = -nnuc(n) enddo do n=1,nptot potparm(2,n) = -1.d0 enddo open (unit=1,file=potlib,status='old',readonly,form='formatted') ntcur = 0 npcur = 0 call getline(1,nline,line) nn = 0 do while ((nline.gt.-1).AND.(ntcur.lt.nttot).AND.(npcur.lt.nptot)) if ((line(1:1).eq.'#').or.(nline.eq.0)) goto 100 ndlm = index(line(1:nline),':') if (ndlm.gt.0) then line(ndlm:ndlm) = ' ' nptype = 2 call GtLbl(nn,nline,line,label1) call SkipWt(nn,nline,line) call GtLbl(nn,nline,line,label2) nz(1) = GtAtNo(label1) nz(2) = GtAtNo(label2) call skipwt(nn,nline,line) call GtReal(nn,nline,line,rstar) call skipwt(nn,nline,line) call GtReal(nn,nline,line,alpha) call skipwt(nn,nline,line) call GtReal(nn,nline,line,beta) call skipwt(nn,nline,line) call GtReal(nn,nline,line,B) B = B/EtoeV call skipwt(nn,nline,line) if (nn.lt.nline) then call GtReal(nn,nline,line,C) C = C/EtoeV else C = B endif else nptype = 1 call GtLbl(nn,nline,line,label1) nz(1) = GtAtNo(label1) call skipwt(nn,nline,line) call GtReal(nn,nline,line,cchi) cchi = cchi/EtoeV call skipwt(nn,nline,line) call GtReal(nn,nline,line,hhard) hhard = hhard/EtoeV call skipwt(nn,nline,line) call GtReal(nn,nline,line,zzeta) call skipwt(nn,nline,line) call GtReal(nn,nline,line,zzeff) call skipwt(nn,nline,line) call GtReal(nn,nline,line,A) A = A/EtoeV call skipwt(nn,nline,line) call GtReal(nn,nline,line,eta) endif c is this line used in the current code? UFLAG = .FALSE. if (nptype.eq.1) then do n=1,nttot TEST1 = (nnuc(n).EQ.-nz(1)) if (TEST1) then nn1 = n UFLAG = .TRUE. endif enddo else do n=1,nttot TEST1 = (nz(1).eq.ABS(nnuc(n))) if (TEST1) then do n2=1,nttot TEST2 = (nz(2).eq.ABS(nnuc(n2))) if (TEST2) then nn1 = n nn2 = n2 if (nn1.eq.nn2) then npos = nn1 else nmin = MIN(nn1,nn2) nmax = MAX(nn1,nn2) npos = nttot + + (nmax-2)*(nmax-1)/2 + nmin endif if (potparm(2,npos).lt.0.d0) then UFLAG = .TRUE. endif endif enddo endif enddo endif c line is used, with atomic numbers given by nn1,nn2 and c npos gives the potential index for pair potentials c now put information in arrays if (UFLAG) then if (nptype.eq.1) then chi(nn1) = cchi hard(nn1) = hhard zeta(nn1) = zzeta zeff(nn1) = zzeff potparm(1,nn1) = A potparm(6,nn1) = eta nnuc(nn1) = ABS(nnuc(nn1)) else potparm(2,npos) = rstar potparm(3,npos) = alpha potparm(4,npos) = beta potparm(5,npos) = B potparm(7,npos) = C endif endif 100 call getline(1,nline,line) nn = 0 enddo close (1) c check to see if all atomic potential parameters have been c read TEST1 = .FALSE. do n=1,nttot if (nnuc(n).lt.0) then call GtAtSy(label1,-nnuc(n)) print *,'RdPot missing atomic potential for ',label1 TEST1 = .TRUE. endif enddo do n=1,nptot if (potparm(2,n).lt.0.d0) then print *,'RdPot missing pair potential: ',n TEST1 = .TRUE. endif enddo if (TEST1) stop return end subroutine esel(NDIM,natom,nbrtot,nbrlist,xnbr,nrtot,rrvec,nktot, & rkvec,Omega,rnuc,zchg,ntype,zeff,zeta,alphew,pvec,xvec,cijkl, & vijkl) implicit double precision (a-h,o-z) double precision derfc c c use rationalized atomic scale units c c time: femtoseconds [10**(-15) sec] c c distance: Angstrom [10**(-10) m] c c mass: 10**(-28) kg c c energy: 10**(-18) J [6.241808 eV] c c velocity A/fs [10**5 m/sec] c c pressure: 10**(13) Pa (10 TPa, 100 Mbar) c parameter (EtoK=72435.42d0) parameter (EtoeV=6.241808d0) parameter (AMUtoM=16.6053d0) parameter (EtoH=0.229387274d0) parameter (AUtoA=0.52917d0) dimension nbrlist(0:2,nbrtot) dimension xnbr(0:3,nbrtot) dimension rrvec(0:3,nrtot) dimension rkvec(0:3,nktot) dimension zchg(natom) dimension rnuc(3,natom),ntype(natom) dimension zeff(*),zeta(*) dimension pvec(6,natom) dimension xvec(6,natom) dimension cijkl(21),vijkl(21) dimension evec(21) dimension e6(6),e21(21) dimension xi6(6),xi21(21) dimension xj6(6),xj21(21) dimension vv6(6),vv21(21) call aClear(21,vijkl) call aClear(21,evec) call aClear(6*natom,xvec) call ewel(NDIM,natom,rnuc,zchg,nrtot,rrvec,nktot,rkvec,Omega, & alphew,pvec,vijkl) do n=1,nbrtot i = nbrlist(1,n) j = nbrlist(2,n) ni = ntype(i) zi = zeta(ni) nj = ntype(j) zj = zeta(nj) do m=1,6 e6(m) = 0.d0 xi6(m) = 0.d0 xj6(m) = 0.d0 vv6(m) = 0.d0 enddo do m=1,21 e21(m) = 0.d0 xi21(m) = 0.d0 xj21(m) = 0.d0 vv21(m) = 0.d0 enddo vsum = 0.d0 dist = xnbr(0,m) call MkEsel(natom,xnbr(0,n),zi,zj,e6,e21,xi6,xi21,xj6,xj21, & vv6,vv21,vsum) do m=1,6 xvec(m,i) = xvec(m,i) + zeff(nj)*xi6(m) xvec(m,j) = xvec(m,j) + zeff(ni)*xj6(m) pvec(m,i) = pvec(m,i) + zchg(j)*vv6(m) pvec(m,j) = pvec(m,j) + zchg(i)*vv6(m) enddo do m=1,21 evec(m) = evec(m) + zeff(nj)*zchg(i)*xi21(m) & + zeff(ni)*zchg(j)*xj21(m) evec(m) = evec(m) + zchg(i)*zchg(j)*vv21(m) enddo enddo do m=1,21 cijkl(m) = vijkl(m) + evec(m) enddo call aScale(6*natom,AUtoA/EtoH,pvec) call aScale(6*natom,AUtoA/EtoH,xvec) do m=1,21 cijkl(m) = cijkl(m)*AUtoA/EtoH vijkl(m) = vijkl(m)*AUtoA/EtoH enddo return end subroutine MkEsel(natom,rr,zi,zj,e6,e21,xi6,xi21,xj6,xj21,vv6, & vv21,vsum) implicit double precision (a-h,o-z) dimension rr(0:3) dimension e6(6),e21(21) dimension xi6(6),xi21(21) dimension xj6(6),xj21(21) dimension vv6(6),vv21(21) dist = rr(0) call CalVij(dist,zi,zj,vij,dvij,d2vij) call CalVn(dist,zi,vi,dvi,d2vi) call CalVn(dist,zj,vj,dvj,d2vj) vsum = vsum + vij de = (dvij-dvi-dvj) d2e = (d2vij-d2vi-d2vj) call MkCijkl(e21,e6,rr,de,d2e) de = (dvi-dvij) d2e = (d2vi-d2vij) call MkCijkl(xi21,xi6,rr,de,d2e) de = (dvj-dvij) d2e = (d2vj-d2vij) call MkCijkl(xj21,xj6,rr,de,d2e) de = dvij d2e = d2vij call MkCijkl(vv21,vv6,rr,de,d2e) return end subroutine ewel(NDIM,natom,rnuc,znuc,nrtot,rrvec,nktot,rkvec, & Omega,alpha,Pvec,vijkl) implicit double precision(a-h,o-z) double precision derfc c calculate Ewald contributions to electrostatic energy c c natom: number of atoms in unit cell c c rnuc(3,natom): nuclear coordinates of atoms c c znuc(natom): nuclear charges c c nrtot: total number of real lattice points c c rrvec(0:3,nrtot): coordinates of lattice sites c c nktot: number of reciprocal lattice points c c rkvec(0:3,nktot): coordinates of reciprocal lattice sites c c Omega: area of unit cell c c alpha: Gaussian exponent for Ewald procedure c dimension rnuc(3,natom), znuc(natom) dimension rrvec(0:3,nrtot) dimension rkvec(0:3,nrtot) dimension Pvec(6,natom), vijkl(21) nherm = natom*(natom+1)/2 call aClear(6*natom,Pvec) call aClear(21,vijkl) call eweld(natom,rnuc,znuc,nrtot,rrvec,alpha,Pvec,vijkl) if (NDIM.eq.2) then call ewel2(natom,rnuc,znuc,nktot,rkvec,Omega,alpha, & Pvec,vijkl) elseif (NDIM.eq.3) then call ewel3(natom,rnuc,znuc,nktot,rkvec,Omega,alpha, & Pvec,vijkl) else print *,'ewel: invalid number of dimensions: ',NDIM stop endif return end subroutine ewel2(natom,rnuc,znuc,nktot,rkvec,Omega,alpha,Pvec, & vijkl) implicit double precision (a-h,o-z) double precision derfc c calculate Ewald contributions from reciprocal lattice c sums to electrostatic pressures c and elastic constants c c natom: number of atoms in unit cell c rnuc(3,natom): nuclear coordinates of atoms c znuc(natom): nuclear charges c nktot: number of reciprocal lattice points c rkvec(0:3,nktot): coordinates of reciprocal lattice sites c Omega: area of unit cell c alpha: Gaussian exponent for Ewald procedure c Pvec(6,natom): pressure components c vijkl(21): elastic constant components dimension rnuc(3,natom),znuc(natom) dimension rkvec(0:3,nktot) dimension rr(0:3) dimension Pvec(6,natom),vijkl(21) dimension cc(21),pp(6) dimension cc2(21),pp2(6) dimension csum(21),psum(6) dimension corig(21),porig(6) dimension volp(6),volpp(21) parameter (TOL=1.d-16) TAU = -DLOG(TOL) c calculate derivatives of inverse volume wrt strains call aClear(6,volp) call aClear(21,volpp) volp(1) = -1.d0 volp(2) = -1.d0 volpp(1) = 2.d0 volpp(2) = 2.d0 - 1.d0 volpp(3) = 2.d0 volpp(10) = -0.25d0 volpp(15) = -0.25d0 volpp(21) = +0.5d0 PI = 3.14159265358979323846 HPI = PI/2.0 TWOPI = 2*PI sqalph = SQRT(ALPHA) y0 = -SQRT(PI/alpha)/Omega alphK = 0.5d0/sqalph sum2 = 0.d0 call aClear(6,porig) call aClear(21,corig) pre = 2*pi/Omega do n=1,nktot rK = rkvec(0,n) ak2 = (alphK*rK)**2 if (ak2.lt.TAU) then xp = EXP(-ak2)/SQRT(alpha*PI) term = derfc(alphK*rk)/rk dFdK = -term/rk - xp/rk d2FdK2 = -2*dFdK/rK + 0.5d0*xp/alpha call MakZK(cc,pp,2*pre*dFdK,2*pre*d2FdK2,0.d0,0.d0,0.d0, & rkvec(0,n),0.d0) sum2 = sum2 + pre*term do m=1,6 porig(m) = porig(m) + pp(m) enddo do m=1,21 corig(m) = corig(m) + cc(m) enddo endif enddo zorig = 2*(y0+sum2) mn = 0 do n=1,6 do m=1,n mn = mn + 1 corig(mn) = corig(mn) + zorig*volpp(mn) + porig(m)*volp(n) & + porig(n)*volp(m) enddo enddo porig(1) = porig(1) - zorig porig(2) = porig(2) - zorig ij = 0 do i=1,natom do j=1,i-1 call aClear(6,psum) call aClear(21,csum) dist = 0.d0 do k=1,3 rr(k) = rnuc(k,i) - rnuc(k,j) dist = dist + rr(k)**2 enddo dist = SQRT(dist) zmn = rr(3) azmn = sqalph*zmn xp1 = 1.d0 - derfc(azmn) xp2 = EXP(-azmn*azmn) xp3 = 4*SQRT(PI*alpha)*xp2/Omega term1 = 2*pi*zmn*xp1/Omega term2 = 2*SQRT(pi/alpha)*xp2/Omega sum0 = -(term1+term2) dFdZ = 2*PI*xp1/Omega + zmn*xp3 - 2*alpha*zmn*term2 d2FdZ2 = 2*(1.d0-azmn*azmn)*xp3 + & ((2*alpha*zmn)**2-2*alpha)*term2 psum(3) = psum(3) - zmn*dFdZ csum(6) = csum(6) - d2FdZ2*zmn*zmn csum(10) = csum(10) + 0.75d0*zmn*dFdZ csum(15) = csum(15) + 0.75d0*zmn*dFdZ sum2 = 0.d0 do n=nktot,1,-1 rK = rkvec(0,n) if (ABS(zmn*rK).lt.TAU) then aa = alphK*rK pre = 2*PI/Omega xp0 = EXP(zmn*rK) f1 = pre*xp0*derfc(aa+azmn)/rK f2 = pre*derfc(aa-azmn)/xp0/rK xp = pre*EXP(-azmn*azmn-aa*aa)/SQRT(alpha*PI) df1dK = zmn*f1 - (f1+xp)/rK df2dK = -zmn*f2 - (f2+xp)/rK dFdK = df1dK + df2dK d2f1dK2 = (zmn-1.d0/rK)*df1dK + f1/rK**2 + & (0.5d0/alpha+1.d0/rK**2)*xp d2f2dK2 = (-zmn-1.d0/rK)*df2dK + f2/rK**2 + & (0.5d0/alpha+1.d0/rK**2)*xp d2FdK2 = d2f1dK2 + d2f2dK2 df1dZ = rK*f1 - 2*alpha*xp/rK df2dZ = -rK*f2 + 2*alpha*xp/rK dFdZ = df1dZ + df2dZ d2f1dZ2 = rK*df1dZ + 4*alpha*alpha*zmn*xp/rK d2f2dZ2 = -rK*df2dZ - 4*alpha*alpha*zmn*xp/rK d2FdZ2 = d2f1dZ2 + d2f2dZ2 d2f1dZdK = f1 + rK*df1dK + xp*(1.d0+0.5d0/aa**2) d2f2dZdK = -f2 - rK*df2dK - xp*(1.d0+0.5d0/aa**2) d2FdZdK = d2f1dZdK + d2f2dZdK call MakZK(cc,pp,dFdK,d2FdK2,dFdZ,d2FdZ2,d2FdZdK, & rkvec(0,n),zmn) c include derivatives of phase factor wrt strain phase = rkvec(1,n)*rr(1) + rkvec(2,n)*rr(2) pcos = cos(phase) psin = cos(phase-HPI) eterm = f1 + f2 call MakPhase(cc2,pp2,rkvec(0,n),rr) mn = 0 do mm=1,6 do nn=1,mm mn = mn + 1 cc(mn) = pcos*(cc(mn)-eterm*pp2(mm)*pp2(nn)) & - psin*(pp(mm)*pp2(nn)+pp(nn)*pp2(mm) & +eterm*cc2(mn)) enddo enddo do m=1,6 pp(m) = pcos*pp(m) - psin*eterm*pp2(m) enddo do m=1,21 csum(m) = csum(m) + cc(m) enddo do m=1,6 psum(m) = psum(m) + pp(m) enddo sum2 = sum2 + eterm*pcos endif enddo c include derivatives of volume wrt strain eterm = sum2 + sum0 mn = 0 do n=1,6 do m=1,n mn = mn + 1 csum(mn) = csum(mn) + eterm*volpp(mn) & + psum(m)*volp(n) + psum(n)*volp(m) enddo enddo psum(1) = psum(1) - eterm psum(2) = psum(2) - eterm do m=1,6 Pvec(m,i) = Pvec(m,i) + znuc(j)*psum(m) Pvec(m,j) = Pvec(m,j) + znuc(i)*psum(m) enddo do m=1,21 vijkl(m) = vijkl(m) + znuc(i)*znuc(j)*csum(m) enddo enddo eterm = zorig do m=1,6 Pvec(m,i) = Pvec(m,i) + znuc(i)*porig(m) enddo do m=1,21 vijkl(m) = vijkl(m) + 0.5d0*corig(m)*znuc(i)**2 enddo enddo return end subroutine MakPhase(cc,pp,kk,rr) implicit double precision (a-h,o-z) dimension cc(21) dimension pp(6) double precision kk(0:3) dimension rr(0:3) pp(1) = 0.d0 pp(2) = 0.d0 pp(3) = 0.d0 pp(4) = kk(2)*rr(3) pp(5) = kk(1)*rr(3) pp(6) = 0.d0 cc(1) = 0.d0 cc(2) = 0.d0 cc(3) = 0.d0 cc(4) = 0.d0 cc(5) = 0.d0 cc(6) = 0.d0 cc(7) = 0.d0 cc(8) = -1.5d0*kk(2)*rr(3) cc(9) = 0.5d0*kk(2)*rr(3) cc(10) = 0.d0 cc(11) = -1.5d0*kk(1)*rr(3) cc(12) = 0.d0 cc(13) = 0.5d0*kk(1)*rr(3) cc(14) = 0.d0 cc(15) = 0.d0 cc(16) = 0.d0 cc(17) = 0.d0 cc(18) = 0.d0 cc(19) = -0.75d0*kk(1)*rr(3) cc(20) = -0.75d0*kk(2)*rr(3) cc(21) = 0.d0 return end subroutine MakZK(cc,pp,dFdK,d2FdK2,dFdZ,d2FdZ2,d2FdZdK,rr,zmn) implicit double precision (a-h,o-z) dimension cc(21),pp(6) dimension rr(0:3) c c Voigt notation vv1 = rr(1)*rr(1) vv2 = rr(2)*rr(2) vv6 = rr(1)*rr(2) c Elastic constants pterm = dFdK/rr(0) cterm = (d2FdK2-pterm)/rr(0)**2 do m=1,21 cc(m) = 0.d0 enddo cc(1) = vv1*vv1*cterm cc(2) = vv1*vv2*cterm cc(3) = vv2*vv2*cterm cc(16) = vv1*vv6*cterm cc(17) = vv2*vv6*cterm cc(21) = vv6*vv6*cterm pp(1) = -pterm*vv1 pp(2) = -pterm*vv2 pp(3) = 0.d0 pp(4) = 0.d0 pp(5) = 0.d0 pp(6) = -pterm*vv6 c do i=1,6 c pp(i) = -pterm*Voigt(i) c enddo c Add delta function terms to elastic constants cc(1) = cc(1) - 3*pp(1) cc(3) = cc(3) - 3*pp(2) cc(10) = cc(10) + 0.25d0*(pp(2)+pp(3)) cc(14) = cc(14) + 0.25d0*pp(6) cc(15) = cc(15) + 0.25d0*pp(1) cc(16) = cc(16) - 1.5d0*pp(6) cc(17) = cc(17) - 1.5d0*pp(6) cc(21) = cc(21) - 0.75d0*(pp(1)+pp(2)) c add contributions with respect to Z coordinate pp(3) = pp(3) + zmn*dFdZ d2 = d2FdZdK/rr(0) cc(4) = cc(4) - rr(1)*rr(1)*zmn*d2 cc(5) = cc(5) - rr(2)*rr(2)*zmn*d2 cc(6) = cc(6) - 2*rr(3)*rr(3)*zmn*d2 + d2FdZ2*zmn*zmn cc(9) = cc(9) - rr(2)*rr(3)*zmn*d2 cc(10) = cc(10) - 0.75d0*zmn*dFdZ cc(13) = cc(13) - rr(1)*rr(3)*zmn*d2 cc(15) = cc(15) - 0.75d0*zmn*dFdZ cc(18) = cc(18) - rr(1)*rr(2)*zmn*d2 return end subroutine ewel3(natom,rnuc,znuc,nktot,rkvec,Omega,alpha,Pvec, & vijkl) implicit double precision (a-h,o-z) double precision derfc c calculate Ewald contributions from reciprocal lattice c sums to electrostatic pressures c and elastic constants c c natom: number of atoms in unit cell c rnuc(3,natom): nuclear coordinates of atoms c znuc(natom): nuclear charges c nktot: number of reciprocal lattice points c rkvec(0:3,nktot): coordinates of reciprocal lattice sites c Omega: area of unit cell c alpha: Gaussian exponent for Ewald procedure c Pvec(6,natom): pressure components c vijkl(21): elastic constant components c common /stest/vtest(21,100),ptest(6,100),vstore(100) c save /stest/ dimension rnuc(3,natom),znuc(natom) dimension rkvec(0:3,nktot) dimension Pvec(6,natom),vijkl(21) dimension cc(21),pp(6) dimension csum(21),psum(6) dimension volp(6),volpp(21) PI = 3.14159265358979323846 alphk = 0.25d0/alpha call aClear(6,psum) call aClear(21,csum) pre = 8*PI/Omega c calculate derivatives of inverse volume wrt strains call aClear(6,volp) call aClear(21,volpp) vv = -1.d0 vv2 = 2.d0 volp(1) = vv volp(2) = vv volp(3) = vv volpp(1) = vv2 volpp(2) = vv + vv2 volpp(3) = vv2 volpp(4) = vv + vv2 volpp(5) = vv + vv2 volpp(6) = vv2 volpp(10) = -vv/2.d0 volpp(15) = -vv/2.d0 volpp(21) = -vv/2.d0 zsum = 0.d0 do n=nktot,1,-1 d2 = rkvec(0,n)**2 add = alphk*d2 v = pre*EXP(-add)/d2 vd = -2*v*(1.d0+add)/rkvec(0,n) vdd = 2*v*(3+add*(3+2*add))/d2 call MkKijkl(cc,pp,rkvec(0,n),vd,vdd) zsum = zsum + v mm = 0 do m1=1,6 do m2=1,m1 mm = mm + 1 cc(mm) = cc(mm) + v*volpp(mm) & + (pp(m1)*volp(m2)+pp(m2)*volp(m1)) csum(mm) = csum(mm) + cc(mm) enddo enddo do m=1,6 pp(m) = +pp(m) + volp(m)*v psum(m) = psum(m) + pp(m) enddo ij = 0 do i=1,natom do j=1,i-1 ij = ij + 1 phase = 0.d0 do k=1,3 phase = phase + rkvec(k,n)*(rnuc(k,i)-rnuc(k,j)) enddo cpre = COS(phase) do m=1,6 Pvec(m,i) = Pvec(m,i) + znuc(j)*cpre*pp(m) Pvec(m,j) = Pvec(m,j) + znuc(i)*cpre*pp(m) enddo do m=1,21 vijkl(m) = vijkl(m) + znuc(i)*znuc(j)*cpre*cc(m) enddo c do m=1,21 c vtest(m,ij) = vtest(m,ij) + cpre*cc(m) c enddo c do m=1,6 c ptest(m,ij) = ptest(m,ij) + cpre*pp(m) c enddo enddo ij = ij + 1 enddo enddo ij = 0 do i=1,natom ij = ij + i do m=1,6 Pvec(m,i) = Pvec(m,i) + znuc(i)*psum(m) enddo do m=1,21 vijkl(m) = vijkl(m) + 0.5d0*csum(m)*znuc(i)**2 enddo c do m=1,21 c vtest(m,ij) = vtest(m,ij) + csum(m) c enddo c do m=1,6 c ptest(m,ij) = ptest(m,ij) + psum(m) c enddo enddo return end subroutine MkKijkl(Cvec,Pvec,Rvec,vd,vdd) implicit double precision (a-h,o-z) c MkKijkl calculates pair potential contribution to pressures c and elastic constants wrt pair potential in K c c Pressure notation c c P(1) 11 c P(2) 22 c P(3) 33 c P(4) 23 c P(5) 13 c P(6) 12 c c Elastic constant notation c c C(1) 11 11 c C(2) 11 22 c C(3) 22 22 c C(4) 11 33 c C(5) 22 33 c C(6) 33 33 c C(7) 11 23 c C(8) 22 23 c C(9) 33 23 c C(10) 23 23 c C(11) 11 13 c C(12) 22 13 c C(13) 33 13 c C(14) 23 13 c C(15) 13 13 c C(16) 11 12 c C(17) 22 12 c C(18) 33 12 c C(19) 23 12 c C(20) 13 12 c C(21) 12 12 dimension Cvec(21) dimension Pvec(6) dimension Rvec(0:3) dimension Voigt(6) c Voigt notation Voigt(1) = Rvec(1)*Rvec(1) Voigt(2) = Rvec(2)*Rvec(2) Voigt(3) = Rvec(3)*Rvec(3) Voigt(4) = Rvec(2)*Rvec(3) Voigt(5) = Rvec(1)*Rvec(3) Voigt(6) = Rvec(1)*Rvec(2) c Elastic constants pterm = vd/Rvec(0) cterm = (vdd-pterm)/Rvec(0)**2 ij = 0 do i=1,6 do j=1,i ij = ij + 1 Cvec(ij) = Voigt(i)*Voigt(j)*cterm enddo enddo do i=1,6 Pvec(i) = -pterm*Voigt(i) enddo c Add delta function terms to elastic constants Cvec(1) = Cvec(1) - 3*Pvec(1) Cvec(3) = Cvec(3) - 3*Pvec(2) Cvec(6) = Cvec(6) - 3*Pvec(3) Cvec(8) = Cvec(8) - 1.5d0*Pvec(4) Cvec(9) = Cvec(9) - 1.5d0*Pvec(4) Cvec(10) = Cvec(10) - 0.75d0*(Pvec(2)+Pvec(3)) Cvec(11) = Cvec(11) - 1.5d0*Pvec(5)/2.d0 Cvec(13) = Cvec(13) - 1.5d0*Pvec(5)/2.d0 Cvec(14) = Cvec(14) - 0.75d0*Pvec(6) Cvec(15) = Cvec(15) - 0.75d0*(Pvec(1)+Pvec(3)) Cvec(16) = Cvec(16) - 1.5d0*Pvec(6) Cvec(17) = Cvec(17) - 1.5d0*Pvec(6) Cvec(19) = Cvec(19) - 1.5d0*Pvec(5) Cvec(20) = Cvec(20) - 1.5d0*Pvec(4) Cvec(21) = Cvec(21) - 0.75d0*(Pvec(1)+Pvec(2)) return end subroutine eweld(natom,rnuc,znuc,nrtot,rrvec,alpha,Pvec,vijkl) implicit double precision (a-h,o-z) double precision derfc c calculate direct Ewald contributions to electrostatic pressure c and elastic constants c c natom: number of atoms in unit cell c rnuc(3,natom): nuclear coordinates of atoms c znuc(natom): nuclear charges c nrtot: total number of real lattice points c rrvec(0:3,nrtot): coordinates of lattice sites c alpha: Gaussian exponent for Ewald procedure c Pvec(natom,6): pressure components c vijkl(21): elastic constant components dimension rnuc(3,natom),znuc(natom) dimension rrvec(0:3,nrtot) dimension rr(0:3),rr1(0:3),rr2(0:3) dimension Pvec(6,natom),vijkl(21) dimension pp(6),cc(21) parameter (TOL=1.d-16) TAU = -DLOG(TOL) PI = 3.14159265358979323846 sqalph = SQRT(alpha) sqa2pi = SQRT(alpha/PI) do i=1,natom do j=1,i-1 call aClear(6,pp) call aClear(21,cc) dd = 0.d0 do k=1,3 rr(k) = rnuc(k,i) - rnuc(k,j) dd = dd + rr(k)**2 enddo dist = SQRT(dd) rr(0) = dist do n=nrtot,1,-1 dd1 = 0.d0 dd2 = 0.d0 do k=1,3 rr1(k) = rr(k) - rrvec(k,n) rr2(k) = rr(k) + rrvec(k,n) dd1 = dd1 + rr1(k)**2 dd2 = dd2 + rr2(k)**2 enddo if (alpha*dd1.lt.TAU) then dist1 = SQRT(dd1) rr1(0) = dist1 f1 = derfc(sqalph*dist1)/dd1 g1 = 2*sqa2pi*EXP(-alpha*dd1) v1d = -f1 - g1/dist1 v1dd = 2*f1/dist1 + 2*g1*(alpha+1.d0/dd1) call MkCijkl(cc,pp,rr1,v1d,v1dd) endif if (alpha*dd2.lt.TAU) then dist2 = SQRT(dd2) rr2(0) = dist2 f2 = derfc(sqalph*dist2)/dd2 g2 = 2*sqa2pi*EXP(-alpha*dd2) v2d = -f2 - g2/dist2 v2dd = 2*f2/dist2 + 2*g2*(alpha+1.d0/dd2) call MkCijkl(cc,pp,rr2,v2d,v2dd) endif enddo f = derfc(sqalph*dist)/dd g = 2*sqa2pi*EXP(-alpha*dd) vd = -f - g/dist vdd = 2*f/dist + 2*g*(alpha+1.d0/dd) call MkCijkl(cc,pp,rr,vd,vdd) c add to current vectors do m=1,21 vijkl(m) = vijkl(m) + cc(m)*znuc(i)*znuc(j) enddo do m=1,6 Pvec(m,i) = Pvec(m,i) + pp(m)*znuc(j) Pvec(m,j) = Pvec(m,j) + pp(m)*znuc(i) enddo enddo enddo call aClear(6,pp) call aClear(21,cc) fsum = 0.d0 do n=nrtot,1,-1 dist = rrvec(0,n) dd = dist*dist if (alpha*dd.lt.TAU) then f = derfc(sqalph*dist)/dd g = 2*sqa2pi*EXP(-alpha*dd) vd = -f - g/dist vdd = 2*f/dist + 2*g*(alpha+1.d0/dd) call MkCijkl(cc,pp,rrvec(0,n),vd,vdd) fsum = fsum + f*dist endif enddo do i=1,natom do m=1,6 Pvec(m,i) = Pvec(m,i) + 2*pp(m)*znuc(i) enddo do m=1,21 vijkl(m) = vijkl(m) + cc(m)*znuc(i)**2 enddo enddo return end subroutine MkCijkl(Cijkl,Pij,Rvec,vd,vdd) implicit double precision (a-h,o-z) c MkCijkl calculates pair potential contribution to stresses c and elastic constants c c Pressure notation c c P(1) 11 c P(2) 22 c P(3) 33 c P(4) 23 c P(5) 13 c P(6) 12 c c Elastic constant notation c c C(1) 11 11 c C(2) 11 22 c C(3) 22 22 c C(4) 11 33 c C(5) 22 33 c C(6) 33 33 c C(7) 11 23 c C(8) 22 23 c C(9) 33 23 c C(10) 23 23 c C(11) 11 13 c C(12) 22 13 c C(13) 33 13 c C(14) 23 13 c C(15) 13 13 c C(16) 11 12 c C(17) 22 12 c C(18) 33 12 c C(19) 23 12 c C(20) 13 12 c C(21) 12 12 dimension Cvec(21),Cijkl(21) dimension Pvec(6),Pij(6) dimension Rvec(0:3) dimension Voigt(6) c Voigt notation Voigt(1) = Rvec(1)*Rvec(1) Voigt(2) = Rvec(2)*Rvec(2) Voigt(3) = Rvec(3)*Rvec(3) Voigt(4) = Rvec(2)*Rvec(3) Voigt(5) = Rvec(1)*Rvec(3) Voigt(6) = Rvec(1)*Rvec(2) c Pressures pterm = vd/Rvec(0) do i=1,6 Pvec(i) = pterm*Voigt(i) enddo do i=1,6 Pij(i) = Pij(i) + Pvec(i) enddo c Elastic constants cterm = (vdd-pterm)/Rvec(0)**2 ij = 0 do i=1,6 do j=1,i ij = ij + 1 Cvec(ij) = Voigt(i)*Voigt(j)*cterm enddo enddo c Add delta function terms to elastic constants Cvec(1) = Cvec(1) + Pvec(1) Cvec(3) = Cvec(3) + Pvec(2) Cvec(6) = Cvec(6) + Pvec(3) Cvec(8) = Cvec(8) + Pvec(4)/2.d0 Cvec(9) = Cvec(9) + Pvec(4)/2.d0 Cvec(10) = Cvec(10) + (Pvec(2)+Pvec(3))/4.d0 Cvec(11) = Cvec(11) + Pvec(5)/2.d0 Cvec(13) = Cvec(13) + Pvec(5)/2.d0 Cvec(14) = Cvec(14) + Pvec(6)/4.d0 Cvec(15) = Cvec(15) + (Pvec(1)+Pvec(3))/4.d0 Cvec(16) = Cvec(16) + Pvec(6)/2.d0 Cvec(17) = Cvec(17) + Pvec(6)/2.d0 Cvec(19) = Cvec(19) + Pvec(5)/4.d0 Cvec(20) = Cvec(20) + Pvec(4)/4.d0 Cvec(21) = Cvec(21) + (Pvec(1)+Pvec(2))/4.d0 do i=1,21 Cijkl(i) = Cijkl(i) + Cvec(i) enddo return end subroutine elpair(natom,nptot,ntype,nbrtot,nbrlist,xnbr,potparm, & ppair,cpair,pden,eamterm) implicit double precision (a-h,o-z) dimension ntype(natom) dimension nbrlist(0:2,nbrtot) dimension xnbr(0:3,nbrtot) dimension potparm(10,3) dimension ppair(6) dimension cpair(21) dimension pden(0:3,natom) dimension eamterm(6,natom) dimension voigt(6) dimension npot(2,2) data npot/1,3,3,2/ call aClear(6,ppair) call aClear(21,cpair) call aClear(6*natom,eamterm) call GetDen(natom,nbrtot,nbrlist,xnbr,ntype,potparm,pden) do m=1,nbrtot i = nbrlist(1,m) j = nbrlist(2,m) ni = ntype(i) etai = potparm(6,ni) betai = potparm(4,ni) ri = potparm(2,ni) nj = ntype(j) etaj = potparm(6,nj) betaj = potparm(4,nj) rj = potparm(2,nj) nij = npot(nj,ni) dist = xnbr(0,m) call eamppot(dist,potparm(1,nij),phi,dphi,d2phi) rhoi = -etai*betai*EXP(-betai*(dist-ri)) rhoj = -etaj*betaj*EXP(-betaj*(dist-rj)) rho2i = -betai*rhoi rho2j = -betaj*rhoj de = dphi + rhoj*pden(2,i) + rhoi*pden(2,j) d2e = d2phi + rho2j*pden(2,i) + rho2i*pden(2,j) Voigt(1) = xnbr(1,m)*xnbr(1,m) Voigt(2) = xnbr(2,m)*xnbr(2,m) Voigt(3) = xnbr(3,m)*xnbr(3,m) Voigt(4) = xnbr(2,m)*xnbr(3,m) Voigt(5) = xnbr(1,m)*xnbr(3,m) Voigt(6) = xnbr(1,m)*xnbr(2,m) pj = rhoj/xnbr(0,m) pi = rhoi/xnbr(0,m) do k=1,6 eamterm(k,i) = eamterm(k,i) + pj*Voigt(k) eamterm(k,j) = eamterm(k,j) + pi*Voigt(k) enddo call MkCijkl(cpair,ppair,xnbr(0,m),de,d2e) enddo c add unique EAM terms mn = 0 do m=1,6 do n=1,m mn = mn + 1 sum = 0.d0 do i=1,natom sum = sum + pden(3,i)*eamterm(m,i)*eamterm(n,i) enddo cpair(mn) = cpair(mn) + sum enddo enddo return end subroutine MkPij(Pvec,Rvec,vd) implicit double precision (a-h,o-z) c MkPij calculates pair potential contribution to stresses c and elastic constants c c Pressure notation c c P(1) 11 c P(2) 22 c P(3) 33 c P(4) 23 c P(5) 13 c P(6) 12 c dimension Pvec(6) dimension Rvec(0:3) dimension Voigt(6) c Voigt notation Voigt(1) = Rvec(1)*Rvec(1) Voigt(2) = Rvec(2)*Rvec(2) Voigt(3) = Rvec(3)*Rvec(3) Voigt(4) = Rvec(2)*Rvec(3) Voigt(5) = Rvec(1)*Rvec(3) Voigt(6) = Rvec(1)*Rvec(2) c Elastic constants pterm = vd/Rvec(0) do i=1,6 Pvec(i) = pvec(i) + pterm*Voigt(i) enddo return end subroutine update(natom,nrtot,nbrtot,nbrlist,rnuc,rrvec,xnbr) implicit double precision (a-h,o-z) dimension nbrlist(0:2,nbrtot) dimension rnuc(3,natom) dimension rrvec(0:3,nrtot) dimension xnbr(0:3,nbrtot) c call timex(tstart) do m=1,nbrtot n = nbrlist(0,m) if (n.eq.0) then do k=1,3 xnbr(k,m) = rnuc(k,nbrlist(1,m)) & - rnuc(k,nbrlist(2,m)) enddo else do k=1,3 xnbr(k,m) = rnuc(k,nbrlist(1,m)) & - rnuc(k,nbrlist(2,m)) xnbr(k,m) = xnbr(k,m) - rrvec(k,n) enddo endif enddo do m=1,nbrtot dist = 0.d0 do k=1,3 dist = dist + xnbr(k,m)**2 enddo xnbr(0,m) = SQRT(dist) enddo #ifdef DEBUG rmin=xnbr(0,1) imin=1 do m = 2, nbrtot if(xnbr(0,m).lt.rmin) then rmin=xnbr(0,m) imin=m endif enddo c call timex(tend) print 10000,nbrtot,rmin,nbrlist(0,imin),nbrlist(1,imin),nbrlist(2,imin) #endif call flush(6) 10000 format('update: ',i10,f10.4,3i6) return end subroutine ascale(n,s,vector) c c MULTIPLIES A VECTOR ARRAY TIME A SCALAR S c implicit double precision (a-h,o-z) dimension vector(n) do k=1,n vector(k) = vector(k)*s enddo return end subroutine aclr(n,vector) c c SETS VECTOR = 0 c implicit double precision (a-h,o-z) dimension vector(n) do k=1,n vector(k) = 0.d0 enddo return end subroutine axpy(n,vector1,vector2,s) c c SUMS SCALAR S TIMES VECTOR2 INTO VECTOR1 c implicit double precision (a-h,o-z) dimension vector1(n),vector2(n) do k=1,n vector1(k) = s*vector2(k) + vector1(k) enddo return end subroutine across(a,b,c) c c CALCULATES CROSS PRODUCT A = B X C c implicit double precision (a-h,o-z) dimension a(3),b(3),c(3) a(1) = b(2)*c(3) - b(3)*c(2) a(2) = b(3)*c(1) - b(1)*c(3) a(3) = b(1)*c(2) - b(2)*c(1) return end subroutine agemv(n,array,vector1,vector2) c c STORES ARRAY TIMES VECTOR1 IN VECTOR 2 c implicit double precision (a-h,o-z) dimension array(n,n),vector1(n),vector2(n) do k=1,n vector2(k) = 0.d0 do j=1,n vector2(k) = array(k,j)*vector1(j) + vector2(k) enddo enddo return end subroutine agemm(m,n,a,b,c) c c performs matrix multiplication c= a x b c implicit double precision (a-h,o-z) dimension a(m,n),b(n,m),c(m,m) do k=1,m do j=1,m do i=1,n c(j,k) = c(j,k) + a(j,i)*b(i,k) enddo enddo enddo return end double precision function anrm2(n,vector,nstep) c c takes norm of vector, assumes nstep=1 c implicit double precision (a-h,o-z) dimension vector(n) anrm2 = sqrt(adot(n,vector,vector)) return end double precision function adot(n,vector1,vector2) c c takes dot product of vector1 and vector2 c implicit double precision (a-h,o-z) dimension vector1(n),vector2(n) sum = 0.d0 do j=1,n sum = vector1(j)*vector2(j) + sum enddo adot = sum return end subroutine acopy(n,vector1,nstep,vector2) c c copies vector1 (using stride nstep) into vector2 c implicit double precision (a-h,o-z) dimension vector1(*),vector2(n) k = 1 do j=1,n k = (j-1)*nstep + 1 vector2(j) = vector1(k) enddo return end subroutine acrs(a,b,c) implicit double precision (a-h,o-z) dimension c(3),a(3),b(3) a(1) = b(2)*c(3) - b(3)*c(2) a(2) = b(3)*c(1) - b(1)*c(3) a(3) = b(1)*c(2) - b(2)*c(1) return end subroutine strainedge(natom,nbot,ntop,rnuc,dstrain) implicit double precision (a-h,o-z) dimension rnuc(3,natom),dstrain(3),scale(3) rtop = 0.d0 rbot = 0.d0 do k=1,nbot rbot = rbot + rnuc(3,k) rnuc(1,k) = rnuc(1,k) - 0.5d0*dstrain(1) rnuc(2,k) = rnuc(2,k) - 0.5d0*dstrain(2) rnuc(3,k) = rnuc(3,k) - 0.5d0*dstrain(3) enddo do k=natom-ntop+1,natom rtop = rtop + rnuc(3,k) rnuc(1,k) = rnuc(1,k) + 0.5d0*dstrain(1) rnuc(2,k) = rnuc(2,k) + 0.5d0*dstrain(2) rnuc(3,k) = rnuc(3,k) + 0.5d0*dstrain(3) enddo rbot = rbot/float(nbot) rtop = rtop/float(ntop) scale(1) = dstrain(1)/(rtop-rbot) scale(2) = dstrain(2)/(rtop-rbot) scale(3) = dstrain(3)/(rtop-rbot) origin = (rtop+rbot)/2.d0 do k=nbot+1,natom-ntop rnuc(1,k) = rnuc(1,k) + scale(1)*(rnuc(3,k)-origin) rnuc(2,k) = rnuc(2,k) + scale(2)*(rnuc(3,k)-origin) rnuc(3,k) = rnuc(3,k) + scale(3)*(rnuc(3,k)-origin) enddo return end subroutine shearpull(natom,nbot,ntop,vnuc,halfvel) implicit double precision (a-h,o-z) dimension vnuc(3,natom),halfvel(3) c c NOTE THAT THE Z COMPONENT OF THE VELOCITY IS NOT CHANGED ! c do k=1,nbot vnuc(1,k) = - halfvel(1) vnuc(2,k) = - halfvel(2) enddo do k=natom-ntop+1,natom vnuc(1,k) = + halfvel(1) vnuc(2,k) = + halfvel(2) enddo return end subroutine zeroforce(natom,nbot,ntop,force,dif,add) implicit double precision (a-h,o-z) dimension force(3,natom) dimension dif(3),add(3) do i=1,3 bot = 0.d0 top = 0.d0 do k=1,nbot bot = bot + force(i,k) force(i,k) = 0.d0 enddo do k=natom-ntop+1,natom top = top + force(i,k) force(i,k) = 0.d0 enddo dif(i) = bot add(i) = top enddo return end subroutine assignforce(natom,nbot,ntop,force,halfzforce,sum1,sum2) implicit double precision (a-h,o-z) dimension force(3,natom) do k=1,nbot force(3,k) = halfzforce enddo do k=natom-ntop+1,natom force(3,k) = -halfzforce enddo return end subroutine align(NDIM,natom,ucell,kcell,rnuc) implicit double precision (a-h,o-z) c align: aligns nuclear coordinates inside unit cell c c NDIM: number of periodic boundary condiations c natom: number of atoms c ucell: PBC's c rnuc: nuclear coordinates double precision ucell(3,NDIM), kcell(3,NDIM) dimension rnuc(3,natom) TWOPI = 2*DACOS(-1.d0) do k=1,NDIM do i=1,natom ravg = 0.d0 do m=1,3 ravg = ravg + kcell(m,k)*rnuc(m,i) enddo nshift = NINT(ravg/TWOPI) if (nshift.ne.0) then print *, 'shifting atom',i,' by ', nshift do m=1,3 rnuc(m,i) = rnuc(m,i) - nshift*ucell(m,k) enddo endif enddo enddo return end subroutine applystress(p,r,h,n) implicit double precision (a-h,o-z) dimension p(3,3),h(3,3),hinv(3,3),hnew(3,3) dimension r(3,n),s(3) call matinv(h,hinv,det) call agemm(3,3,p,h,hnew) do i = 1, n call agemv(3,hinv,r(1,i),s) call agemv(3,hnew,s,r(1,i)) enddo call acopy(9,hnew,h) return end SUBROUTINE SORT3(N,RA,RB,RC) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION RA(N) DIMENSION RB(N) DIMENSION RC(N) L = N/2 + 1 IR = N 100 CONTINUE IF (L.GT.1) THEN L = L - 1 RRA = RA(L) RRB = RB(L) RRC = RC(L) ELSE RRA = RA(IR) RRB = RB(IR) RRC = RC(IR) RA(IR) = RA(1) RB(IR) = RB(1) RC(IR) = RC(1) IR = IR - 1 IF (IR.EQ.1) THEN RA(1) = RRA RB(1) = RRB RC(1) = RRC RETURN ENDIF ENDIF I = L J = L + L 200 IF (J.LE.IR) THEN IF (J.LT.IR) THEN IF (RA(J).LT.RA(J+1)) J = J + 1 ENDIF IF (RRA.LT.RA(J)) THEN RA(I) = RA(J) RB(I) = RB(J) RC(I) = RC(J) I = J J = J + J ELSE J = IR + 1 ENDIF GOTO 200 ENDIF RA(I) = RRA RB(I) = RRB RC(I) = RRC GOTO 100 END subroutine GetConfig(jconf,nstep,tprint,itemp,jtest,istyle,pmass, & nconf,ilist,iconf,dt,teq,gamma,nbot,ntop,dstrain,intstrn,pq, & papp,hv,hf) implicit double precision (a-h,o-z) dimension papp(3,3),ppack(6),hv(3),dstrain(3) character*72 title open (unit=4,file='config',status='old') read (4,10000) title read (4,*) nstep,tprint,itemp,jtest,istyle,pmass,nconf do j=1,jconf if (istyle.eq.1) then read (4,*) iconf,dt,teq,gamma,ilist else if (istyle.eq.2) then read (4,*) iconf,dt,teq,gamma,nbot,ntop,dstrain,intstrn else if (istyle.eq.3) then read (4,*) iconf,dt,teq,pq,gamma if (pq.eq.-99.0) then read (4,*) ppack peq = ppack(1) + ppack(3) + ppack(6) papp(1,1)=ppack(1)+1.0 papp(2,1)=ppack(6)/2.0 papp(3,1)=ppack(5)/2.0 papp(1,2)=ppack(6)/2.0 papp(2,2)=ppack(2)+1.0 papp(2,3)=ppack(4)/2.0 papp(3,1)=ppack(5)/2.0 papp(3,2)=ppack(4)/2.0 papp(3,3)=ppack(3)+1.0 else peq = pq endif else if (istyle.eq.4) then read (4,*) iconf,dt,teq,gamma,nbot,ntop,hv,hf,ilist hv(1)=hv(1)*0.5 hv(2)=hv(2)*0.5 hv(3)=hv(3)*0.5 hf=hf/20000.0 else print *,'incorrect istyle entry: ',istyle stop endif enddo close (4) c c PRINT HEADERS IN FILES... c if (jconf.eq.1) then open (1,file='zoutput',status='unknown') open (20,file='zhmat',status='unknown') open (21,file='zstrain',status='unknown') open (24,file='zmod',status='unknown') write (1,10000) title if (itemp.eq.1) write (1,*) 'Velocity scaling' if (itemp.eq.2) write (1,*) 'Pseudo-steepest descent' if (itemp.eq.3) write (1,*) & 'Pseudo-steepest descent and Velocity scaling' if (jtest.ne.0) then write (1,*) 'Dynamic equilibrium checking' else write (1,*) ' No equilibrium checking' endif if (istyle.eq.1) write (1,*) 'Thermalization run' if (istyle.eq.2) write (1,*) 'Constant strain run' if (istyle.eq.3) write (1,*) 'Stress profile run' write (1,11000) nstep,tprint,ilist,dt,nconf,teq,gamma,peq, & pmass,dstrain,intstrn write (20,10000) "# file of h-matrix" write (21,10000) '# file of compacted strain matrix' write (24,10000) '# file of moduli' close(20) close (21) close (24) endif return 10000 format(a) 11000 format(' number of steps',t25,'=',i6,/ & ' print interval (fs)',t25,'=',f10.4,/' list interval',t25,'=',i6,/ & ' initial timestep',t25,'=',f10.4,/' number of configs',t25,'=',i6,/ & ' initial temperature',t25,'=',f10.4,/ & ' damping constant',t25,'=',f10.4,/ & ' desired pressure',t25,'=',f10.4,/ & ' piston mass',t25,'=',f10.4,/' delta strain',t25,'=',3f10.4/ & ' strain interval',t25,'=',i6) end subroutine fflush(n) #ifndef cray call flush(n) #endif return end SUBROUTINE INDEXX(N2,ARRIN,INDX) implicit double precision (a-h,o-z) DIMENSION ARRIN(N2) dimension indx(N2) DO J=1,N2 INDX(J) = J enddo L = N2/2 + 1 IR = N2 100 CONTINUE IF (L.GT.1) THEN L = L - 1 INDXT = INDX(L) Q = ARRIN(INDXT) ELSE INDXT = INDX(IR) Q = ARRIN(INDXT) INDX(IR) = INDX(1) IR = IR - 1 IF (IR.EQ.1) THEN INDX(1) = INDXT RETURN ENDIF ENDIF I = L J = L + L 200 IF (J.LE.IR) THEN IF (J.LT.IR) THEN IF (ARRIN(INDX(J)).LT.ARRIN(INDX(J+1))) J = J + 1 ENDIF IF (Q.LT.ARRIN(INDX(J))) THEN INDX(I) = INDX(J) I = J J = J + J ELSE J = IR + 1 ENDIF GOTO 200 ENDIF INDX(I) = INDXT GOTO 100 END subroutine kinet(nmol,nttot,v,ia,prs,amass,rk,boxk) implicit double precision (a-h,o-z) c c routine to compute the kinetic energy c c prs is in Voight notation: c prs(1) = vx*vx c prs(2) = vy*vy c prs(3) = vz*vz c prs(4) = vy*vz c prs(5) = vz*vx c prs(6) = vx*vy c parameter (HARTREE=27.21078584d0) dimension v(3,nmol),prs(6),amass(nmol) dimension ia(nmol) integer iherm(6) data iherm/1,6,2,5,4,3/ rk = 0.0 do k=1,nmol rm = amass(k) ik = 0 do ii=1,3 do jj=1,ii ik = ik + 1 prs(iherm(ik)) = rm*v(jj,k)*v(ii,k) + prs(iherm(ik)) enddo rk = rk + 0.5d0*rm*v(ii,k)*v(ii,k) enddo enddo boxk = 0.d0 return end subroutine matinv(a,ainv,det) implicit double precision (a-h,o-z) c.... c.... routine to compute the inverse and determinant of a 3x3 matrix c.... ainv = (cofactor(a))transpose / det c... dimension a(3,3),ainv(3,3) d11 = a(2,2)*a(3,3) - a(2,3)*a(3,2) d22 = a(3,3)*a(1,1) - a(3,1)*a(1,3) d33 = a(1,1)*a(2,2) - a(1,2)*a(2,1) d12 = a(2,3)*a(3,1) - a(2,1)*a(3,3) d23 = a(3,1)*a(1,2) - a(3,2)*a(1,1) d31 = a(1,2)*a(2,3) - a(1,3)*a(2,2) d13 = a(2,1)*a(3,2) - a(3,1)*a(2,2) d21 = a(3,2)*a(1,3) - a(1,2)*a(3,3) d32 = a(1,3)*a(2,1) - a(2,3)*a(1,1) det = a(1,1)*d11 + a(1,2)*d12 + a(1,3)*d13 ainv(1,1) = d11/det ainv(2,2) = d22/det ainv(3,3) = d33/det ainv(1,2) = d21/det ainv(2,1) = d12/det ainv(1,3) = d31/det ainv(3,1) = d13/det ainv(2,3) = d32/det ainv(3,2) = d23/det det = abs(det) return end subroutine summations(natom,rnuc,rold,force,fvold,du,dusum) implicit double precision (a-h,o-z) parameter (EtoeV=6.241808d0) dimension rnuc(3,natom) dimension rold(3,natom) dimension force(3,natom) dimension fvold(3,natom) dusum = 0.d0 do kk=1,natom dusum = dusum - (rnuc(1,kk)-rold(1,kk))* & (fvold(1,kk)+force(1,kk))/2.d0 dusum = dusum - (rnuc(2,kk)-rold(2,kk))* & (fvold(2,kk)+force(2,kk))/2.d0 dusum = dusum - (rnuc(3,kk)-rold(3,kk))* & (fvold(3,kk)+force(3,kk))/2.d0 enddo dusum = dusum*EtoeV diff = du - dusum write (*,10000) istep,dusum,du,diff return 10000 format(i6,3(2x,1pe13.5)) end subroutine tcheck(SIMPLE,natom,en,teq,LENGTH,RTOL,ANSWER) implicit double precision (a-h,o-z) logical SIMPLE,ANSWER parameter (MAXLENGTH=10000) dimension vector(MAXLENGTH) data kount/0/,old/0.d0/ data sqrtth/1.22474487139158904909/ save kount,old,vector ANSWER = .FALSE. if (SIMPLE) then diff = abs(en-old) TOL = ABS(RTOL*old) old = en if (diff.lt.TOL) then kount = kount + 1 else kount = 0 endif if (kount.ge.LENGTH) then ANSWER = .TRUE. kount = 0 endif return else kount = kount + 1 if (kount.le.LENGTH) then vector(kount) = en return else im = mod(kount,LENGTH) + 1 vector(im) = en sumx = 0.d0 sumy = 0.d0 sumxx = 0.d0 sumxy = 0.d0 sumyy = 0.d0 do k=1,LENGTH sumy = sumy + vector(k) sumyy = sumyy + vector(k)**2 sumx = sumx + k sumxx = sumx + k*k sumxy = sumxy + k*vector(k) enddo slope = (LENGTH*sumxy-sumx*sumy)/(LENGTH*sumxx-sumx*sumx) varsq = (sumyy-sumy*sumy/LENGTH)/(LENGTH-1) TOL = 3*teq/(2*natom*LENGTH) TOL = max(1.0d-6, TOL) if (varsq.lt.TOL) then ANSWER = .TRUE. kount = 0 endif endif return endif end subroutine thermal(nmol,itemp,v,hv,fv,rk,teq,temp,gamma) implicit double precision (a-h,o-z) include "PARAMETERS.h" parameter (EtoeV=6.241808d0) dimension fv(3,nmol),hv(3,3) dimension v(3,nmol) parameter (EtoK=72435.42d0) if (temp.NE.0.0) then if (gamma.NE.0.0) then if (teq.EQ.0.0) then gstar = 1.0 - 0.5*gamma else gstar = sqrt(1.0+(teq/temp-1.0)*gamma) endif if (teq.EQ.0.0) then bstar = 1.0 - 0.5*gamma else bstar = sqrt(1.0+(teq/temp-1.0)*gamma) endif if (itemp.NE.3) then call ascale(3*nmol,gstar,v) call ascale(9,bstar,hv) endif endif if (itemp.NE.1) then do j=1,nmol if(adot(3,fv(1,j),v(1,j)).LE.0.0) then v(1,j) = 0.0 v(2,j) = 0.0 v(3,j) = 0.0 endif enddo endif endif return end subroutine verlet(nmol,coord,vel,force,ia,amass,dt,iflag, &rmove,qup) implicit double precision (a-h,o-z) include "PARAMETERS.h" dimension coord(3,nmol),vel(3,nmol),force(3,nmol),amass(*) dimension ia(nmol) dimension rmove(3,nmol) logical QUP,UPX,UPY,UPZ c c if iflag.eq.1 then we are starting the verlet step... c QUP=.FALSE. do k=1,nmol fff = dt/(2.d0*amass(k)) vel(1,k) = vel(1,k) + fff*force(1,k) vel(2,k) = vel(2,k) + fff*force(2,k) vel(3,k) = vel(3,k) + fff*force(3,k) enddo if (iflag.eq.1) then do k=1,nmol dx = dt*vel(1,k) dy = dt*vel(2,k) dz = dt*vel(3,k) coord(1,k) = coord(1,k) + dx coord(2,k) = coord(2,k) + dy coord(3,k) = coord(3,k) + dz rmove(1,k) = rmove(1,k) + dx rmove(2,k) = rmove(2,K) + dy rmove(3,k) = rmove(3,k) + dz enddo do k = 1, nmol UPX = (abs(rmove(1,k)).GT.TOLMOVE) UPY = (abs(rmove(2,k)).GT.TOLMOVE) UPZ = (abs(rmove(3,k)).GT.TOLMOVE) if(UPX.OR.UPY.OR.UPZ) QUP=.TRUE. enddo endif return end subroutine writcn(NDIM,ntype,atsym,r,v,h0,nmol,ia,di,time, & filename) implicit double precision (a-h,o-z) character*9 filename dimension r(3,nmol),v(3,nmol) dimension h0(3,3) dimension ia(nmol) dimension di(nmol) character*2 atsym(ntype) c.... c.... routine to write final output configuration to unit 99 c.... parameter (ncn=99) open (unit=ncn,file=filename,status='unknown') write (ncn,10000) time,ndim write (ncn,11000) h0 do i=1,nmol write (ncn,12000) atsym(ia(i)),(r(k,i),k=1,3),di(i),(v(j,i), & j=1,3) enddo close (ncn) return 10000 format('#',f12.3/i3) 11000 format(3f20.15) 12000 format(a2,7(2x,1pe21.14)) end subroutine writxyz(nmol,ntype,r,v,nnuc,ia,di,time,filename) implicit double precision (a-h,o-z) character*7 filename dimension r(3,nmol) dimension v(3,nmol) dimension di(nmol) dimension ia(nmol) dimension nnuc(ntype) c.... c.... routine to write final output configuration to unit 99 c.... parameter (ncn=8) open (unit=ncn,file=filename,status='unknown',access='append') write (ncn,*) nmol write (ncn,*) 'time=',time do i=1,nmol write (ncn,10000) nnuc(ia(i)),r(1,i),r(2,i),r(3,i),di(i), & v(1,i),v(2,i),v(3,i) enddo close (8) return 10000 format(i2,7(1x,1pe13.6)) end subroutine MkStrain(NDIM,natom,strain,ucell,kcell,Omega,rnuc, & nrtot,rrvec,nktot,rkvec) implicit double precision (a-h,o-z) double precision ucell(3,3),kcell(3,3) dimension rnuc(3,natom) dimension rrvec(0:3,nrtot) dimension rkvec(0:3,nktot) dimension strain(6) dimension a(3),RotMat(3,3) dimension eps(6) dimension trans(3,3) PI = 3.14159265358979323846 TWOPI = 2*PI if (NDIM.eq.2) then eps(1) = strain(1) eps(2) = strain(2) eps(3) = strain(3) eps(4) = 0.d0 eps(5) = 0.d0 eps(6) = strain(6) call MkTrans(3,eps,trans) else if (NDIM.eq.3) then call MkTrans(3,strain,trans) else print *,'MkStrain-invalid value of NDIM' stop endif call Rotate(3,trans,ucell) call Rotate(natom,trans,rnuc) call RecLatt(NDIM,Omega,ucell,kcell) call RecLatt(3,dummy,trans,RotMat) do n=1,nrtot do k=1,3 a(k) = 0.d0 enddo sum = 0.d0 do k=1,3 do j=1,3 a(k) = a(k) + trans(k,j)*rrvec(j,n) enddo sum = sum + a(k)**2 enddo do k=1,3 rrvec(k,n) = a(k) enddo rrvec(0,n) = SQRT(sum) enddo do n=1,nktot do k=1,3 a(k) = 0.d0 enddo sum = 0.d0 do k=1,3 do j=1,3 a(k) = a(k) + RotMat(j,k)*rkvec(j,n) enddo enddo do k=1,3 rkvec(k,n) = a(k)/TWOPI sum = sum + a(k)**2 enddo rkvec(0,n) = SQRT(sum)/TWOPI enddo return end subroutine MkTrans(n,strain,trans) implicit double precision (a-h,o-z) dimension trans(3,3) dimension strain(6) trans(1,1) = 1.d0 + strain(1) trans(2,2) = 1.d0 + strain(2) trans(3,3) = 1.d0 + strain(3) trans(2,3) = strain(4)/2.d0 trans(3,2) = strain(4)/2.d0 trans(1,3) = strain(5)/2.d0 trans(3,1) = strain(5)/2.d0 trans(1,2) = strain(6)/2.d0 trans(2,1) = strain(6)/2.d0 return end subroutine Rotate(n,RotMat,a) implicit double precision (a-h,o-z) dimension RotMat(3,3),a(3,n) dimension b(3) do i=1,n do k=1,3 b(k) = 0.d0 do j=1,3 b(k) = b(k) + RotMat(k,j)*a(j,i) enddo enddo do k=1,3 a(k,i) = b(k) enddo enddo return end subroutine OrAxis(a,b,RotMat) implicit double precision (a-h,o-z) parameter (TOLER=1.d-5) c rotate vectors so that vector a(*) now lies along vector b(*) c and return appropriate rotation matrix dimension RotMat(3,3) dimension a(3),b(3),c(3) call Nabla(c,a,b) sum1 = SQRT(a(1)*a(1)+a(2)*a(2)+a(3)*a(3)) sum2 = SQRT(b(1)*b(1)+b(2)*b(2)+b(3)*b(3)) sum3 = SQRT(c(1)*c(1)+c(2)*c(2)+c(3)*c(3)) dotp = a(1)*b(1) + a(2)*b(2) + a(3)*b(3) if (dotp.lt.0.d0) then dotp = -dotp sum3 = -sum3 endif Angle = DACOS(dotp/sum1/sum2) if (ABS(sum3).gt.TOLER) then do i=1,3 c(i) = c(i)/sum3 enddo endif if (ABS(Angle).gt.TOLER) then call RotAxis(c,Angle,.FALSE.,RotMat) else call MakUnit(RotMat) endif return end subroutine RotAxis(Axis,Angle,Reflct,RotMat) implicit double precision (a-h,o-z) dimension Axis(3),RotMat(3,3) dimension uvec(3) logical Reflct parameter (TOLER=1.d-5) if (ABS(Angle).lt.TOLER) then call MakUnit(RotMat) if (Reflct) then sum = SQRT(Axis(1)**2+Axis(2)**2+Axis(3)**2) do i=1,3 uvec(i) = Axis(i)/sum enddo do i=1,3 do j=1,3 RotMat(j,i) = RotMat(j,i) - 2*uvec(i)*uvec(j) enddo enddo endif return endif sum = SQRT(Axis(1)**2+Axis(2)**2+Axis(3)**2) do i=1,3 uvec(i) = Axis(i)/sum enddo COF = COS(Angle) do i=1,3 do j=1,3 RotMat(j,i) = uvec(i)*uvec(j)*(1.d0-COF) enddo RotMat(i,i) = RotMat(i,i) + COF enddo COF = SIN(Angle) RotMat(1,2) = RotMat(1,2) - uvec(3)*COF RotMat(2,3) = RotMat(2,3) - uvec(1)*COF RotMat(3,1) = RotMat(3,1) - uvec(2)*COF RotMat(1,3) = RotMat(1,3) + uvec(2)*COF RotMat(2,1) = RotMat(2,1) + uvec(3)*COF RotMat(3,2) = RotMat(3,2) + uvec(1)*COF if (Reflct) then do i=1,3 do j=1,3 RotMat(j,i) = RotMat(j,i) - 2*uvec(i)*uvec(j) enddo enddo endif return end subroutine AlnAxis(RotMat,natoms,coord,RotStd) implicit double precision (a-h,o-z) dimension RotMat(3,3),RotStd(3,3) dimension coord(3,natoms) call Rotate(natoms,RotMat,coord) call Rotate(3,RotMat,RotStd) return end subroutine MakUnit(RotMat) implicit double precision (a-h,o-z) dimension RotMat(3,3) do i=1,3 do j=1,i-1 RotMat(j,i) = 0.d0 RotMat(i,j) = 0.d0 enddo RotMat(i,i) = 1.d0 enddo return end subroutine MakCart(n,vec) implicit double precision (a-h,o-z) dimension vec(3) do i=1,3 vec(i) = 0.d0 enddo vec(n) = 1.d0 return end subroutine Orient(natom,ucell,rnuc) implicit double precision (a-h,o-z) dimension ucell(3,3),rnuc(3,natom) dimension a(3),b(3),RotMat(3,3) call Nabla(b,ucell(1,1),ucell(1,2)) call VNorm(b) call MakCart(3,a) call OrAxis(b,a,RotMat) call Rotate(3,RotMat,ucell) call Rotate(natom,RotMat,rnuc) return end subroutine GetInpV(NDIM,MAXATM,MAXTYP,natom,ucell,kcell,Omega, & rnuc,zchg,vnuc,amass,nttot,ntype,nnuc,hard,chi,zeff,zeta, & potparm,time) implicit double precision (a-h,o-z) c c use rationalized atomic scale units c c time: femtoseconds [10**(-15) sec] c c distance: Angstrom [10**(-10) m] c c mass: 10**(-28) kg c c energy: 10**(-18) J [62.41808 eV] c c velocity A/fs [10**5 m/sec] c c pressure: 10**(12) Pa (1 TPa, 10 Mbar) c parameter (EtoK=72435.42d0) parameter (EtoeV=6.241808d0) parameter (AMUtoM=16.6053d0) parameter (EtoH=0.229387274d0) parameter (AUtoA=0.52917d0) double precision ucell(3,3),kcell(3,3),Omega dimension rnuc(3,MAXATM) dimension zchg(MAXATM) dimension vnuc(3,MAXATM) dimension amass(MAXATM) dimension ntype(MAXATM) dimension nnuc(MAXTYP) dimension hard(MAXTYP) dimension chi(MAXTYP) dimension zeff(MAXTYP) dimension zeta(MAXTYP) dimension potparm(10,3) dimension snuc(3) dimension pnuc(3) character*240 line c if 1st line begins with #, use RdInpt() call getline(5,nline,line) if (line(1:1).eq.'#') then if(nline.gt.1) then read(line(2:nline),*) time endif call RdInpV(NDIM,MAXATM,MAXTYP,natom,ucell,kcell,Omega,rnuc, & zchg,vnuc,amass,nttot,ntype,nnuc,hard,chi,zeff,zeta,potparm) return endif c read in number of dimensions read (line,*) NDIM c read in periodic boundary conditions read (5,*) ucell c read in atomic positions natom = 0 nttot = 0 do i=1,MAXATM read (5,*) nn,snuc, pnuc if (nn.le.0) goto 200 natom = natom + 1 do k=1,3 rnuc(k,i) = 0.d0 do j=1,3 rnuc(k,i) = rnuc(k,i) + snuc(j)*ucell(k,j) vnuc(k,i) = vnuc(k,i) + pnuc(j)*ucell(k,j) enddo enddo do n=1,nttot if (nn.eq.nnuc(n)) then ntype(i) = n goto 100 endif enddo nttot = nttot + 1 if (nttot.gt.MAXTYP) then print *,'Too many atom types used: ',MAXTYP stop endif nnuc(nttot) = nn ntype(i) = nttot 100 continue enddo 200 continue c read in values of hard(*),chi(*),zeff(*),zeta(*) do n=1,nttot nnuc(n) = -nnuc(n) enddo do n=1,nttot read (5,*,END=300,ERR=300) nn,cc,hh,xx,zz do i=1,nttot if (nn.eq.-nnuc(i)) then nnuc(i) = nn hard(i) = hh/EtoeV chi(i) = cc/EtoeV zeff(i) = zz zeta(i) = 1.5d0/xx goto 400 endif enddo 300 print *,'Error in reading Coulomb parameters' stop 400 continue enddo c set atomic masses do i=1,natom amass(i) = 0.d0 enddo c read in values of potential parameters c c potparm(1,*): A c potparm(2,*): r* c potparm(3,*): alpha c potparm(4,*): beta c potparm(5,*): B c potparm(6,*): eta c potparm(7,*): C c nptot = nttot*(nttot+1)/2 do n=1,nptot read (5,*) Ea,rstar,al,beta,Eb,eta,Ec potparm(1,n) = Ea/EtoeV potparm(2,n) = rstar potparm(3,n) = al potparm(4,n) = beta potparm(5,n) = Eb/EtoeV potparm(6,n) = eta potparm(7,n) = Ec/EtoeV enddo c set up reciprocal lattice vectors call RecLatt(NDIM,Omega,ucell,kcell) c get atomic masses do n=1,natom amass(n) = GetMass(nnuc(ntype(n))) enddo return end subroutine RdInpV(NDIM,MAXATM,MAXTYP,natom,ucell,kcell,Omega,rnuc, & zchg,vnuc,amass,nttot,ntype,nnuc,hard,chi,zeff,zeta,potparm) implicit double precision (a-h,o-z) c c use rationalized atomic scale units c c time: femtoseconds [10**(-15) sec] c c distance: Angstrom [10**(-10) m] c c mass: 10**(-28) kg c c energy: 10**(-18) J [62.41808 eV] c c velocity A/fs [10**5 m/sec] c c pressure: 10**(12) Pa (1 TPa, 10 Mbar) c parameter (EtoK=72435.42d0) parameter (EtoeV=6.241808d0) parameter (AMUtoM=16.6053d0) parameter (EtoH=0.229387274d0) parameter (AUtoA=0.52917d0) double precision ucell(3,3),kcell(3,3),Omega dimension rnuc(3,MAXATM) dimension vnuc(3,MAXATM) dimension zchg(MAXATM) dimension amass(MAXATM) dimension ntype(MAXATM) dimension nnuc(MAXTYP) dimension hard(MAXTYP) dimension chi(MAXTYP) dimension zeff(MAXTYP) dimension zeta(MAXTYP) dimension potparm(10,3) integer GtAtNo character*2 symbol character*240 line logical isdigit parameter (HARTREE=27.21078584d0) parameter (AU=0.52917d0) c read in number of dimensions read (5,*) NDIM c read in periodic boundary conditions read (5,*) ucell c read in atomic positions natom = 0 nttot = 0 call getline(5,nline,line) nn = 0 do while ((nline.gt.0).and.(natom.lt.MAXATM)) natom = natom + 1 if (isdigit(line(1:1)).or.isdigit(line(2:2))) then if(isspace(line(1:1))) nn=1 call GtIntg(nn,nline,line,iatno) call GtAtSy(symbol,iatno) else call GtLbl(nn,nline,line,symbol) endif do k=1,3 call GtReal(nn,nline,line,rnuc(k,natom)) enddo if(nn.lt.nline) then call GtReal(nn,nline,line,zchg(natom)) do k=1,3 call GtReal(nn,nline,line,vnuc(k,natom)) enddo else zchg(natom)=0.0 vnuc(1,natom)=0.0 vnuc(2,natom)=0.0 vnuc(3,natom)=0.0 endif nz = GtAtNo(symbol) do n=1,nttot if (nz.eq.nnuc(n)) then ntype(natom) = n goto 100 endif enddo nttot = nttot + 1 if (nttot.gt.MAXTYP) then print *,'Too many atom types used: ',MAXTYP stop endif nnuc(nttot) = nz ntype(natom) = nttot 100 continue call getline(5,nline,line) nn = 0 if ((natom.ge.MAXATM).and.(nline.gt.0)) then print *,'Number of atoms limit exceeded' stop endif enddo c read in potential parameters call RdPot(nttot,nnuc,chi,hard,zeta,zeff,potparm) c set up reciprocal lattice vectors call RecLatt(NDIM,Omega,ucell,kcell) c get atomic masses do n=1,natom amass(n) = GetMass(nnuc(ntype(n))) enddo return end subroutine escharge(NDIM,natom,nbrtot,nbrlist,xnbr,nrtot,rrvec, & nktot,rkvec,Omega,rnuc,zchg,ntype,hard,chi,zeff,zeta,alphew, & ees,eew,vewald,catom,vwork,work) implicit double precision (a-h,o-z) dimension nbrlist(0:2,nbrtot) dimension xnbr(0:3,nbrtot) dimension rrvec(0:3,nrtot) dimension rkvec(0:3,nktot) dimension zchg(natom) dimension rnuc(3,natom) dimension ntype(natom) dimension hard(*),chi(*),zeff(*),zeta(*) dimension vewald(*),catom(natom) dimension vwork(*) dimension work(natom,3) double precision xlambda dimension cpu(0:4) parameter (EtoK=72435.42d0) parameter (EtoeV=6.241808d0) parameter (AMUtoM=16.6053d0) parameter (EtoH=0.229387274d0) parameter (AUtoA=0.52917d0) c call timex(t1) do i=1,natom catom(i) = 0.d0 enddo nherm = natom*(natom+1)/2 call ewald(NDIM,natom,nherm,rnuc,nbrtot,nbrlist,xnbr,nrtot,rrvec, & nktot,rkvec,Omega,alphew,vewald) c call timex(t2) call elocal(natom,rnuc,nbrtot,nbrlist,xnbr,ntype,zeta,zeff,e0, & catom,vewald) c call timex(t3) e0 = e0*AutoA/EtoH call aScale(natom,AutoA/EtoH,catom) call aScale(nherm,AutoA/EtoH,vewald) ij = 0 do i=1,natom ni = ntype(i) ij = ij + i catom(i) = catom(i) + chi(ni) vewald(ij) = vewald(ij) + hard(ni) enddo do i=1,natom work(i,1) = -catom(i) work(i,2) = 1 enddo do i=1,nherm vwork(i) = vewald(i) enddo c call timex(t1p) call DSPSV('U',natom,2,vwork,work(1,3),work,natom,info) c call timex(t2p) sum1 = 0.d0 sum2 = 0.d0 do i=1,natom sum1 = sum1 + work(i,1) sum2 = sum2 + work(i,2) enddo xlambda = -sum1/sum2 do i=1,natom zchg(i) = work(i,1) + xlambda*work(i,2) enddo ees = 0.d0 do i=1,natom ees = ees + catom(i)*zchg(i) enddo ees = ees/2.d0 c call timex(t4) do i=1,nherm vwork(i) = vewald(i) enddo ij = 0 do i=1,natom ni = ntype(i) ij = ij + i vwork(ij) = vewald(ij) - hard(ni) enddo eew = ewsum(natom,zchg,vwork) c call timex(t5) cpu(0) = t5 - t1 cpu(1) = t2 - t1 cpu(2) = t3 - t2 cpu(3) = t2p-t1p cpu(4) = t5 - t4 c print 111, cpu 111 format('total: ',f8.3,' ewald:',f8.3,' elocal:',f8.3, & ' DSPSV:',f8.3,' esum:',f8.3) return end subroutine elocal(natom,rnuc,nbrtot,nbrlist,xnbr,ntype,zeta,zeff, & e0,catom,vewald) implicit double precision (a-h,o-z) dimension ntype(*),zeta(*),zeff(*) dimension nbrlist(0:2,nbrtot) dimension xnbr(0:3,nbrtot) dimension rnuc(3,natom) dimension vewald(*) dimension catom(natom) e0 = 0.d0 do m=1,nbrtot i = nbrlist(1,m) j = nbrlist(2,m) ni = ntype(i) nj = ntype(j) dist = xnbr(0,m) call getvij(dist,zeta(ni),zeta(nj),vi,vj,vij) if (i.eq.j) then ij = i*(i+1)/2 vewald(ij) = vewald(ij) + 2*vij catom(i) = catom(i) + 2*zeff(nj)*(vi-vij) else m1 = MIN(i,j) m2 = MAX(i,j) ij = m2*(m2-1)/2 + m1 vewald(ij) = vewald(ij) + vij catom(i) = catom(i) + zeff(nj)*(vi-vij) catom(j) = catom(j) + zeff(ni)*(vj-vij) endif e0 = e0 + 2*zeff(ni)*zeff(nj)*(vij-vi-vj) enddo return end subroutine hermprint(array,num) double precision array dimension array(*) k = 0 is= 1 do i = 1,num k = k + i print 111, (array(j),j=is,k) is = k + 1 enddo 111 format(f11.5,2x) return end subroutine pprint(array,num) double precision array dimension array(*) print 111, (array(j),j=1,num) 111 format(f11.5,2x) return end