diff --git a/source/bath.i b/source/bath.i old mode 100644 new mode 100755 index f06f93ef5..f9fddb7be --- a/source/bath.i +++ b/source/bath.i @@ -35,24 +35,38 @@ c anisotrop logical flag governing use of anisotropic pressure c thermostat choice of temperature control method to be used c barostat choice of pressure control method to be used c volscale choice of scaling method for Monte Carlo barostat +c q_iso1/2 masses of isokinetic thermostats +c v_iso1/2 velocities of isokinetic thermostats +c f_iso1/2 forces acting on isokinetic thermostats +c len_nhc In isokinetic dynamics, this is the Nose-Hoover chain length M. +c In stochastic isokinetic dynamics, this is the chain length L c -c - integer maxnose - parameter (maxnose=4) - integer voltrial + real*8 q_iso1,q_iso2,v_iso1,v_iso2 + real*8 vnh,qnh,gnh real*8 kelvin,atmsph real*8 tautemp,taupres real*8 compress,collide - real*8 vnh,qnh,gnh real*8 vbar,qbar,gbar real*8 eta,volmove + integer maxnose + parameter (maxnose=4) + integer len_nhc,isok_L,isok_M + integer isok_chain + parameter (isok_chain=4) + integer voltrial logical isothermal logical isobaric logical anisotrop character*9 volscale character*11 barostat character*11 thermostat - common /bath/ kelvin,atmsph,tautemp,taupres,compress,collide, + common /bath/ q_iso1(isok_chain,3,maxatm), + & q_iso2(isok_chain,3,maxatm), + & v_iso1(isok_chain,3,maxatm), + & v_iso2(isok_chain,3,maxatm), + & f_iso1(isok_chain,3,maxatm), + & kelvin,atmsph,tautemp,taupres,compress,collide, & vnh(maxnose),qnh(maxnose),gnh(maxnose),vbar,qbar, - & gbar,eta,volmove,voltrial,isothermal,isobaric, - & anisotrop,thermostat,barostat,volscale + & gbar,eta,volmove,voltrial,len_nhc,isothermal, + & isobaric,anisotrop,thermostat, + & barostat,volscale diff --git a/source/dynamic.f b/source/dynamic.f old mode 100644 new mode 100755 index cc5945241..cd3527a3f --- a/source/dynamic.f +++ b/source/dynamic.f @@ -225,6 +225,7 @@ program dynamic end do end if end if + c c initialize any holonomic constraints and setup dynamics c @@ -261,9 +262,17 @@ program dynamic write (iout,390) 390 format (/,' Molecular Dynamics Trajectory via', & ' r-RESPA MTS Algorithm') - else + else if (integrate .eq. 'STOCH-ISOK') then write (iout,400) 400 format (/,' Molecular Dynamics Trajectory via', + & ' Stochastic Isokinetic Algorithm') + else if (integrate .eq. 'ISOKINETIC') then + write (iout,410) + 410 format (/,' Molecular Dynamics Trajectory via', + & ' Isokinetic NHC Algorithm') + else + write (iout,420) + 420 format (/,' Molecular Dynamics Trajectory via', & ' Modified Beeman Algorithm') end if c @@ -284,6 +293,10 @@ program dynamic call rgdstep (istep,dt) else if (integrate .eq. 'RESPA') then call respa (istep,dt) + else if (integrate .eq. 'ISOKINETIC') then + call isokinetic (istep,dt) + else if (integrate .eq. 'STOCH-ISOK') then + call stochisok (istep,dt) else call beeman (istep,dt) end if diff --git a/source/isokinetic.f b/source/isokinetic.f new file mode 100644 index 000000000..b5c5e978a --- /dev/null +++ b/source/isokinetic.f @@ -0,0 +1,480 @@ +c +c +c ################################################### +c ## COPYRIGHT (C) 2011 by Jay William Ponder ## +c ## All Rights Reserved ## +c ################################################### +c +c ############################################################# +c ## ## +c ## subroutine isokinetic ## +c ## ## +c ############################################################# +c +c +c + subroutine isokinetic (istep,dt) + implicit none + include 'sizes.i' + include 'atmtyp.i' + include 'freeze.i' + include 'moldyn.i' + include 'units.i' + include 'virial.i' + include 'atoms.i' + include 'bath.i' + include 'bond.i' + include 'bound.i' + include 'inform.i' + include 'iounit.i' + include 'keys.i' + include 'mdstuf.i' + include 'potent.i' + include 'solute.i' + include 'stodyn.i' + include 'usage.i' + + integer i,j,k + integer istep + integer nalt + real*8 dt,dt_2,tempstor,temp2 + real*8 dta,dta_2,kBT,LkT + real*8 epot,etot + real*8 eksum,eps + real*8 temp,pres + real*8 ealt,dalt + real*8 ekin(3,3) + real*8 stress(3,3) + real*8 viralt(3,3) + real*8, allocatable :: xold(:) + real*8, allocatable :: yold(:) + real*8, allocatable :: zold(:) + real*8, allocatable :: derivs(:,:) +10 format (3f16.14) + +c +c set some values for the dynamics integration +c + dt_2 = 0.5d0 * dt + dta = dt / (dble(nrespa)) + dta_2 = 0.5d0 * dta + + kBT = boltzmann*kelvin + LkT = len_nhc*kBT +c +c perform dynamic allocation of some local arrays +c + allocate (xold(n)) + allocate (yold(n)) + allocate (zold(n)) + allocate (derivs(3,n)) + +c +c initialize virial from fast-evolving potential energy terms +c + do i = 1, 3 + do j = 1, 3 + viralt(j,i) = 0.0d0 + end do + end do + +c +c +c begin inner RESPA loop +c + do k = 1, nrespa + do i = 1, n + if (use(i)) then + xold(i) = x(i) + yold(i) = y(i) + zold(i) = z(i) + end if + end do + +c +c +c nose-like isokinetic operator +c + call appisok (dta) +c +c +c force-like isokinetic operator +c + + call appisok2 (dta) +c +c Evolve positions +c + if (use_rattle) then + do i = 1, n + xold(i) = x(i) + yold(i) = y(i) + zold(i) = z(i) + x(i) = x(i) + v(1,i)*dta + y(i) = y(i) + v(2,i)*dta + z(i) = z(i) + v(3,i)*dta + end do + call rattle (dta,xold,yold,zold) + else + do i = 1, n + if (use(i)) then + x(i) = x(i) + v(1,i)*dta + y(i) = y(i) + v(2,i)*dta + z(i) = z(i) + v(3,i)*dta + end if + end do + end if +c +c get the potential energy and atomic forces +c + if (k.eq.nrespa) then + call gradfast (ealt,derivs) + + if (use_rattle) then + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + call gradslow (epot,derivs) + do i = 1, n + do j = 1, 3 + a(j,i) = a(j,i) - ( nrespa * convert + & * derivs(j,i) / mass(i) ) + end do + end do + call rattle2 (dta) + else + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + call gradslow (epot,derivs) + do i = 1, n + do j = 1, 3 + a(j,i) = a(j,i) - ( nrespa * convert + & * derivs(j,i) / mass(i) ) + end do + end do + end if + + else + call gradfast (ealt,derivs) + if (use_rattle) then + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + call rattle2 (dta) + else + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + end if + end if + +c +c +c force-like isokinetic operator +c + call appisok2 (dta) +c +c +c nose-like isokinetic operator +c + call appisok (dta) + + end do + +c +c find the constraint-corrected full-step velocities +c + if (use_rattle) call rattle2 (dt) +c +c total potential and virial from sum of fast and slow parts +c + epot = epot + ealt + do i = 1, 3 + do j = 1, 3 + vir(j,i) = vir(j,i) + viralt(j,i) + end do + end do + + +c +c perform deallocation of some local arrays +c + + deallocate (derivs) + deallocate (xold) + deallocate (yold) + deallocate (zold) +c +c total energy is sum of kinetic and potential energies +c + etot = eksum + epot + + call kinetic (eksum,ekin) + call pressure (dt,epot,ekin,temp,pres,stress) + temp = 2.0d0 * eksum / (dble(nfree) * gasconst) + + +c tempstor=0.0d0 +c do i = 1, n +c tempstor=tempstor+(mass(i)*v(1,i)*v(1,i)) +c temp2=0.5*q_iso1(1,1,i)*v_iso1(1,1,i)*v_iso1(1,1,i) +c tempstor=tempstor+temp2 +c +c tempstor=tempstor+(mass(i)*v(2,i)*v(2,i)) +c temp2=0.5*q_iso1(1,2,i)*v_iso1(1,2,i)*v_iso1(1,2,i) +c tempstor=tempstor+temp2 +c +c tempstor=tempstor+(mass(i)*v(3,i)*v(3,i)) +c temp2=0.5*q_iso1(1,3,i)*v_iso1(1,3,i)*v_iso1(1,3,i) +c tempstor=tempstor+temp2 +c end do +c write (iout,10) tempstor/(3.0*n*kBT) + + +c +c compute statistics and save trajectory for this step +c + call mdstat (istep,dt,etot,epot,eksum,temp,pres) + call mdsave (istep,dt,epot) + call mdrest (istep) + + + + + return + end + + +c +c ################################################################## +c ## ## +c ## subroutine applyisok ## +c ## ## +c ################################################################## +c +c +c "applyisok" represents the Liuville operator iL_N in Leimkuhler, Margul, and Tuckerman (2013) +c +c + subroutine appisok (dt) + implicit none + include 'sizes.i' + include 'atmtyp.i' + include 'freeze.i' + include 'moldyn.i' + include 'units.i' + include 'virial.i' + include 'atoms.i' + include 'bath.i' + include 'bond.i' + include 'bound.i' + include 'inform.i' + include 'iounit.i' + include 'keys.i' + include 'mdstuf.i' + include 'potent.i' + include 'solute.i' + include 'stodyn.i' + include 'usage.i' + integer nc,ns,i,j,k,iatom,inhc,ixyz,index + integer len_m1,len_p1,len_m2 + real*8 dt,dtc,dts + real*8 dt2,dt4,dt8 + real*8 tempstor,kBT,LkT,temp2,temp3 + real*8 w(3) + + len_m1=len_nhc-1 + len_p1=len_nhc+1 + len_m2=len_nhc-2 + + kBT = boltzmann*kelvin + LkT = len_nhc*kBT + + nc = 5 + ns = 3 + dtc = dt / dble(nc) + w(1) = 1.0d0 / (2.0d0-2.0d0**(1.0d0/3.0d0)) + w(2) = 1.0d0 - 2.0d0*w(1) + w(3) = w(1) + + +c +c use multiple time steps and Suzuki-Yoshida decomposition +c + do k = 1, nc + do j = 1, ns + + dts = w(j) * dtc + dt2 = 0.5d0 * dts + dt4 = 0.25d0 * dts + dt8 = 0.125d0 * dts + + do iatom = 1, n + do ixyz = 1, 3 + do inhc = 2, len_nhc + f_iso1(inhc,ixyz,iatom)=((q_iso1(inhc,ixyz,iatom) + & *v_iso1(inhc-1,ixyz,iatom) + & *v_iso1(inhc-1,ixyz,iatom))-kBT) + end do + end do + end do + + do iatom = 1, n + do ixyz = 1, 3 + v_iso1(len_nhc,ixyz,iatom)=v_iso1(len_nhc,ixyz,iatom) + & +dt4*f_iso1(len_nhc,ixyz,iatom) + & /q_iso1(len_nhc,ixyz,iatom) + end do + end do + + + do iatom = 1, n + do inhc = len_m1, 2, -1 + do ixyz = 1, 3 + index=inhc+1 + tempstor = exp(-dt8*v_iso1(index,ixyz,iatom)) + v_iso1(inhc,ixyz,iatom)=v_iso1(inhc,ixyz,iatom) + & *tempstor*tempstor + v_iso1(inhc,ixyz,iatom)=v_iso1(inhc,ixyz,iatom) + & +dt4*tempstor*f_iso1(inhc,ixyz,iatom) + & /q_iso1(index,ixyz,iatom) + end do + end do + end do + + + do iatom = 1, n + do ixyz = 1, 3 + temp2=v(ixyz,iatom)*v(ixyz,iatom)*mass(iatom) + & +(0.5*q_iso1(1,ixyz,iatom)*v_iso1(1,ixyz,iatom) + & *v_iso1(1,ixyz,iatom) + & *exp(-2.0*dt2*v_iso1(2,ixyz,iatom))) + temp3=sqrt(kBT/temp2) + v(ixyz,iatom)=v(ixyz,iatom)*temp3 + v_iso1(1,ixyz,iatom)=v_iso1(1,ixyz,iatom) + & *temp3*exp(-dt2*v_iso1(2,ixyz,iatom)) + end do + end do + + do iatom = 1, n + do ixyz = 1, 3 + do inhc = 2, len_nhc + f_iso1(inhc,ixyz,iatom)=((q_iso1(inhc,ixyz,iatom) + & *v_iso1(inhc-1,ixyz,iatom) + & *v_iso1(inhc-1,ixyz,iatom))-kBT) + end do + end do + end do + + do iatom = 1, n + do inhc = 2, len_m1 + do ixyz = 1, 3 + index=inhc+1 + tempstor = exp(-dt8*v_iso1(index,ixyz,iatom)) + v_iso1(inhc,ixyz,iatom)=v_iso1(inhc,ixyz,iatom) + & *tempstor*tempstor + v_iso1(inhc,ixyz,iatom)=v_iso1(inhc,ixyz,iatom) + & +dt4*tempstor*f_iso1(inhc,ixyz,iatom) + & /q_iso1(index,ixyz,iatom) + end do + end do + end do + + do iatom = 1, n + do ixyz = 1, 3 + v_iso1(len_nhc,ixyz,iatom)=v_iso1(len_nhc,ixyz,iatom) + & +dt4*f_iso1(len_nhc,ixyz,iatom) + & /q_iso1(len_nhc,ixyz,iatom) + end do + end do + + end do + end do + + + return + end +c +c ################################################################## +c ## ## +c ## subroutine applyisok2 ## +c ## ## +c ################################################################## +c +c +c + subroutine appisok2 (dt) + implicit none + include 'sizes.i' + include 'atmtyp.i' + include 'freeze.i' + include 'moldyn.i' + include 'units.i' + include 'virial.i' + include 'atoms.i' + include 'bath.i' + include 'bond.i' + include 'bound.i' + include 'inform.i' + include 'iounit.i' + include 'keys.i' + include 'mdstuf.i' + include 'potent.i' + include 'solute.i' + include 'stodyn.i' + include 'usage.i' + real*8 kBT,LkT,isoA,isoB + real*8 dt,dt_2,rootB,arg + real*8 iso_s,iso_sdot + real*8 epot,etot + integer iatom,inhc,ixyz,i,j,k + + dt_2 = 0.5d0*dt + kBT = boltzmann*kelvin + LkT = len_nhc*kBT + + + do iatom = 1, n + do ixyz = 1, 3 + isoA = mass(iatom)*a(ixyz,iatom)*v(ixyz,iatom)/kBT + isoB = mass(iatom)*a(ixyz,iatom)*a(ixyz,iatom)/kBT + rootB = sqrt(isoB) + arg = dt_2*rootB + + if (arg > 0.00001) then + iso_s = (1.0/rootB)*sinh(arg) + + & (isoA/isoB)*(cosh(arg)-1.0) + iso_sdot = cosh(arg) + (isoA/rootB)*sinh(arg) + else + iso_s = ((((isoB*isoA/24.0)*dt_2 +isoB/6.0)*dt_2 + + & 0.5*isoA)*dt_2 + 1.0)*dt_2 + iso_sdot = (((isoB*isoA/6.0)*dt_2 +0.5*isoB)*dt_2 + + & isoA)*dt_2 + 1.0 + end if + + v(ixyz,iatom) = (v(ixyz,iatom) + ((a(ixyz,iatom))*iso_s)) + & /iso_sdot + v_iso1(1,ixyz,iatom)=v_iso1(1,ixyz,iatom)/iso_sdot + end do + end do + + + return + end diff --git a/source/mdinit.f b/source/mdinit.f old mode 100644 new mode 100755 index ff73528a2..a2b9eac7a --- a/source/mdinit.f +++ b/source/mdinit.f @@ -39,13 +39,14 @@ subroutine mdinit include 'units.i' include 'uprior.i' include 'usage.i' - integer i,j,k,idyn - integer size,next + integer i,j,k,kk,idyn + integer size,next,len_isok integer lext,freeunit real*8 e,ekt,qterm real*8 maxwell,speed real*8 vec(3) real*8, allocatable :: derivs(:,:) + real*8 tempstor,LkT,kT,random logical exist character*7 ext character*20 keyword @@ -68,6 +69,8 @@ subroutine mdinit use_pred = .false. polpred = 'LSQR' iprint = 100 + nrespa = 1 + len_nhc = 4 c c set default values for temperature and pressure control c @@ -90,6 +93,7 @@ subroutine mdinit voltrial = 20 volmove = 100.0d0 volscale = 'ATOMIC' + c c check for keywords containing any altered parameters c @@ -146,10 +150,16 @@ subroutine mdinit call getword (record,volscale,next) call upcase (volscale) else if (keyword(1:9) .eq. 'PRINTOUT ') then - read (string,*,err=10,end=10) iprint + read (string,*,err=10,end=10) iprint + else if (keyword(1:11) .eq. 'NHC-LENGTH ') then + read (string,*,err=10,end=10) len_nhc + else if (keyword(1:6) .eq. 'NRESPA ') then + read (string,*,err=10,end=10) nrespa end if 10 continue end do + + c c make sure all atoms or groups have a nonzero mass c @@ -167,7 +177,7 @@ subroutine mdinit end do else do i = 1, n - if (use(i) .and. mass(i).le.0.0d0 .and. atomic(i).ne.0) then + if (use(i) .and. mass(i).le.0.0d0) then mass(i) = 1.0d0 totmass = totmass + 1.0d0 write (iout,30) i @@ -361,7 +371,7 @@ subroutine mdinit allocate (derivs(3,n)) call gradslow (e,derivs) do i = 1, n - if (use(i) .and. mass(i).ne.0.0d0) then + if (use(i)) then speed = maxwell (mass(i),kelvin) call ranvec (vec) do j = 1, 3 @@ -372,23 +382,126 @@ subroutine mdinit do j = 1, 3 v(j,i) = 0.0d0 a(j,i) = 0.0d0 + aalt(j,i) = 0.0d0 end do end if end do call gradfast (e,derivs) do i = 1, n - if (use(i) .and. mass(i).ne.0.0d0) then + if (use(i)) then do j = 1, 3 aalt(j,i) = -convert * derivs(j,i) / mass(i) end do - else + end if + end do + deallocate (derivs) + if (nuse .eq. n) call mdrest (0) + +c +c set physical and extended velocities, +c physical accelerations, and extended +c masses for isokinetic dynamics +c + else if (integrate .eq. 'ISOKINETIC') then + kT = boltzmann*kelvin + len_isok=1 + LkT = len_isok*kT + allocate (derivs(3,n)) + call gradient (e,derivs) + do i = 1, n do j = 1, 3 - aalt(j,i) = 0.0d0 + a(j,i) = -convert * derivs(j,i) / mass(i) + do k = 1, len_nhc + q_iso1(k,j,i) = kT*tautemp*tautemp + q_iso2(k,j,i) = kT*tautemp*tautemp + speed = maxwell (q_iso1(k,j,i),kelvin) + call ranvec (vec) + v_iso1(k,j,i) = speed * vec(1) + v_iso2(k,j,i) = speed * vec(2) + end do + + v(j,i) = random () + do k = 1, len_isok + v_iso1(k,j,i) = random () + end do + tempstor = v(j,i) * v(j,i) + do k = 1, len_isok + tempstor = tempstor + (v_iso1(k,j,i)*v_iso1(k,j,i)) + end do + tempstor = sqrt(tempstor) + v(j,i) = v(j,i) / tempstor + do k = 1, len_isok + v_iso1(k,j,i) = v_iso1(k,j,i) / tempstor + end do + v(j,i) = v(j,i) / (sqrt(mass(i)/LkT)) + do k = 1, len_isok + tempstor=(sqrt(q_iso1(k,j,i) + & /(kT*dble(len_isok+1)))) + v_iso1(k,j,i) = v_iso1(k,j,i)/tempstor + end do + tempstor = mass(i)*v(j,i)*v(j,i) + do k = 1, len_isok + tempstor = tempstor + dble(len_isok)*q_iso1(k,j,i)* + & v_iso1(k,j,i)*v_iso1(k,j,i)/dble(len_isok+1) + end do + do k = 1, len_nhc + end do + + end do + end do + deallocate (derivs) + call prtdyn + if (nuse .eq. n) call mdrest (0) + +c +c set physical and extended velocities, +c physical accelerations, and extended +c masses for stochastic isokinetic dynamics +c + else if (integrate .eq. 'STOCH-ISOK') then + kT = boltzmann*kelvin + LkT = len_nhc*kT + allocate (derivs(3,n)) + call gradient (e,derivs) + do i = 1, n + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + do k = 1, len_nhc + q_iso1(k,j,i) = kT*tautemp*tautemp + q_iso2(k,j,i) = kT*tautemp*tautemp + v_iso2(k,j,i) = kT*random ()/q_iso2(k,j,i) + end do + v(j,i) = random () + do k = 1, len_nhc + v_iso1(k,j,i) = random () + end do + tempstor = v(j,i) * v(j,i) + do k = 1, len_nhc + tempstor = tempstor + (v_iso1(k,j,i)*v_iso1(k,j,i)) + end do + tempstor = sqrt(tempstor) + v(j,i) = v(j,i) / tempstor + do k = 1, len_nhc + v_iso1(k,j,i) = v_iso1(k,j,i) / tempstor + end do + v(j,i) = v(j,i) / (sqrt(mass(i)/LkT)) + do k = 1, len_nhc + tempstor=(sqrt(q_iso1(k,j,i)/(kT*dble(len_nhc+1)))) + v_iso1(k,j,i) = v_iso1(k,j,i)/tempstor + end do + tempstor = mass(i)*v(j,i)*v(j,i) + do k = 1, len_nhc + tempstor = tempstor + dble(len_nhc)*q_iso1(k,j,i)* + & v_iso1(k,j,i)*v_iso1(k,j,i)/dble(len_nhc+1) + end do end do - end if end do deallocate (derivs) + call prtdyn if (nuse .eq. n) call mdrest (0) + + + c c set velocities and accelerations for Cartesian dynamics c @@ -396,7 +509,7 @@ subroutine mdinit allocate (derivs(3,n)) call gradient (e,derivs) do i = 1, n - if (use(i) .and. mass(i).ne.0.0d0) then + if (use(i)) then speed = maxwell (mass(i),kelvin) call ranvec (vec) do j = 1, 3 @@ -415,6 +528,11 @@ subroutine mdinit deallocate (derivs) if (nuse .eq. n) call mdrest (0) end if + + + + + c c check for any prior dynamics coordinate sets c diff --git a/source/mdstuf.i b/source/mdstuf.i index 2db0c60f1..6b21932ef 100644 --- a/source/mdstuf.i +++ b/source/mdstuf.i @@ -20,14 +20,14 @@ c velsave logical flag to save velocity vector components c frcsave logical flag to save force vector components c uindsave logical flag to save induced atomic dipoles c integrate type of molecular dynamics integration algorithm +c nrespa number of small time steps per large time step in RESPA c -c - integer nfree,irest + integer nfree,irest,nrespa integer bmnmix logical dorest logical velsave logical frcsave logical uindsave character*11 integrate - common /mdstuf/ nfree,irest,bmnmix,dorest,velsave,frcsave, + common /mdstuf/ nfree,irest,nrespa,bmnmix,dorest,velsave,frcsave, & uindsave,integrate diff --git a/source/respa.f b/source/respa.f index cff107d20..aa28c3d88 100644 --- a/source/respa.f +++ b/source/respa.f @@ -38,6 +38,7 @@ subroutine respa (istep,dt) include 'atmtyp.i' include 'atoms.i' include 'freeze.i' + include 'mdstuf.i' include 'moldyn.i' include 'units.i' include 'usage.i' @@ -62,12 +63,8 @@ subroutine respa (istep,dt) c c set some time values for the dynamics integration c - eps = 0.00000001d0 - dalt = 0.00025d0 - nalt = int(dt/(dalt+eps)) + 1 - dalt = dble(nalt) dt_2 = 0.5d0 * dt - dta = dt / dalt + dta = dt / (dble(nrespa)) dta_2 = 0.5d0 * dta c c make half-step temperature and pressure corrections @@ -102,7 +99,7 @@ subroutine respa (istep,dt) c c find fast-evolving velocities and positions via Verlet recursion c - do k = 1, nalt + do k = 1, nrespa if (use_rattle) then do i = 1, n if (use(i)) then diff --git a/source/stochisok.f b/source/stochisok.f new file mode 100755 index 000000000..f5ed1d1fa --- /dev/null +++ b/source/stochisok.f @@ -0,0 +1,467 @@ +c +c +c ################################################### +c ## COPYRIGHT (C) 2011 by Jay William Ponder ## +c ## All Rights Reserved ## +c ################################################### +c +c ############################################################# +c ## ## +c ## subroutine stochisok ## +c ## ## +c ############################################################# +c +c +c + subroutine stochisok (istep,dt) + implicit none + include 'sizes.i' + include 'atmtyp.i' + include 'freeze.i' + include 'moldyn.i' + include 'units.i' + include 'virial.i' + include 'atoms.i' + include 'bath.i' + include 'bond.i' + include 'bound.i' + include 'inform.i' + include 'iounit.i' + include 'keys.i' + include 'mdstuf.i' + include 'potent.i' + include 'solute.i' + include 'stodyn.i' + include 'usage.i' + integer i,j,k + integer istep + real*8 dt,dt_2 + real*8 dta,dta_2 + real*8 epot,etot + real*8 eksum,eps + real*8 temp,pres + real*8 ealt,dalt + real*8 kBT,LkT,tempstor,len_fac + real*8 ekin(3,3) + real*8 stress(3,3) + real*8 viralt(3,3) + real*8, allocatable :: xold(:) + real*8, allocatable :: yold(:) + real*8, allocatable :: zold(:) + real*8, allocatable :: derivs(:,:) + + +c +c +c set some values for the dynamics integration +c + dt_2 = 0.5d0 * dt + dta = dt / (dble(nrespa)) + dta_2 = 0.5d0 * dta + + kBT = boltzmann*kelvin + LkT = len_nhc*kBT + len_fac=dble(len_nhc)/dble(len_nhc+1) +c +c +c perform dynamic allocation of some local arrays +c + allocate (derivs(3,n)) + allocate (xold(n)) + allocate (yold(n)) + allocate (zold(n)) + +c +c +c initialize virial from fast-evolving potential energy terms +c + do i = 1, 3 + do j = 1, 3 + viralt(j,i) = 0.0d0 + end do + end do + +c +c +c begin inner RESPA loop +c + + do k = 1, nrespa + do i = 1, n + if (use(i)) then + xold(i) = x(i) + yold(i) = y(i) + zold(i) = z(i) + end if + end do + +c +c +c nose-like isokinetic operator +c + call applyisok (dta) +c +c +c stochastic noise operator +c + call v2random (dta) +c +c +c force-like isokinetic operator +c + + call applyisok2 (dta) +c +c Evolve positions +c + if (use_rattle) then + do i = 1, n + xold(i) = x(i) + yold(i) = y(i) + zold(i) = z(i) + x(i) = x(i) + v(1,i)*dta + y(i) = y(i) + v(2,i)*dta + z(i) = z(i) + v(3,i)*dta + end do + call rattle (dta,xold,yold,zold) + else + do i = 1, n + if (use(i)) then + x(i) = x(i) + v(1,i)*dta + y(i) = y(i) + v(2,i)*dta + z(i) = z(i) + v(3,i)*dta + end if + end do + end if +c +c get the potential energy and atomic forces +c + if (k.eq.nrespa) then + call gradfast (ealt,derivs) + + if (use_rattle) then + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + call gradslow (epot,derivs) + do i = 1, n + do j = 1, 3 + a(j,i) = a(j,i) - ( nrespa * convert + & * derivs(j,i) / mass(i) ) + end do + end do + call rattle2 (dta) + else + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + call gradslow (epot,derivs) + do i = 1, n + do j = 1, 3 + a(j,i) = a(j,i) - ( nrespa * convert + & * derivs(j,i) / mass(i) ) + end do + end do + end if + + else + call gradfast (ealt,derivs) + if (use_rattle) then + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + call rattle2 (dta) + else + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + end if + end if + +c +c +c force-like isokinetic operator +c + call applyisok2 (dta) +c +c +c nose-like isokinetic operator +c + call applyisok (dta) + + end do + +c +c find the constraint-corrected full-step velocities +c + if (use_rattle) call rattle2 (dt) +c +c total potential and virial from sum of fast and slow parts +c + epot = epot + ealt + do i = 1, 3 + do j = 1, 3 + vir(j,i) = vir(j,i) + viralt(j,i) + end do + end do + + +c +c perform deallocation of some local arrays +c + + deallocate (derivs) + deallocate (xold) + deallocate (yold) + deallocate (zold) +c +c total energy is sum of kinetic and potential energies +c + etot = eksum + epot + + call kinetic (eksum,ekin) + call pressure (dt,epot,ekin,temp,pres,stress) + temp = 2.0d0 * eksum / (dble(nfree) * gasconst) + +c +c compute statistics and save trajectory for this step +c + call mdstat (istep,dt,etot,epot,eksum,temp,pres) + call mdsave (istep,dt,epot) + call mdrest (istep) + return + end + + +c +c ################################################################## +c ## ## +c ## subroutine applyisok ## +c ## ## +c ################################################################## +c +c +c "applyisok" represents the Liouville operator iL_N in Leimkuhler, Margul, and Tuckerman (2013) +c +c + subroutine applyisok (dt) + implicit none + include 'sizes.i' + include 'atmtyp.i' + include 'atoms.i' + include 'bath.i' + include 'freeze.i' + include 'moldyn.i' + include 'units.i' + include 'usage.i' + include 'virial.i' + integer nc,ns,i,j,k,iatom,inhc,ixyz + real*8 dt,dtc,dts + real*8 dt2,dt4,dt8 + real*8 tempstor,kBT,LkT + real*8 w(3) + real*8 len_fac + + len_fac = dble(len_nhc)/dble(len_nhc+1) + + kBT = boltzmann*kelvin + LkT = (dble(len_nhc))*kBT + + nc = 5 + ns = 3 + dtc = dt / dble(nc) + w(1) = 1.0d0 / (2.0d0-2.0d0**(1.0d0/3.0d0)) + w(2) = 1.0d0 - 2.0d0*w(1) + w(3) = w(1) + +c +c use multiple time steps and Suzuki-Yoshida decomposition +c + do k = 1, nc + do j = 1, ns + + dts = w(j) * dtc + dt2 = 0.5d0 * dts + dt4 = 0.25d0 * dts + dt8 = 0.125d0 * dts + + do iatom = 1, n + do inhc = 1, len_nhc + do ixyz = 1, 3 + tempstor = q_iso1(inhc,ixyz,iatom) + & *v_iso1(inhc,ixyz,iatom)*v_iso1(inhc,ixyz,iatom) + & - kBT + v_iso2(inhc,ixyz,iatom) = v_iso2(inhc,ixyz,iatom) + & + dt4*tempstor/q_iso2(inhc,ixyz,iatom) + end do + end do + end do + + do iatom = 1, n + do ixyz = 1, 3 + tempstor = mass(iatom)*v(ixyz,iatom)*v(ixyz,iatom) + do inhc = 1, len_nhc + tempstor = tempstor + (len_fac * + & q_iso1(inhc,ixyz,iatom) * + & v_iso1(inhc,ixyz,iatom) * + & v_iso1(inhc,ixyz,iatom) * + & exp(-2.0 * dt2 * v_iso2(inhc,ixyz,iatom))) + end do + tempstor = sqrt(LkT/tempstor) + v(ixyz,iatom) = v(ixyz,iatom) * tempstor + do inhc = 1, len_nhc + v_iso1(inhc,ixyz,iatom)=v_iso1(inhc,ixyz,iatom) * + & tempstor*exp(-dt2*v_iso2(inhc,ixyz,iatom)) + end do + end do + end do + + do iatom = 1, n + do inhc = 1, len_nhc + do ixyz = 1, 3 + tempstor = q_iso1(inhc,ixyz,iatom) + & *v_iso1(inhc,ixyz,iatom)*v_iso1(inhc,ixyz,iatom) + & - kBT + v_iso2(inhc,ixyz,iatom) = v_iso2(inhc,ixyz,iatom) + & + dt4*tempstor/q_iso2(inhc,ixyz,iatom) + end do + end do + end do + + end do + end do +c do i = 1, 1 +c do j = 1, 1 +c write (*,10) v(j,i) +c do k = 1, len_nhc +c write (*,10) v_iso1(k,j,i) +c end do +c end do +c end do +10 format (3f16.6) + return + end +c +c ################################################################## +c ## ## +c ## subroutine applyisok2 ## +c ## ## +c ################################################################## +c +c +c "applyisok2" represents the Liuville operator iL_v in Leimkuhler, Margul, and Tuckerman (2013) +c +c + subroutine applyisok2 (dt) + implicit none + include 'sizes.i' + include 'atmtyp.i' + include 'atoms.i' + include 'bath.i' + include 'freeze.i' + include 'moldyn.i' + include 'units.i' + include 'usage.i' + include 'virial.i' + real*8 kBT,LkT,isoA,isoB + real*8 dt,dt_2,rootB,arg + real*8 iso_s,iso_sdot + real*8 epot,etot + integer iatom,inhc,ixyz,i,j,k + + dt_2 = 0.5d0*dt + kBT = boltzmann*kelvin + LkT = (dble(len_nhc))*kBT + + + do iatom = 1, n + do ixyz = 1, 3 + isoA = mass(iatom)*a(ixyz,iatom)*v(ixyz,iatom)/LkT + isoB = mass(iatom)*a(ixyz,iatom)*a(ixyz,iatom)/LkT + rootB = sqrt(isoB) + arg = dt_2*rootB +c write (*,10) isoA +c write (*,10) isoB +c write (*,10) LkT + if (arg > 0.00001) then + iso_s = (1.0/rootB)*sinh(arg) + + & (isoA/isoB)*(cosh(arg)-1.0) + iso_sdot = cosh(arg) + (isoA/rootB)*sinh(arg) + else + iso_s = ((((isoB*isoA/24.0)*dt_2 +isoB/6.0)*dt_2 + + & 0.5*isoA)*dt_2 + 1.0)*dt_2 + iso_sdot = (((isoB*isoA/6.0)*dt_2 +0.5*isoB)*dt_2 + + & isoA)*dt_2 + 1.0 + end if + v(ixyz,iatom) = (v(ixyz,iatom) + (a(ixyz,iatom))*iso_s) + & /iso_sdot + do inhc = 1, len_nhc + v_iso1(inhc,ixyz,iatom)=v_iso1(inhc,ixyz,iatom)/iso_sdot + end do + end do + end do + + + return + end +c +c ################################################################## +c ## ## +c ## subroutine v2random ## +c ## ## +c ################################################################## +c +c +c "v2random" represents the Liuville operator iL_OU in Leimkuhler, Margul, and Tuckerman (2013) +c +c + subroutine v2random (dt) + implicit none + include 'sizes.i' + include 'atmtyp.i' + include 'atoms.i' + include 'bath.i' + include 'freeze.i' + include 'moldyn.i' + include 'units.i' + include 'usage.i' + include 'virial.i' + real*8 dt,dt_2,kBT + real*8 gamma,stoch_e,sigma + real*8 normal + integer iatom,inhc,ixyz + + kBT = boltzmann*kelvin + gamma = 1.0d0 / (200.0d0*dt) + sigma = sqrt(kBT* (1.0d0-exp(-2.0d0*gamma*dt) )/q_iso2(1,1,1) ) + stoch_e = exp(-gamma*dt) + + do iatom = 1, n + do ixyz = 1, 3 + do inhc = 1, len_nhc + v_iso2(inhc,ixyz,iatom) = v_iso2(inhc,ixyz,iatom) * + & stoch_e + v_iso2(inhc,ixyz,iatom) = v_iso2(inhc,ixyz,iatom) + + & (sigma * normal () ) + end do + end do + end do + + return + end \ No newline at end of file