From f861b654f6ea3bc81576cc4653c4d949946c93f6 Mon Sep 17 00:00:00 2001 From: "Daniel T. Margul" Date: Fri, 7 Feb 2014 01:31:58 -0500 Subject: [PATCH 1/8] Single timestep stochastic isokinetic dynamics --- source/bath.i | 25 ++- source/dynamic.f | 9 +- source/mdinit.f | 84 ++++++++- source/stochisok.f | 446 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 546 insertions(+), 18 deletions(-) mode change 100644 => 100755 source/bath.i mode change 100644 => 100755 source/dynamic.f mode change 100644 => 100755 source/mdinit.f create mode 100755 source/stochisok.f diff --git a/source/bath.i b/source/bath.i old mode 100644 new mode 100755 index f06f93ef5..76c7b30c9 --- a/source/bath.i +++ b/source/bath.i @@ -36,23 +36,32 @@ 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 -c - integer maxnose - parameter (maxnose=4) - integer voltrial + real*8 q_isok,v_isok + 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), + & 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,isok_L,isok_M, + & 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..f89459986 --- 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,13 @@ 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 + write (iout,410) + 410 format (/,' Molecular Dynamics Trajectory via', & ' Modified Beeman Algorithm') end if c @@ -284,6 +289,8 @@ program dynamic call rgdstep (istep,dt) else if (integrate .eq. 'RESPA') then call respa (istep,dt) + else if (integrate .eq. 'STOCH-ISOK') then + call stochisok (istep,dt) else call beeman (istep,dt) end if diff --git a/source/mdinit.f b/source/mdinit.f old mode 100644 new mode 100755 index ff73528a2..6764fb20a --- 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 @@ -90,6 +91,7 @@ subroutine mdinit voltrial = 20 volmove = 100.0d0 volscale = 'ATOMIC' + c c check for keywords containing any altered parameters c @@ -147,9 +149,18 @@ subroutine mdinit call upcase (volscale) else if (keyword(1:9) .eq. 'PRINTOUT ') then read (string,*,err=10,end=10) iprint + else if (keyword(1:6) .eq. 'ISOKL ') then + read (string,*,err=10,end=10) isok_L + else if (keyword(1:6) .eq. 'ISOKM ') then + read (string,*,err=10,end=10) isok_M + else if (keyword(1:11) .eq. 'NHC-LENGTH ') then + read (string,*,err=10,end=10) len_nhc + end if 10 continue end do + + c c make sure all atoms or groups have a nonzero mass c @@ -167,7 +178,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 +372,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 +383,73 @@ 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 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 - 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 + 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 + write (iout,94) tempstor/LkT 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 +457,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 +476,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/stochisok.f b/source/stochisok.f new file mode 100755 index 000000000..985a1c87c --- /dev/null +++ b/source/stochisok.f @@ -0,0 +1,446 @@ +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 + integer nalt + real*8 dt,dt_2,tempstor + 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(:,:) +c +c +c set some time values for the dynamics integration +c + dt_2 = 0.5d0 * dt + 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 save atom positions +c + 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 print some velocities and v_1j +c + +c do i = 1, 2 +c do j = 1, 3 +c write (iout,10) v(j,i) +c do k = 1, len_nhc +c write (iout,10) v_iso1(k,j,i) +c end do +c end do +c end do +10 format (3f16.6) + + do i = 1, 1 + do j = 1, 1 + 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 + write (iout,10) tempstor/LkT + end do + end do + +c +c {v1-v2} evolution and Nose-like isokinetic term +c + call applyisok (dt_2) +c +c Random noise +c + call v2random (dt) +c +c Force-like isokinetic term +c + call applyisok2 (dt) +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)*dt + y(i) = y(i) + v(2,i)*dt + z(i) = z(i) + v(3,i)*dt + end do + call rattle (dta,xold,yold,zold) + else + do i = 1, n + if (use(i)) then + x(i) = x(i) + v(1,i)*dt + y(i) = y(i) + v(2,i)*dt + z(i) = z(i) + v(3,i)*dt + end if + end do + end if + +c +c get the potential energy and atomic forces +c + call gradient (epot,derivs) + +c +c use Newton's second law to get the next accelerations +c + do i = 1, n + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end do +c +c Force-like isokinetic term +c + call applyisok2 (dt) +c +c {v1-v2} evolution and Nose-like isokinetic term +c + call applyisok (dt_2) + +c +c perform deallocation of some local arrays +c + deallocate (xold) + deallocate (yold) + deallocate (zold) + deallocate (derivs) +c +c find the constraint-corrected full-step velocities +c + if (use_rattle) call rattle2 (dt) + + call kinetic (eksum,ekin) + 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 Liuville 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 = len_nhc*kBT +c write (*,40) dt +40 format (3f16.6) + + + 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 = 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 + +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 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.0 / (200.0*dt) + sigma = sqrt(kBT* (1.0-exp(-2.0*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 +c +c ################################################################## +c ## ## +c ## subroutine checkisok ## +c ## ## +c ################################################################## +c +c +c +c + subroutine checkisok + 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.0 / (20.0*dt) + sigma = sqrt(kBT* (1.0-exp(-2.0*gamma*dt) )/q_iso2(1,1,1) ) + stoch_e = exp(-gamma*dt) + + return + end \ No newline at end of file From 75dfb85335dfc0b34eab6d12286491aa541c75cb Mon Sep 17 00:00:00 2001 From: "Daniel T. Margul" Date: Fri, 7 Feb 2014 10:17:18 -0500 Subject: [PATCH 2/8] Single time step stochastic isokinetic dynamics. REQUIRED keywords integrator = STOCH-ISOK, nhc-length = [1-4] --- source/bath.i | 25 ++- source/dynamic.f | 9 +- source/mdinit.f | 84 ++++++++- source/stochisok.f | 422 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 522 insertions(+), 18 deletions(-) mode change 100644 => 100755 source/bath.i mode change 100644 => 100755 source/dynamic.f mode change 100644 => 100755 source/mdinit.f create mode 100755 source/stochisok.f diff --git a/source/bath.i b/source/bath.i old mode 100644 new mode 100755 index f06f93ef5..76c7b30c9 --- a/source/bath.i +++ b/source/bath.i @@ -36,23 +36,32 @@ 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 -c - integer maxnose - parameter (maxnose=4) - integer voltrial + real*8 q_isok,v_isok + 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), + & 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,isok_L,isok_M, + & 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..f89459986 --- 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,13 @@ 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 + write (iout,410) + 410 format (/,' Molecular Dynamics Trajectory via', & ' Modified Beeman Algorithm') end if c @@ -284,6 +289,8 @@ program dynamic call rgdstep (istep,dt) else if (integrate .eq. 'RESPA') then call respa (istep,dt) + else if (integrate .eq. 'STOCH-ISOK') then + call stochisok (istep,dt) else call beeman (istep,dt) end if diff --git a/source/mdinit.f b/source/mdinit.f old mode 100644 new mode 100755 index ff73528a2..6764fb20a --- 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 @@ -90,6 +91,7 @@ subroutine mdinit voltrial = 20 volmove = 100.0d0 volscale = 'ATOMIC' + c c check for keywords containing any altered parameters c @@ -147,9 +149,18 @@ subroutine mdinit call upcase (volscale) else if (keyword(1:9) .eq. 'PRINTOUT ') then read (string,*,err=10,end=10) iprint + else if (keyword(1:6) .eq. 'ISOKL ') then + read (string,*,err=10,end=10) isok_L + else if (keyword(1:6) .eq. 'ISOKM ') then + read (string,*,err=10,end=10) isok_M + else if (keyword(1:11) .eq. 'NHC-LENGTH ') then + read (string,*,err=10,end=10) len_nhc + end if 10 continue end do + + c c make sure all atoms or groups have a nonzero mass c @@ -167,7 +178,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 +372,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 +383,73 @@ 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 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 - 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 + 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 + write (iout,94) tempstor/LkT 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 +457,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 +476,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/stochisok.f b/source/stochisok.f new file mode 100755 index 000000000..31007b9da --- /dev/null +++ b/source/stochisok.f @@ -0,0 +1,422 @@ +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 + integer nalt + real*8 dt,dt_2,tempstor + 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(:,:) +c +c +c set some time values for the dynamics integration +c + dt_2 = 0.5d0 * dt + 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 save atom positions +c + 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 {v1-v2} evolution and Nose-like isokinetic term +c + call applyisok (dt) +c +c Random noise +c + call v2random (dt) +c +c Force-like isokinetic term +c + call applyisok2 (dt) +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)*dt + y(i) = y(i) + v(2,i)*dt + z(i) = z(i) + v(3,i)*dt + end do + call rattle (dta,xold,yold,zold) + else + do i = 1, n + if (use(i)) then + x(i) = x(i) + v(1,i)*dt + y(i) = y(i) + v(2,i)*dt + z(i) = z(i) + v(3,i)*dt + end if + end do + end if + +c +c get the potential energy and atomic forces +c + call gradient (epot,derivs) + +c +c use Newton's second law to get the next accelerations +c + do i = 1, n + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end do +c +c Force-like isokinetic term +c + call applyisok2 (dt) +c +c {v1-v2} evolution and Nose-like isokinetic term +c + call applyisok (dt) + +c +c perform deallocation of some local arrays +c + deallocate (xold) + deallocate (yold) + deallocate (zold) + deallocate (derivs) +c +c find the constraint-corrected full-step velocities +c + if (use_rattle) call rattle2 (dt) + + call kinetic (eksum,ekin) + temp = 2.0d0 * eksum / (dble(nfree) * gasconst) + etot = eksum + epot +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 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 = len_nhc*kBT +c write (*,40) dt +40 format (3f16.6) + + + 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 = 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 + +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 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.0 / (200.0*dt) + sigma = sqrt(kBT* (1.0-exp(-2.0*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 +c +c ################################################################## +c ## ## +c ## subroutine checkisok ## +c ## ## +c ################################################################## +c +c +c +c + subroutine checkisok + 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.0 / (20.0*dt) + sigma = sqrt(kBT* (1.0-exp(-2.0*gamma*dt) )/q_iso2(1,1,1) ) + stoch_e = exp(-gamma*dt) + + return + end From 7ab9c25900de803a2a9b481c1e9d22164c755e81 Mon Sep 17 00:00:00 2001 From: dmargul Date: Fri, 7 Feb 2014 10:22:43 -0500 Subject: [PATCH 3/8] A few changes Cleaned up comments, and fixed calls to applyisok (dt_2) which are now applyisok (dt). --- source/stochisok.f | 90 ++++------------------------------------------ 1 file changed, 6 insertions(+), 84 deletions(-) diff --git a/source/stochisok.f b/source/stochisok.f index 985a1c87c..0173b6098 100755 --- a/source/stochisok.f +++ b/source/stochisok.f @@ -76,35 +76,11 @@ subroutine stochisok (istep,dt) end if end do -c -c print some velocities and v_1j -c - -c do i = 1, 2 -c do j = 1, 3 -c write (iout,10) v(j,i) -c do k = 1, len_nhc -c write (iout,10) v_iso1(k,j,i) -c end do -c end do -c end do -10 format (3f16.6) - - do i = 1, 1 - do j = 1, 1 - 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 - write (iout,10) tempstor/LkT - end do - end do c c {v1-v2} evolution and Nose-like isokinetic term c - call applyisok (dt_2) + call applyisok (dt) c c Random noise c @@ -157,7 +133,7 @@ subroutine stochisok (istep,dt) c c {v1-v2} evolution and Nose-like isokinetic term c - call applyisok (dt_2) + call applyisok (dt) c c perform deallocation of some local arrays @@ -173,7 +149,7 @@ subroutine stochisok (istep,dt) call kinetic (eksum,ekin) temp = 2.0d0 * eksum / (dble(nfree) * gasconst) - + etot = eksum + epot c c compute statistics and save trajectory for this step c @@ -216,8 +192,6 @@ subroutine applyisok (dt) kBT = boltzmann*kelvin LkT = len_nhc*kBT -c write (*,40) dt -40 format (3f16.6) nc = 5 @@ -283,15 +257,7 @@ subroutine applyisok (dt) 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 @@ -333,9 +299,7 @@ subroutine applyisok2 (dt) 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) @@ -346,6 +310,7 @@ subroutine applyisok2 (dt) 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 @@ -354,16 +319,6 @@ subroutine applyisok2 (dt) 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 @@ -411,36 +366,3 @@ subroutine v2random (dt) return end -c -c ################################################################## -c ## ## -c ## subroutine checkisok ## -c ## ## -c ################################################################## -c -c -c -c - subroutine checkisok - 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.0 / (20.0*dt) - sigma = sqrt(kBT* (1.0-exp(-2.0*gamma*dt) )/q_iso2(1,1,1) ) - stoch_e = exp(-gamma*dt) - - return - end \ No newline at end of file From 285d9b2b6bb36f69dcbfeb5f7b9cd26167fa8992 Mon Sep 17 00:00:00 2001 From: "Daniel T. Margul" Date: Fri, 7 Feb 2014 11:33:00 -0500 Subject: [PATCH 4/8] Added variable nrespa to mdstuf.i: number of small time steps per large time step in RESPA schemes --- source/mdstuf.i | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) 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 From 6797cc5c4554c02975b2be48a60664a6eabc7062 Mon Sep 17 00:00:00 2001 From: "Daniel T. Margul" Date: Fri, 7 Feb 2014 11:50:46 -0500 Subject: [PATCH 5/8] Modified mdinit.f and respa.f so that any small or large time step is possible with the keyword nrespa in .key file. nrespa is the number of small time steps per large time step. --- source/mdinit.f | 3 ++- source/respa.f | 9 +++------ 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/source/mdinit.f b/source/mdinit.f index 6764fb20a..7c8edc2ef 100755 --- a/source/mdinit.f +++ b/source/mdinit.f @@ -155,7 +155,8 @@ subroutine mdinit read (string,*,err=10,end=10) isok_M 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 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 From bc825657c2c82bea78b36d1740b44892eddf4d55 Mon Sep 17 00:00:00 2001 From: "Daniel T. Margul" Date: Fri, 7 Feb 2014 11:56:03 -0500 Subject: [PATCH 6/8] cleaned up stochisok.f --- source/stochisok.f | 60 +++------------------------------------------- 1 file changed, 3 insertions(+), 57 deletions(-) diff --git a/source/stochisok.f b/source/stochisok.f index 31007b9da..8e338fefc 100755 --- a/source/stochisok.f +++ b/source/stochisok.f @@ -192,8 +192,6 @@ subroutine applyisok (dt) kBT = boltzmann*kelvin LkT = len_nhc*kBT -c write (*,40) dt -40 format (3f16.6) nc = 5 @@ -259,15 +257,7 @@ subroutine applyisok (dt) 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 @@ -309,9 +299,7 @@ subroutine applyisok2 (dt) 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) @@ -322,6 +310,7 @@ subroutine applyisok2 (dt) 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 @@ -330,16 +319,6 @@ subroutine applyisok2 (dt) 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 @@ -387,36 +366,3 @@ subroutine v2random (dt) return end -c -c ################################################################## -c ## ## -c ## subroutine checkisok ## -c ## ## -c ################################################################## -c -c -c -c - subroutine checkisok - 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.0 / (20.0*dt) - sigma = sqrt(kBT* (1.0-exp(-2.0*gamma*dt) )/q_iso2(1,1,1) ) - stoch_e = exp(-gamma*dt) - - return - end From f0cc122322951f37fac75b29f08b2eb9fce05366 Mon Sep 17 00:00:00 2001 From: "Daniel T. Margul" Date: Tue, 18 Feb 2014 22:48:58 -0500 Subject: [PATCH 7/8] Isokinetic and stochastic isokinetic MD schemes now support multiple time steps. --- source/bath.i | 13 ++- source/dynamic.f | 8 +- source/mdinit.f | 63 +++++++++++-- source/stochisok.f | 226 ++++++++++++++++++++++++++++++++------------- 4 files changed, 236 insertions(+), 74 deletions(-) diff --git a/source/bath.i b/source/bath.i index 76c7b30c9..f9fddb7be 100755 --- a/source/bath.i +++ b/source/bath.i @@ -35,8 +35,12 @@ 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 - real*8 q_isok,v_isok real*8 q_iso1,q_iso2,v_iso1,v_iso2 real*8 vnh,qnh,gnh real*8 kelvin,atmsph @@ -60,8 +64,9 @@ c & 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,isok_L,isok_M, - & len_nhc,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 index f89459986..cd3527a3f 100755 --- a/source/dynamic.f +++ b/source/dynamic.f @@ -266,9 +266,13 @@ program dynamic write (iout,400) 400 format (/,' Molecular Dynamics Trajectory via', & ' Stochastic Isokinetic Algorithm') - else + 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 @@ -289,6 +293,8 @@ 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 diff --git a/source/mdinit.f b/source/mdinit.f index 7c8edc2ef..a2b9eac7a 100755 --- a/source/mdinit.f +++ b/source/mdinit.f @@ -69,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 @@ -148,11 +150,7 @@ 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 - else if (keyword(1:6) .eq. 'ISOKL ') then - read (string,*,err=10,end=10) isok_L - else if (keyword(1:6) .eq. 'ISOKM ') then - read (string,*,err=10,end=10) isok_M + 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 @@ -399,7 +397,61 @@ subroutine mdinit 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 + 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, @@ -442,7 +494,6 @@ subroutine mdinit 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 - write (iout,94) tempstor/LkT end do end do deallocate (derivs) diff --git a/source/stochisok.f b/source/stochisok.f index 064a4f4a3..f5ed1d1fa 100755 --- a/source/stochisok.f +++ b/source/stochisok.f @@ -33,16 +33,15 @@ subroutine stochisok (istep,dt) include 'solute.i' include 'stodyn.i' include 'usage.i' - integer i,j,k integer istep - integer nalt - real*8 dt,dt_2,tempstor - real*8 dta,dta_2,kBT,LkT + 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) @@ -50,106 +49,196 @@ subroutine stochisok (istep,dt) real*8, allocatable :: yold(:) real*8, allocatable :: zold(:) real*8, allocatable :: derivs(:,:) + + c c -c set some time values for the dynamics integration +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)) - allocate (derivs(3,n)) c -c save atom positions c - do i = 1, n - if (use(i)) then - xold(i) = x(i) - yold(i) = y(i) - zold(i) = z(i) - end if +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 {v1-v2} evolution and Nose-like isokinetic term c - call applyisok (dt) +c nose-like isokinetic operator c -c Random noise + call applyisok (dta) c - call v2random (dt) c -c Force-like isokinetic term +c stochastic noise operator +c + call v2random (dta) c - call applyisok2 (dt) c -c Evolve positions +c force-like isokinetic operator 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)*dt - y(i) = y(i) + v(2,i)*dt - z(i) = z(i) + v(3,i)*dt - end do - call rattle (dta,xold,yold,zold) - else - do i = 1, n - if (use(i)) then - x(i) = x(i) + v(1,i)*dt - y(i) = y(i) + v(2,i)*dt - z(i) = z(i) + v(3,i)*dt + 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 - end do - end if - c c get the potential energy and atomic forces c - call gradient (epot,derivs) + 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 use Newton's second law to get the next accelerations c - do i = 1, n - do j = 1, 3 - a(j,i) = -convert * derivs(j,i) / mass(i) - end do +c force-like isokinetic operator +c + call applyisok2 (dta) +c +c +c nose-like isokinetic operator +c + call applyisok (dta) + end do + c -c Force-like isokinetic term +c find the constraint-corrected full-step velocities c - call applyisok2 (dt) + if (use_rattle) call rattle2 (dt) c -c {v1-v2} evolution and Nose-like isokinetic term +c total potential and virial from sum of fast and slow parts c - call applyisok (dt) + 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) - deallocate (derivs) + deallocate (zold) c -c find the constraint-corrected full-step velocities +c total energy is sum of kinetic and potential energies c - if (use_rattle) call rattle2 (dt) + etot = eksum + epot call kinetic (eksum,ekin) + call pressure (dt,epot,ekin,temp,pres,stress) temp = 2.0d0 * eksum / (dble(nfree) * gasconst) - etot = eksum + epot + c c compute statistics and save trajectory for this step c @@ -159,6 +248,7 @@ subroutine stochisok (istep,dt) return end + c c ################################################################## c ## ## @@ -167,7 +257,7 @@ subroutine stochisok (istep,dt) c ################################################################## c c -c "applyisok" represents the Liuville operator iL_N in Leimkuhler, Margul, and Tuckerman (2013) +c "applyisok" represents the Liouville operator iL_N in Leimkuhler, Margul, and Tuckerman (2013) c c subroutine applyisok (dt) @@ -191,8 +281,7 @@ subroutine applyisok (dt) len_fac = dble(len_nhc)/dble(len_nhc+1) kBT = boltzmann*kelvin - LkT = len_nhc*kBT - + LkT = (dble(len_nhc))*kBT nc = 5 ns = 3 @@ -257,6 +346,15 @@ subroutine applyisok (dt) 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 @@ -289,7 +387,7 @@ subroutine applyisok2 (dt) dt_2 = 0.5d0*dt kBT = boltzmann*kelvin - LkT = len_nhc*kBT + LkT = (dble(len_nhc))*kBT do iatom = 1, n @@ -298,7 +396,9 @@ subroutine applyisok2 (dt) 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) @@ -309,7 +409,6 @@ subroutine applyisok2 (dt) 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 @@ -318,6 +417,7 @@ subroutine applyisok2 (dt) end do end do + return end c @@ -348,8 +448,8 @@ subroutine v2random (dt) integer iatom,inhc,ixyz kBT = boltzmann*kelvin - gamma = 1.0 / (200.0*dt) - sigma = sqrt(kBT* (1.0-exp(-2.0*gamma*dt) )/q_iso2(1,1,1) ) + 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 @@ -364,4 +464,4 @@ subroutine v2random (dt) end do return - end + end \ No newline at end of file From 15825fd524a0432acc10a154a8e2037d5c7da729 Mon Sep 17 00:00:00 2001 From: "Daniel T. Margul" Date: Fri, 11 Apr 2014 10:06:42 -0400 Subject: [PATCH 8/8] isokinetic.f was not added, fixed --- source/isokinetic.f | 480 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 480 insertions(+) create mode 100644 source/isokinetic.f 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