subroutine MD c------------------------------------------------ c The driver for molecular dynamics subroutines c------------------------------------------------ implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' double precision cm(3),L(3),vcm(3) double precision energia(0:n_ene) integer ilen,rstcount external ilen character*50 tytul common /gucio/ cm,energia integer itime c t_MDsetup=0.0d0 t_langsetup=0.0d0 t_MD=0.0d0 t_enegrad=0.0d0 t_sdsetup=0.0d0 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started" tt0 = tcpu() c Determine the inverse of the inertia matrix. call setup_MD_matrices c Initialize MD call init_MD t_MDsetup = tcpu()-tt0 rstcount=0 c Entering the MD loop tt0 = tcpu() if (lang.eq.2 .or. lang.eq.3) then call setup_fricmat if (lang.eq.2) then call sd_verlet_p_setup else call sd_verlet_ciccotti_setup endif do i=1,dimen do j=1,dimen pfric0_mat(i,j,0)=pfric_mat(i,j) afric0_mat(i,j,0)=afric_mat(i,j) vfric0_mat(i,j,0)=vfric_mat(i,j) prand0_mat(i,j,0)=prand_mat(i,j) vrand0_mat1(i,j,0)=vrand_mat1(i,j) vrand0_mat2(i,j,0)=vrand_mat2(i,j) enddo enddo flag_stoch(0)=.true. do i=1,maxflag_stoch flag_stoch(i)=.false. enddo else if (lang.eq.1 .or. lang.eq.4) then call setup_fricmat endif t_langsetup=tcpu()-tt0 tt0=tcpu() do itime=1,n_timestep rstcount=rstcount+1 if (lang.gt.0 .and. surfarea .and. & mod(itime,reset_fricmat).eq.0) then if (lang.eq.2 .or. lang.eq.3) then call setup_fricmat if (lang.eq.2) then call sd_verlet_p_setup else call sd_verlet_ciccotti_setup endif do i=1,dimen do j=1,dimen pfric0_mat(i,j,0)=pfric_mat(i,j) afric0_mat(i,j,0)=afric_mat(i,j) vfric0_mat(i,j,0)=vfric_mat(i,j) prand0_mat(i,j,0)=prand_mat(i,j) vrand0_mat1(i,j,0)=vrand_mat1(i,j) vrand0_mat2(i,j,0)=vrand_mat2(i,j) enddo enddo flag_stoch(0)=.true. do i=1,maxflag_stoch flag_stoch(i)=.false. enddo else if (lang.eq.1 .or. lang.eq.4) then call setup_fricmat endif write (iout,'(a,i10)') & "Friction matrix reset based on surface area, itime",itime endif if (reset_vel .and. tbf .and. lang.eq.0 & .and. mod(itime,count_reset_vel).eq.0) then call random_vel write(iout,'(a,f20.2)') & "Velocities reset to random values, time",totT do i=0,2*nres do j=1,3 d_t_old(j,i)=d_t(j,i) enddo enddo endif if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then call inertia_tensor call vcm_vel(vcm) do j=1,3 d_t(j,0)=d_t(j,0)-vcm(j) enddo call kinetic(EK) kinetic_T=2.0d0/(dimen*Rb)*EK scalfac=dsqrt(T_bath/kinetic_T) write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT do i=0,2*nres do j=1,3 d_t_old(j,i)=scalfac*d_t(j,i) enddo enddo endif if (lang.ne.4) then if (RESPA) then c Time-reversible RESPA algorithm c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992) call RESPA_step(itime) else c Variable time step algorithm. call velverlet_step(itime) endif else call brown_step(itime) endif if (mod(itime,ntwe).eq.0) call statout(itime) if (mod(itime,ntwx).eq.0) then write (tytul,'("time",f8.2)') totT if(mdpdb) then call pdbout(potE,tytul,ipdb) else call cartout(totT) endif endif if (rstcount.eq.1000.or.itime.eq.n_timestep) then open(irest2,file=rest2name,status='unknown') write(irest2,*) totT,EK,potE,totE do i=1,2*nres write (irest2,'(3e15.5)') (d_t(j,i),j=1,3) enddo do i=1,2*nres write (irest2,'(3e15.5)') (dc(j,i),j=1,3) enddo close(irest2) rstcount=0 endif enddo t_MD=tcpu()-tt0 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') & ' Timing ', & 'MD calculations setup:',t_MDsetup, & 'Energy & gradient evaluation:',t_enegrad, & 'Stochastic MD setup:',t_langsetup, & 'Stochastic MD step setup:',t_sdsetup, & 'MD steps:',t_MD write (iout,'(/28(1h=),a25,27(1h=))') & ' End of MD calculation ' return end c------------------------------------------------------------------------------- subroutine brown_step(itime) c------------------------------------------------ c Perform a single Euler integration step of Brownian dynamics c------------------------------------------------ implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' double precision energia(0:n_ene),zapas(MAXRES6) integer ilen,rstcount external ilen double precision stochforcvec(MAXRES6) double precision Bmat(MAXRES6,MAXRES2),Cmat(maxres2,maxres2), & Cinv(maxres2,maxres2),GBmat(MAXRES6,MAXRES2), & Tmat(MAXRES6,MAXRES2),Pmat(maxres6,maxres6),Td(maxres6), & ppvec(maxres2) common /stochcalc/ stochforcvec common /gucio/ cm, energia integer itime logical lprn /.false./,lprn1 /.false./ integer maxiter /5/ double precision difftol /1.0d-5/ nbond=nct-nnt do i=nnt,nct if (itype(i).ne.10) nbond=nbond+1 enddo c if (lprn1) then write (iout,*) "Generalized inverse of fricmat" call matout(dimen,dimen,MAXRES6,MAXRES6,fricmat) endif do i=1,dimen do j=1,nbond Bmat(i,j)=0.0d0 enddo enddo ind=3 ind1=0 do i=nnt,nct-1 ind1=ind1+1 do j=1,3 Bmat(ind+j,ind1)=dC_norm(j,i) enddo ind=ind+3 enddo do i=nnt,nct if (itype(i).ne.10) then ind1=ind1+1 do j=1,3 Bmat(ind+j,ind1)=dC_norm(j,i+nres) enddo ind=ind+3 endif enddo if (lprn1) then write (iout,*) "Matrix Bmat" call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Bmat) endif do i=1,dimen do j=1,nbond GBmat(i,j)=0.0d0 do k=1,dimen GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j) enddo enddo enddo if (lprn1) then write (iout,*) "Matrix GBmat" call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Gbmat) endif do i=1,nbond do j=1,nbond Cmat(i,j)=0.0d0 do k=1,dimen Cmat(i,j)=Cmat(i,j)+Bmat(k,i)*GBmat(k,j) enddo enddo enddo if (lprn1) then write (iout,*) "Matrix Cmat" call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat) endif call matinvert(nbond,MAXRES2,Cmat,Cinv) if (lprn1) then write (iout,*) "Matrix Cinv" call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cinv) endif do i=1,dimen do j=1,nbond Tmat(i,j)=0.0d0 do k=1,nbond Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j) enddo enddo enddo if (lprn1) then write (iout,*) "Matrix Tmat" call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Tmat) endif do i=1,dimen do j=1,dimen if (i.eq.j) then Pmat(i,j)=1.0d0 else Pmat(i,j)=0.0d0 endif do k=1,nbond Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k) enddo enddo enddo if (lprn1) then write (iout,*) "Matrix Pmat" call MATOUT(dimen,dimen,MAXRES6,MAXRES6,Pmat) endif do i=1,dimen Td(i)=0.0d0 ind=0 do k=nnt,nct-1 ind=ind+1 Td(i)=Td(i)+vbl*Tmat(i,ind) enddo do k=nnt,nct if (itype(k).ne.10) then ind=ind+1 Td(i)=Td(i)+vbldsc0(itype(k))*Tmat(i,ind) endif enddo enddo if (lprn1) then write (iout,*) "Vector Td" do i=1,dimen write (iout,'(i5,f10.5)') i,Td(i) enddo endif call stochastic_force(stochforcvec) if (lprn) then write (iout,*) "stochforcvec" do i=1,dimen write (iout,*) i,stochforcvec(i) enddo endif do j=1,3 zapas(j)=-gcart(j,0)+stochforcvec(j) d_t_work(j)=d_t(j,0) dC_work(j)=dC_old(j,0) enddo ind=3 do i=nnt,nct-1 do j=1,3 ind=ind+1 zapas(ind)=-gcart(j,i)+stochforcvec(ind) dC_work(ind)=dC_old(j,i) enddo enddo do i=nnt,nct if (itype(i).ne.10) then do j=1,3 ind=ind+1 zapas(ind)=-gxcart(j,i)+stochforcvec(ind) dC_work(ind)=dC_old(j,i+nres) enddo endif enddo if (lprn) then write (iout,*) "Initial d_t_work" do i=1,dimen write (iout,*) i,d_t_work(i) enddo endif do i=1,dimen d_t_work(i)=0.0d0 do j=1,dimen d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j) enddo enddo do i=1,dimen zapas(i)=Td(i) do j=1,dimen zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time) enddo enddo if (lprn1) then write (iout,*) "Final d_t_work and zapas" do i=1,dimen write (iout,*) i,d_t_work(i),zapas(i) enddo endif do j=1,3 d_t(j,0)=d_t_work(j) dc(j,0)=zapas(j) dc_work(j)=dc(j,0) enddo ind=3 do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_work(i) dc(j,i)=zapas(ind+j) dc_work(ind+j)=dc(j,i) enddo ind=ind+3 enddo do i=nnt,nct do j=1,3 d_t(j,i+nres)=d_t_work(ind+j) dc(j,i+nres)=zapas(ind+j) dc_work(ind+j)=dc(j,i+nres) enddo ind=ind+3 enddo if (lprn) then call chainbuild_cart write (iout,*) "Before correction for rotational lengthening" write (iout,*) "New coordinates", & " and differences between actual and standard bond lengths" ind=0 do i=nnt,nct-1 ind=ind+1 xx=vbld(i+1)-vbl write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i),j=1,3),xx enddo do i=nnt,nct if (itype(i).ne.10) then ind=ind+1 xx=vbld(i+nres)-vbldsc0(itype(i)) write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i+nres),j=1,3),xx endif enddo endif c Second correction (rotational lengthening) c do iter=1,maxiter diffmax=0.0d0 ind=0 do i=nnt,nct-1 ind=ind+1 blen2 = scalar(dc(1,i),dc(1,i)) ppvec(ind)=2*vbl**2-blen2 diffbond=dabs(vbl-dsqrt(blen2)) if (diffbond.gt.diffmax) diffmax=diffbond if (ppvec(ind).gt.0.0d0) then ppvec(ind)=dsqrt(ppvec(ind)) else ppvec(ind)=0.0d0 endif if (lprn) then write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind) endif enddo do i=nnt,nct if (itype(i).ne.10) then ind=ind+1 blen2 = scalar(dc(1,i+nres),dc(1,i+nres)) ppvec(ind)=2*vbldsc0(itype(i))**2-blen2 diffbond=dabs(vbldsc0(itype(i))-dsqrt(blen2)) if (diffbond.gt.diffmax) diffmax=diffbond if (ppvec(ind).gt.0.0d0) then ppvec(ind)=dsqrt(ppvec(ind)) else ppvec(ind)=0.0d0 endif if (lprn) then write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind) endif endif enddo if (lprn) write (iout,*) "iter",iter," diffmax",diffmax if (diffmax.lt.difftol) goto 10 do i=1,dimen Td(i)=0.0d0 do j=1,nbond Td(i)=Td(i)+ppvec(j)*Tmat(i,j) enddo enddo do i=1,dimen zapas(i)=Td(i) do j=1,dimen zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j) enddo enddo do j=1,3 dc(j,0)=zapas(j) dc_work(j)=zapas(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 dc(j,i)=zapas(ind+j) dc_work(ind+j)=zapas(ind+j) enddo ind=ind+3 enddo do i=nnt,nct if (itype(i).ne.10) then do j=1,3 dc(j,i+nres)=zapas(ind+j) dc_work(ind+j)=zapas(ind+j) enddo ind=ind+3 endif enddo c Building the chain from the newly calculated coordinates call chainbuild_cart if (large.and. mod(itime,ntwe).eq.0) then write (iout,*) "Cartesian and internal coordinates: step 1" call cartprint call intout write (iout,'(a)') "Potential forces" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3), & (-gxcart(j,i),j=1,3) enddo write (iout,'(a)') "Stochastic forces" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3), & (stochforc(j,i+nres),j=1,3) enddo write (iout,'(a)') "Velocities" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), & (d_t(j,i+nres),j=1,3) enddo endif if (lprn) then write (iout,*) "After correction for rotational lengthening" write (iout,*) "New coordinates", & " and differences between actual and standard bond lengths" ind=0 do i=nnt,nct-1 ind=ind+1 xx=vbld(i+1)-vbl write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i),j=1,3),xx enddo do i=nnt,nct if (itype(i).ne.10) then ind=ind+1 xx=vbld(i+nres)-vbldsc0(itype(i)) write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i+nres),j=1,3),xx endif enddo endif c ENDDO c write (iout,*) "Too many attempts at correcting the bonds" c stop 10 continue tt0 = tcpu() c Calculate energy and forces call zerograd call etotal(energia) potE=energia(0)-energia(20) call cartgrad totT=totT+d_time c Calculate the kinetic and total energy and the kinetic temperature call kinetic(EK) t_enegrad=t_enegrad+tcpu()-tt0 totE=EK+potE kinetic_T=2.0d0/(dimen*Rb)*EK return end c------------------------------------------------------------------------------- subroutine velverlet_step(itime) c------------------------------------------------------------------------------- c Perform a single velocity Verlet step; the time step can be rescaled if c increments in accelerations exceed the threshold c------------------------------------------------------------------------------- implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' include 'COMMON.MUCA' double precision energia(0:n_ene),vcm(3),incr(3) double precision cm(3),L(3) integer ilen,count,rstcount external ilen character*50 tytul integer maxcount_scale /20/ common /gucio/ cm, energia double precision stochforcvec(MAXRES6) common /stochcalc/ stochforcvec integer itime logical scale double precision HNose1,HNose,HNose_nh,H,vtnp(maxres6) double precision vtnp_(maxres6),vtnp_a(maxres6) c scale=.true. icount_scale=0 if (lang.eq.1) then call sddir_precalc else if (lang.eq.2 .or. lang.eq.3) then call stochastic_force(stochforcvec) endif itime_scal=0 do while (scale) icount_scale=icount_scale+1 if (icount_scale.gt.maxcount_scale) then write (iout,*) & "ERROR: too many attempts at scaling down the time step. ", & "amax=",amax,"epdrift=",epdrift, & "damax=",damax,"edriftmax=",edriftmax, & "d_time=",d_time stop endif c First step of the velocity Verlet algorithm if (lang.eq.2) then call sd_verlet1 else if (lang.eq.3) then call sd_verlet1_ciccotti else if (lang.eq.1) then call sddir_verlet1 else if (tnp1) then call tnp1_step1 else if (tnp) then call tnp_step1 else if (tnh) then call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8) do i=0,2*nres do j=1,3 d_t_old(j,i)=d_t_old(j,i)*scale_nh enddo enddo endif call verlet1 endif c Build the chain from the newly calculated coordinates call chainbuild_cart if (rattle) call rattle1 if (large.and. mod(itime,ntwe).eq.0) then write (iout,*) "Cartesian and internal coordinates: step 1" call cartprint call intout write (iout,*) "Accelerations" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), & (d_a(j,i+nres),j=1,3) enddo write (iout,*) "Velocities, step 1" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), & (d_t(j,i+nres),j=1,3) enddo endif tt0 = tcpu() c Calculate energy and forces call zerograd E_old=potE call etotal(energia) potE=energia(0)-energia(20) call cartgrad c Get the new accelerations call lagrangian t_enegrad=t_enegrad+tcpu()-tt0 c Determine maximum acceleration and scale down the timestep if needed call max_accel call predict_edrift(epdrift) if (amax.gt.damax .or. epdrift.gt.edriftmax) then c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step scale=.true. ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) & /dlog(2.0d0)+1 itime_scal=itime_scal+ifac_time c fac_time=dmin1(damax/amax,0.5d0) fac_time=0.5d0**ifac_time d_time=d_time*fac_time if (lang.eq.2 .or. lang.eq.3) then c write (iout,*) "Calling sd_verlet_setup: 1" c Rescale the stochastic forces and recalculate or restore c the matrices of tinker integrator if (itime_scal.gt.maxflag_stoch) then if (large) write (iout,'(a,i5,a)') & "Calculate matrices for stochastic step;", & " itime_scal ",itime_scal if (lang.eq.2) then call sd_verlet_p_setup else call sd_verlet_ciccotti_setup endif write (iout,'(2a,i3,a,i3,1h.)') & "Warning: cannot store matrices for stochastic", & " integration because the index",itime_scal, & " is greater than",maxflag_stoch write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct", & " integration Langevin algorithm for better efficiency." else if (flag_stoch(itime_scal)) then if (large) write (iout,'(a,i5,a,l1)') & "Restore matrices for stochastic step; itime_scal ", & itime_scal," flag ",flag_stoch(itime_scal) do i=1,dimen do j=1,dimen pfric_mat(i,j)=pfric0_mat(i,j,itime_scal) afric_mat(i,j)=afric0_mat(i,j,itime_scal) vfric_mat(i,j)=vfric0_mat(i,j,itime_scal) prand_mat(i,j)=prand0_mat(i,j,itime_scal) vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal) vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal) enddo enddo else if (large) write (iout,'(2a,i5,a,l1)') & "Calculate & store matrices for stochastic step;", & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal) if (lang.eq.2) then call sd_verlet_p_setup else call sd_verlet_ciccotti_setup endif flag_stoch(ifac_time)=.true. do i=1,dimen do j=1,dimen pfric0_mat(i,j,itime_scal)=pfric_mat(i,j) afric0_mat(i,j,itime_scal)=afric_mat(i,j) vfric0_mat(i,j,itime_scal)=vfric_mat(i,j) prand0_mat(i,j,itime_scal)=prand_mat(i,j) vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j) vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j) enddo enddo endif fac_time=1.0d0/dsqrt(fac_time) do i=1,dimen stochforcvec(i)=fac_time*stochforcvec(i) enddo else if (lang.eq.1) then c Rescale the accelerations due to stochastic forces fac_time=1.0d0/dsqrt(fac_time) do i=1,dimen d_as_work(i)=d_as_work(i)*fac_time enddo endif if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') & "itime",itime," Timestep scaled down to ", & d_time," ifac_time",ifac_time," itime_scal",itime_scal else c Second step of the velocity Verlet algorithm if (lang.eq.2) then call sd_verlet2 else if (lang.eq.3) then call sd_verlet2_ciccotti else if (lang.eq.1) then call sddir_verlet2 else if (tnp1) then call tnp1_step2 else if (tnp) then call tnp_step2 else call verlet2 if (tnh) then call kinetic(EK) call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8) do i=0,2*nres do j=1,3 d_t(j,i)=d_t(j,i)*scale_nh enddo enddo endif endif if (rattle) call rattle2 totT=totT+d_time if (d_time.ne.d_time0) then d_time=d_time0 if (lang.eq.2 .or. lang.eq.3) then if (large) write (iout,'(a)') & "Restore original matrices for stochastic step" c write (iout,*) "Calling sd_verlet_setup: 2" c Restore the matrices of tinker integrator if the time step has been restored do i=1,dimen do j=1,dimen pfric_mat(i,j)=pfric0_mat(i,j,0) afric_mat(i,j)=afric0_mat(i,j,0) vfric_mat(i,j)=vfric0_mat(i,j,0) prand_mat(i,j)=prand0_mat(i,j,0) vrand_mat1(i,j)=vrand0_mat1(i,j,0) vrand_mat2(i,j)=vrand0_mat2(i,j,0) enddo enddo endif endif scale=.false. endif enddo c Calculate the kinetic and the total energy and the kinetic temperature if (tnp .or. tnp1) then do i=0,2*nres do j=1,3 d_t_old(j,i)=d_t(j,i) d_t(j,i)=d_t(j,i)/s_np enddo enddo endif call kinetic(EK) totE=EK+potE c diagnostics c call kinetic1(EK1) c write (iout,*) "step",itime," EK",EK," EK1",EK1 c end diagnostics c Couple the system to Berendsen bath if needed if (tbf .and. lang.eq.0) then call verlet_bath endif kinetic_T=2.0d0/(dimen*Rb)*EK c Backup the coordinates, velocities, and accelerations do i=0,2*nres do j=1,3 dc_old(j,i)=dc(j,i) if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i) d_a_old(j,i)=d_a(j,i) enddo enddo c if (mod(itime,ntwe).eq.0 .and. large) then if (mod(itime,ntwe).eq.0) then if(tnp .or. tnp1) then HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen) H=(HNose1-H0)*s_np cd write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0 cd & ,EK+potE+pi_np**2/(2*Q_np)+dimen*0.001986d0*t_bath*log(s_np) write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0 endif if(tnh) then HNose1=Hnose_nh(EK,potE) H=HNose1-H0 write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0 endif if (large) then itnp=0 do j=1,3 itnp=itnp+1 vtnp(itnp)=d_t(j,0) enddo do i=nnt,nct-1 do j=1,3 itnp=itnp+1 vtnp(itnp)=d_t(j,i) enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 itnp=itnp+1 vtnp(itnp)=d_t(j,inres) enddo endif enddo c Transform velocities from UNRES coordinate space to cartesian and Gvec c eigenvector space do i=1,dimen vtnp_(i)=0.0d0 vtnp_a(i)=0.0d0 do j=1,dimen vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j) vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j) enddo vtnp_(i)=vtnp_(i)*dsqrt(geigen(i)) enddo do i=1,dimen write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i) enddo endif endif return end c------------------------------------------------------------------------------- subroutine RESPA_step(itime) c------------------------------------------------------------------------------- c Perform a single RESPA step. c------------------------------------------------------------------------------- implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' double precision energia(0:n_ene),energia_short(0:n_ene), & energia_long(0:n_ene) double precision cm(3),L(3),vcm(3),incr(3) integer ilen,count,rstcount external ilen character*50 tytul integer maxcount_scale /10/ common /gucio/ cm double precision stochforcvec(MAXRES6) double precision grad_tmp(3,0:maxres2) common /stochcalc/ stochforcvec integer itime logical scale double precision vtnp(maxres6), vtnp_(maxres6), vtnp_a(maxres6) if (large.and. mod(itime,ntwe).eq.0) then write (iout,*) "***************** RESPA itime",itime write (iout,*) "Cartesian and internal coordinates: step 0" call cartprint call intout write (iout,*) "Accelerations" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), & (d_a(j,i+nres),j=1,3) enddo write (iout,*) "Velocities, step 0" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), & (d_t(j,i+nres),j=1,3) enddo endif c c Perform the initial RESPA step (increment velocities) c write (iout,*) "*********************** RESPA ini" if (tnp1) then creview call tnp1_respa_step1 call tnp_respa_step1 else if (tnp) then call tnp_respa_step1 else if (tnh.and..not.xiresp) then call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8) do i=0,2*nres do j=1,3 d_t(j,i)=d_t(j,i)*scale_nh enddo enddo endif call RESPA_vel endif cd if(tnp .or. tnp1) then cd write (iout,'(a,3f)') "EE1 NP S, pi",totT, s_np, pi_np cd endif if (mod(itime,ntwe).eq.0 .and. large) then write (iout,*) "Velocities, end" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), & (d_t(j,i+nres),j=1,3) enddo endif c cr if (tnp) then cr do i=0,nres cr do j=1,3 cr grad_tmp(j,i)=gcart(j,i) cr grad_tmp(j,i+nres)=gxcart(j,i) cr enddo cr enddo cr endif c c Compute the short-range forces call zerograd call etotal_short(energia_short) if (tnp.or.tnp1) potE=energia_short(0) call cartgrad call lagrangian do i=0,2*nres do j=1,3 dc_old(j,i)=dc(j,i) if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i) d_a_old(j,i)=d_a(j,i) enddo enddo d_time0=d_time c Split the time step d_time=d_time/ntime_split ctest E_long2=E_long c Perform the short-range RESPSA steps (velocity Verlet increments of c positions and velocities using short-range forces) c write (iout,*) "*********************** RESPA split" creview E_old=potE do itsplit=1,ntime_split if (lang.eq.1) then call sddir_precalc else if (lang.eq.2 .or. lang.eq.3) then call stochastic_force(stochforcvec) endif c First step of the velocity Verlet algorithm if (lang.eq.2) then call sd_verlet1 else if (lang.eq.3) then call sd_verlet1_ciccotti else if (lang.eq.1) then call sddir_verlet1 else if (tnp1) then call tnp1_respa_i_step1 cr if(itsplit.eq.1)then cr d_time_s12=d_time0*0.5*s12_np cr do j=1,3 cr d_t_half(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12 cr enddo cr do i=nnt,nct-1 cr do j=1,3 cr d_t_half(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12 cr enddo cr enddo cr do i=nnt,nct cr if (itype(i).ne.10) then cr inres=i+nres cr do j=1,3 cr d_t_half(j,inres)=d_t_old(j,inres) cr & +d_a_old(j,inres)*d_time_s12 cr enddo cr endif cr enddo cr endif else if (tnp) then call tnp_respa_i_step1 else if (tnh.and.xiresp) then call kinetic(EK) call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8) do i=0,2*nres do j=1,3 d_t_old(j,i)=d_t_old(j,i)*scale_nh enddo enddo cd write(iout,*) "SSS1",itsplit,EK,scale_nh endif call verlet1 endif c Build the chain from the newly calculated coordinates call chainbuild_cart if (rattle) call rattle1 if (large.and. mod(itime,ntwe).eq.0) then write (iout,*) "Cartesian and internal coordinates: step 1" call cartprint call intout write (iout,*) "Accelerations" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), & (d_a(j,i+nres),j=1,3) enddo write (iout,*) "Velocities, step 1" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), & (d_t(j,i+nres),j=1,3) enddo endif tt0 = tcpu() c Calculate energy and forces c c E_long aproximation cr if (tnp .or. tnp1) then cr dtmp=0.5*d_time*(1.0/s12_np+1.0/s_np) cr do i=0,2*nres cr do j=1,3 cr E_long=E_long+d_t_new(j,i)*dtmp*grad_tmp(j,i) cr enddo cr enddo cr endif c------------------------------------- c test of reviewer's comment cr E_long=0 c------------------------------------- c ctest call etotal_long(energia_long) ctest E_long=energia_long(0) ctest call zerograd call etotal_short(energia_short) E_old=potE potE=energia_short(0) c if(tnp .or. tnp1) then c write (iout,*) "kkk",E_long2,E_long c endif call cartgrad c Get the new accelerations call lagrangian t_enegrad=t_enegrad+tcpu()-tt0 c Second step of the velocity Verlet algorithm if (lang.eq.2) then call sd_verlet2 else if (lang.eq.3) then call sd_verlet2_ciccotti else if (lang.eq.1) then call sddir_verlet2 else if (tnp1) then call tnp1_respa_i_step2 else if (tnp) then call tnp_respa_i_step2 else call verlet2 if (tnh.and.xiresp) then call kinetic(EK) call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8) do i=0,2*nres do j=1,3 d_t(j,i)=d_t(j,i)*scale_nh enddo enddo cd write(iout,*) "SSS2",itsplit,EK,scale_nh endif endif if (rattle) call rattle2 c Backup the coordinates, velocities, and accelerations if (tnp .or. tnp1) then do i=0,2*nres do j=1,3 d_t_old(j,i)=d_t(j,i) if (tnp) d_t(j,i)=d_t(j,i)/s_np if (tnp1) d_t(j,i)=d_t(j,i)/s_np enddo enddo endif do i=0,2*nres do j=1,3 dc_old(j,i)=dc(j,i) if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i) d_a_old(j,i)=d_a(j,i) enddo enddo cd if(tnp .or. tnp1) then cd call kinetic(EK) cd HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen) cd H=(HNose1-H0)*s_np cd write (iout,*) "jjj",EK,potE cd write (iout,*) "iii H=",H,abs(HNose1-H0)/H0 cd write (iout,'(a,3f)') cd & "III NP S, pi",totT+itsplit*d_time, s_np, pi_np cd endif enddo c Restore the time step d_time=d_time0 c Compute long-range forces call zerograd call etotal_long(energia_long) E_long=energia_long(0) potE=energia_short(0)+energia_long(0) call cartgrad c Compute accelerations from long-range forces call lagrangian if (large.and. mod(itime,ntwe).eq.0) then write (iout,*) "Cartesian and internal coordinates: step 2" call cartprint call intout write (iout,*) "Accelerations" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), & (d_a(j,i+nres),j=1,3) enddo write (iout,*) "Velocities, step 2" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), & (d_t(j,i+nres),j=1,3) enddo endif c Compute the final RESPA step (increment velocities) c write (iout,*) "*********************** RESPA fin" if (tnp1) then creview call tnp1_respa_step2 call tnp_respa_step2 else if (tnp) then call tnp_respa_step2 else call RESPA_vel if (tnh.and..not.xiresp) then call kinetic(EK) call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8) do i=0,2*nres do j=1,3 d_t(j,i)=d_t(j,i)*scale_nh enddo enddo endif endif if (tnp .or. tnp1) then do i=0,2*nres do j=1,3 d_t(j,i)=d_t_old(j,i)/s_np enddo enddo endif c Compute the complete potential energy cc potE=energia_short(0)+energia_long(0) totT=totT+d_time c Calculate the kinetic and the total energy and the kinetic temperature call kinetic(EK) totE=EK+potE c Couple the system to Berendsen bath if needed if (tbf .and. lang.eq.0) then call verlet_bath endif kinetic_T=2.0d0/(dimen*Rb)*EK c Backup the coordinates, velocities, and accelerations if (mod(itime,ntwe).eq.0 .and. large) then write (iout,*) "Velocities, end" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), & (d_t(j,i+nres),j=1,3) enddo endif c-----review c if(tnp .or. tnp1) then c HNose1=Hnose(EK,s_np,energia_short(0),pi_np,Q_np,t_bath,dimen) c_new_var_csplit Csplit=H0-E_long c Csplit=H0-energia_short(0) c endif c---------- if (mod(itime,ntwe).eq.0) then if(tnp .or. tnp1) then write (iout,'(a3,7f)') "TTT",EK,s_np,potE,pi_np,Csplit, & E_long,energia_short(0) HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen) H=(HNose1-H0)*s_np cd write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0 cd & ,EK+potE+pi_np**2/(2*Q_np)+dimen*0.001986d0*t_bath*log(s_np) write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0 cd write (iout,'(a,3f)') "EE2 NP S, pi",totT, s_np, pi_np endif if(tnh) then HNose1=Hnose_nh(EK,potE) H=HNose1-H0 write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0 endif if (large) then itnp=0 do j=1,3 itnp=itnp+1 vtnp(itnp)=d_t(j,0) enddo do i=nnt,nct-1 do j=1,3 itnp=itnp+1 vtnp(itnp)=d_t(j,i) enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 itnp=itnp+1 vtnp(itnp)=d_t(j,inres) enddo endif enddo c Transform velocities from UNRES coordinate space to cartesian and Gvec c eigenvector space do i=1,dimen vtnp_(i)=0.0d0 vtnp_a(i)=0.0d0 do j=1,dimen vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j) vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j) enddo vtnp_(i)=vtnp_(i)*dsqrt(geigen(i)) enddo do i=1,dimen write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i) enddo endif endif return end c--------------------------------------------------------------------- subroutine RESPA_vel c First and last RESPA step (incrementing velocities using long-range c forces). implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' do j=1,3 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time enddo do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time enddo endif enddo return end c----------------------------------------------------------------- subroutine verlet1 c Applying velocity Verlet algorithm - step 1 to coordinates implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision adt,adt2 do j=1,3 adt=d_a_old(j,0)*d_time adt2=0.5d0*adt dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time d_t_new(j,0)=d_t_old(j,0)+adt2 d_t(j,0)=d_t_old(j,0)+adt enddo do i=nnt,nct-1 do j=1,3 adt=d_a_old(j,i)*d_time adt2=0.5d0*adt dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time d_t_new(j,i)=d_t_old(j,i)+adt2 d_t(j,i)=d_t_old(j,i)+adt enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 adt=d_a_old(j,inres)*d_time adt2=0.5d0*adt dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time d_t_new(j,inres)=d_t_old(j,inres)+adt2 d_t(j,inres)=d_t_old(j,inres)+adt enddo endif enddo return end c--------------------------------------------------------------------- subroutine verlet2 c Step 2 of the velocity Verlet algorithm: update velocities implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' do j=1,3 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time enddo do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time enddo endif enddo return end subroutine sddir_precalc c Applying velocity Verlet algorithm - step 1 to coordinates implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision stochforcvec(MAXRES6) common /stochcalc/ stochforcvec c c Compute friction and stochastic forces c call friction_force call stochastic_force(stochforcvec) c c Compute the acceleration due to friction forces (d_af_work) and stochastic c forces (d_as_work) c do i=1,dimen d_af_work(i)=0.0d0 d_as_work(i)=0.0d0 do j=1,dimen d_af_work(i)=d_af_work(i)+Ginv(i,j)*fric_work(j) d_as_work(i)=d_as_work(i)+Ginv(i,j)*stochforcvec(j) enddo enddo return end c--------------------------------------------------------------------- subroutine sddir_verlet1 c Applying velocity Verlet algorithm - step 1 to velocities implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' c Revised 3/31/05 AL: correlation between random contributions to c position and velocity increments included. double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3) double precision adt,adt2 c c Add the contribution from BOTH friction and stochastic force to the c coordinates, but ONLY the contribution from the friction forces to velocities c do j=1,3 adt=(d_a_old(j,0)+d_af_work(j))*d_time adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt d_t(j,0)=d_t_old(j,0)+adt enddo ind=3 do i=nnt,nct-1 do j=1,3 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt d_t(j,i)=d_t_old(j,i)+adt enddo ind=ind+3 enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt d_t(j,inres)=d_t_old(j,inres)+adt enddo ind=ind+3 endif enddo return end c--------------------------------------------------------------------- subroutine sddir_verlet2 c Calculating the adjusted velocities for accelerations implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6) double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/ c Revised 3/31/05 AL: correlation between random contributions to c position and velocity increments included. c The correlation coefficients are calculated at low-friction limit. c Also, friction forces are now not calculated with new velocities. c call friction_force call stochastic_force(stochforcvec) c c Compute the acceleration due to friction forces (d_af_work) and stochastic c forces (d_as_work) c do i=1,dimen c d_af_work(i)=0.0d0 d_as_work1(i)=0.0d0 do j=1,dimen c d_af_work(i)=d_af_work(i)+Ginv(i,j)*fric_work(j) d_as_work1(i)=d_as_work1(i)+Ginv(i,j)*stochforcvec(j) enddo enddo c c Update velocities c do j=1,3 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time enddo ind=3 do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time enddo ind=ind+3 enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) & +d_af_work(ind+j))+sin60*d_as_work(ind+j) & +cos60*d_as_work1(ind+j))*d_time enddo ind=ind+3 endif enddo return end c--------------------------------------------------------------------- subroutine max_accel c c Find the maximum difference in the accelerations of the the sites c at the beginning and the end of the time step. c implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' double precision aux(3),accel(3) do j=1,3 aux(j)=d_a(j,0)-d_a_old(j,0) enddo amax=0.0d0 do i=nnt,nct c Backbone if (i.lt.nct) then do j=1,3 accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i)) if (dabs(accel(j)).gt.amax) amax=dabs(accel(j)) enddo endif c Side chains do j=1,3 accel(j)=aux(j) enddo if (itype(i).ne.10) then do j=1,3 accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres) enddo endif do j=1,3 if (dabs(accel(j)).gt.amax) amax=dabs(accel(j)) enddo do j=1,3 aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i) enddo enddo return end c--------------------------------------------------------------------- subroutine predict_edrift(epdrift) c c Predict the drift of the potential energy c implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.MUCA' double precision epdrift,epdriftij c Drift of the potential energy epdrift=0.0d0 do i=nnt,nct c Backbone if (i.lt.nct) then do j=1,3 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i)) if (lmuca) epdriftij=epdriftij*factor c write (iout,*) "back",i,j,epdriftij if (epdriftij.gt.epdrift) epdrift=epdriftij enddo endif c Side chains if (itype(i).ne.10) then do j=1,3 epdriftij= & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i)) if (lmuca) epdriftij=epdriftij*factor c write (iout,*) "side",i,j,epdriftij if (epdriftij.gt.epdrift) epdrift=epdriftij enddo endif enddo epdrift=0.5d0*epdrift*d_time*d_time c write (iout,*) "epdrift",epdrift return end c----------------------------------------------------------------------- subroutine verlet_bath c c Coupling to the thermostat by using the Berendsen algorithm c implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision T_half,fact c T_half=2.0d0/(dimen*Rb)*EK fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0)) c write(iout,*) "T_half", T_half c write(iout,*) "EK", EK c write(iout,*) "fact", fact do j=1,3 d_t(j,0)=fact*d_t(j,0) enddo do i=nnt,nct-1 do j=1,3 d_t(j,i)=fact*d_t(j,i) enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t(j,inres)=fact*d_t(j,inres) enddo endif enddo return end c--------------------------------------------------------- subroutine init_MD c Set up the initial conditions of a MD simulation implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MP include 'mpif.h' include 'COMMON.INFO' include 'COMMON.SETUP' character*4 liczba character*16 form #endif include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' real*8 energia(0:n_ene),energia_long(0:n_ene), & energia_short(0:n_ene),vcm(3),incr(3),E_short double precision cm(3),L(3),xv,sigv,lowb,highb character*256 qstr integer ilen external ilen character*50 tytul common /gucio/ cm d_time0=d_time write(iout,*) "d_time", d_time c Compute the standard deviations of stochastic forces for Langevin dynamics c if the friction coefficients do not depend on surface area if (lang.gt.0 .and. .not.surfarea) then do i=nnt,nct-1 stdforcp(i)=stdfp*dsqrt(gamp) enddo do i=nnt,nct stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i))) enddo endif c Open the pdb file for snapshotshots #ifdef MPI if (nprocs.eq.1) then npos=3 else npos = dlog10(dfloat(nprocs-1))+1 endif if (npos.lt.3) npos=3 write (liczba,'(i1)') npos form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) & //')' write (liczba,form) myrank if(mdpdb) then open(ipdb, & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) & //".pdb") else cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) & //".cx" endif #else if(mdpdb) then open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb") else cartname=prefix(:ilen(prefix))//"_MD.cx" endif #endif if (usampl) then write (qstr,'(256(1h ))') ipos=1 do i=1,nfrag iq = qinfrag(i)*10 iw = wfrag(i)/100 if (iw.gt.0) then write (iout,*) "Frag",qinfrag(i),wfrag(i),iq,iw write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw ipos=ipos+7 endif enddo do i=1,npair iq = qinpair(i)*10 iw = wpair(i)/100 if (iw.gt.0) then write (iout,*) "Pair",i,qinpair(i),wpair(i),iq,iw write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw ipos=ipos+7 endif enddo pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb' cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x' statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat' endif icg=1 if (rest) then write(iout,*) "Initial state will be read from file ", & rest2name(:ilen(rest2name)) call readrst else c Generate initial velocities write(iout,*) "Initial velocities randomly generated" call random_vel totT=0.0d0 endif c rest2name = prefix(:ilen(prefix))//'.rst' write(iout,*) "Initial backbone velocities" do i=nnt,nct-1 write(iout,*) (d_t(j,i),j=1,3) enddo write(iout,*) "Initial side-chain velocities" do i=nnt,nct write(iout,*) (d_t(j,i+nres),j=1,3) enddo c Zeroing the total angular momentum of the system write(iout,*) "Calling the zero-angular & momentum subroutine" call inertia_tensor c Getting the potential energy and forces and velocities and accelerations call vcm_vel(vcm) c write (iout,*) "velocity of the center of the mass:" c write (iout,*) (vcm(j),j=1,3) do j=1,3 d_t(j,0)=d_t(j,0)-vcm(j) enddo c Removing the velocity of the center of mass call vcm_vel(vcm) write (iout,*) "vcm right after adjustment:" write (iout,*) (vcm(j),j=1,3) if (.not.rest) then call chainbuild endif call chainbuild_cart call kinetic(EK) if (tbf) then call verlet_bath(EK) endif kinetic_T=2.0d0/(dimen*Rb)*EK call cartprint call intout call zerograd call etotal(energia) potE=energia(0) if(tnp .or. tnp1) then s_np=1.0 pi_np=0.0 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen) H0=Hnose1 write(iout,*) 'H0= ',H0 endif if(tnh) then HNose1=Hnose_nh(EK,potE) H0=HNose1 write (iout,*) 'H0= ',H0 endif call cartgrad call lagrangian call max_accel if (amax*d_time .gt. dvmax) d_time=d_time*dvmax/amax write(iout,*) "Potential energy" write(iout,*) (energia(i),i=0,n_ene) potE=energia(0)-energia(20) totE=EK+potE itime=0 call statout(itime) write (iout,*) "Initial:", & " Kinetic energy",EK," potential energy",potE, & " total energy",totE," maximum acceleration ", & amax if (large) then write (iout,*) "Initial velocities" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), & (d_t(j,i+nres),j=1,3) enddo write (iout,*) "Initial accelerations" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), & (d_a(j,i+nres),j=1,3) enddo endif do i=0,2*nres do j=1,3 dc_old(j,i)=dc(j,i) d_t_old(j,i)=d_t(j,i) d_a_old(j,i)=d_a(j,i) enddo c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3) enddo if (RESPA) then call zerograd call etotal_long(energia_long) E_long=energia_long(0) if(tnp .or. tnp1) then call etotal_short(energia_short) E_short=energia_short(0) HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen) Csplit=Hnose1 c Csplit =110 c_new_var_csplit Csplit=H0-E_long c Csplit = H0-energia_short(0) write(iout,*) 'Csplit= ',Csplit endif call cartgrad call lagrangian c call etotal_short(energia_short) c write (iout,*) "energia_long",energia_long(0), c & " energia_short",energia_short(0), c & " total",energia_long(0)+energia_short(0) endif return end c----------------------------------------------------------- subroutine random_vel implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' double precision xv,sigv,lowb,highb c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m c First generate velocities in the eigenspace of the G matrix c write (iout,*) "Calling random_vel" xv=0.0d0 do i=1,dimen sigv=dsqrt((Rb*t_bath)/geigen(i)) lowb=-5*sigv highb=5*sigv d_t_work_new(i)=anorm_distr(xv,sigv,lowb,highb) enddo c Ek1=0.0d0 c do i=1,dimen c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(i)**2 c enddo c Transform velocities to UNRES coordinate space do i=1,dimen d_t_work(i)=0.0d0 do j=1,dimen d_t_work(i)=d_t_work(i)+Gvec(i,j)*d_t_work_new(j) enddo enddo c Transfer to the d_t vector do j=1,3 d_t(j,0)=d_t_work(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 ind=ind+1 d_t(j,i)=d_t_work(ind) enddo enddo do i=nnt,nct if (itype(i).ne.10) then do j=1,3 ind=ind+1 d_t(j,i+nres)=d_t_work(ind) enddo endif enddo c call kinetic(EK) c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature", c & 2.0d0/(dimen*Rb)*EK,2.0d0/(dimen*Rb)*EK1 return end c----------------------------------------------------------- subroutine sd_verlet_p_setup c Sets up the parameters of stochastic Verlet algorithm implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' double precision emgdt(MAXRES6), & pterm,vterm,rho,rhoc,vsig, & pfric_vec(MAXRES6),vfric_vec(MAXRES6), & afric_vec(MAXRES6),prand_vec(MAXRES6), & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6) logical lprn /.false./ double precision zero /1.0d-8/, gdt_radius /0.05d0/ double precision ktm tt0 = tcpu() c c AL 8/17/04 Code adapted from tinker c c Get the frictional and random terms for stochastic dynamics in the c eigenspace of mass-scaled UNRES friction matrix c do i = 1, dimen gdt = fricgam(i) * d_time c c Stochastic dynamics reduces to simple MD for zero friction c if (gdt .le. zero) then pfric_vec(i) = 1.0d0 vfric_vec(i) = d_time afric_vec(i) = 0.5d0 * d_time * d_time prand_vec(i) = 0.0d0 vrand_vec1(i) = 0.0d0 vrand_vec2(i) = 0.0d0 c c Analytical expressions when friction coefficient is large c else if (gdt .ge. gdt_radius) then egdt = dexp(-gdt) pfric_vec(i) = egdt vfric_vec(i) = (1.0d0-egdt) / fricgam(i) afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i) pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt vterm = 1.0d0 - egdt**2 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm) c c Use series expansions when friction coefficient is small c else gdt2 = gdt * gdt gdt3 = gdt * gdt2 gdt4 = gdt2 * gdt2 gdt5 = gdt2 * gdt3 gdt6 = gdt3 * gdt3 gdt7 = gdt3 * gdt4 gdt8 = gdt4 * gdt4 gdt9 = gdt4 * gdt5 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 & - gdt5/120.0d0 + gdt6/720.0d0 & - gdt7/5040.0d0 + gdt8/40320.0d0 & - gdt9/362880.0d0) / fricgam(i)**2 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i) pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i) pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 & + 127.0d0*gdt9/90720.0d0 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 & - 17.0d0*gdt2/1280.0d0 & + 17.0d0*gdt3/6144.0d0 & + 40967.0d0*gdt4/34406400.0d0 & - 57203.0d0*gdt5/275251200.0d0 & - 1429487.0d0*gdt6/13212057600.0d0) end if c c Compute the scaling factors of random terms for the nonzero friction case c ktm = 0.5d0*d_time/fricgam(i) psig = dsqrt(ktm*pterm) / fricgam(i) vsig = dsqrt(ktm*vterm) rhoc = dsqrt(1.0d0 - rho*rho) prand_vec(i) = psig vrand_vec1(i) = vsig * rho vrand_vec2(i) = vsig * rhoc end if end do if (lprn) then write (iout,*) & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,", & " vrand_vec2" do i=1,dimen write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i), & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i) enddo endif c c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables c call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat) call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat) call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat) call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat) call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec1,vrand_mat1) call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2) c call eigtransf1(dimen,maxres6,mt3mt2,pfric_vec,pfric_mat) c call eigtransf1(dimen,maxres6,mt3mt2,vfric_vec,vfric_mat) c call eigtransf1(dimen,maxres6,mt3mt2,afric_vec,afric_mat) c call eigtransf1(dimen,maxres6,mt3mt1,prand_vec,prand_mat) c call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec1,vrand_mat1) c call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec2,vrand_mat2) t_sdsetup=t_sdsetup+tcpu()-tt0 return end c------------------------------------------------------------- subroutine eigtransf1(n,ndim,ab,d,c) implicit none integer n,ndim double precision ab(ndim,ndim,n),c(ndim,n),d(ndim) integer i,j,k do i=1,n do j=1,n c(i,j)=0.0d0 do k=1,n c(i,j)=c(i,j)+ab(k,j,i)*d(k) enddo enddo enddo return end c------------------------------------------------------------- subroutine eigtransf(n,ndim,a,b,d,c) implicit none integer n,ndim double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim) integer i,j,k do i=1,n do j=1,n c(i,j)=0.0d0 do k=1,n c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k) enddo enddo enddo return end c------------------------------------------------------------- subroutine sd_verlet1 c Applying stochastic velocity Verlet algorithm - step 1 to velocities implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision stochforcvec(MAXRES6) common /stochcalc/ stochforcvec logical lprn /.false./ c write (iout,*) "dc_old" c do i=0,nres c write (iout,'(i5,3f10.5,5x,3f10.5)') c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3) c enddo do j=1,3 dc_work(j)=dc_old(j,0) d_t_work(j)=d_t_old(j,0) d_a_work(j)=d_a_old(j,0) enddo ind=3 do i=nnt,nct-1 do j=1,3 dc_work(ind+j)=dc_old(j,i) d_t_work(ind+j)=d_t_old(j,i) d_a_work(ind+j)=d_a_old(j,i) enddo ind=ind+3 enddo do i=nnt,nct if (itype(i).ne.10) then do j=1,3 dc_work(ind+j)=dc_old(j,i+nres) d_t_work(ind+j)=d_t_old(j,i+nres) d_a_work(ind+j)=d_a_old(j,i+nres) enddo ind=ind+3 endif enddo if (lprn) then write (iout,*) & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,", & " vrand_mat2" do i=1,dimen do j=1,dimen write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j), & vfric_mat(i,j),afric_mat(i,j), & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j) enddo enddo endif do i=1,dimen ddt1=0.0d0 ddt2=0.0d0 do j=1,dimen dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j) ddt1=ddt1+pfric_mat(i,j)*d_t_work(j) ddt2=ddt2+vfric_mat(i,j)*d_a_work(j) enddo d_t_work_new(i)=ddt1+0.5d0*ddt2 d_t_work(i)=ddt1+ddt2 enddo do j=1,3 dc(j,0)=dc_work(j) d_t(j,0)=d_t_work(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 dc(j,i)=dc_work(ind+j) d_t(j,i)=d_t_work(ind+j) enddo ind=ind+3 enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 dc(j,inres)=dc_work(ind+j) d_t(j,inres)=d_t_work(ind+j) enddo ind=ind+3 endif enddo return end c-------------------------------------------------------------------------- subroutine sd_verlet2 c Calculating the adjusted velocities for accelerations implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6) common /stochcalc/ stochforcvec c c Compute the stochastic forces which contribute to velocity change c call stochastic_force(stochforcvecV) do i=1,dimen ddt1=0.0d0 ddt2=0.0d0 do j=1,dimen ddt1=ddt1+vfric_mat(i,j)*d_a_work(j) ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ & vrand_mat2(i,j)*stochforcvecV(j) enddo d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2 enddo do j=1,3 d_t(j,0)=d_t_work(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_work(ind+j) enddo ind=ind+3 enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_work(ind+j) enddo ind=ind+3 endif enddo return end c----------------------------------------------------------- subroutine sd_verlet_ciccotti_setup c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's c version implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' double precision emgdt(MAXRES6), & pterm,vterm,rho,rhoc,vsig, & pfric_vec(MAXRES6),vfric_vec(MAXRES6), & afric_vec(MAXRES6),prand_vec(MAXRES6), & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6) logical lprn /.false./ double precision zero /1.0d-8/, gdt_radius /0.05d0/ double precision ktm tt0 = tcpu() c c AL 8/17/04 Code adapted from tinker c c Get the frictional and random terms for stochastic dynamics in the c eigenspace of mass-scaled UNRES friction matrix c do i = 1, dimen write (iout,*) "i",i," fricgam",fricgam(i) gdt = fricgam(i) * d_time c c Stochastic dynamics reduces to simple MD for zero friction c if (gdt .le. zero) then pfric_vec(i) = 1.0d0 vfric_vec(i) = d_time afric_vec(i) = 0.5d0*d_time*d_time prand_vec(i) = afric_vec(i) vrand_vec2(i) = vfric_vec(i) c c Analytical expressions when friction coefficient is large c else egdt = dexp(-gdt) pfric_vec(i) = egdt vfric_vec(i) = dexp(-0.5d0*gdt)*d_time afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time prand_vec(i) = afric_vec(i) vrand_vec2(i) = vfric_vec(i) c c Compute the scaling factors of random terms for the nonzero friction case c c ktm = 0.5d0*d_time/fricgam(i) c psig = dsqrt(ktm*pterm) / fricgam(i) c vsig = dsqrt(ktm*vterm) c prand_vec(i) = psig*afric_vec(i) c vrand_vec2(i) = vsig*vfric_vec(i) end if end do if (lprn) then write (iout,*) & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,", & " vrand_vec2" do i=1,dimen write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i), & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i) enddo endif c c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables c call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat) call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat) call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat) call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat) call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2) t_sdsetup=t_sdsetup+tcpu()-tt0 return end c------------------------------------------------------------- subroutine sd_verlet1_ciccotti c Applying stochastic velocity Verlet algorithm - step 1 to velocities implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision stochforcvec(MAXRES6) common /stochcalc/ stochforcvec logical lprn /.false./ c write (iout,*) "dc_old" c do i=0,nres c write (iout,'(i5,3f10.5,5x,3f10.5)') c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3) c enddo do j=1,3 dc_work(j)=dc_old(j,0) d_t_work(j)=d_t_old(j,0) d_a_work(j)=d_a_old(j,0) enddo ind=3 do i=nnt,nct-1 do j=1,3 dc_work(ind+j)=dc_old(j,i) d_t_work(ind+j)=d_t_old(j,i) d_a_work(ind+j)=d_a_old(j,i) enddo ind=ind+3 enddo do i=nnt,nct if (itype(i).ne.10) then do j=1,3 dc_work(ind+j)=dc_old(j,i+nres) d_t_work(ind+j)=d_t_old(j,i+nres) d_a_work(ind+j)=d_a_old(j,i+nres) enddo ind=ind+3 endif enddo if (lprn) then write (iout,*) & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,", & " vrand_mat2" do i=1,dimen do j=1,dimen write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j), & vfric_mat(i,j),afric_mat(i,j), & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j) enddo enddo endif do i=1,dimen ddt1=0.0d0 ddt2=0.0d0 do j=1,dimen dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j) ddt1=ddt1+pfric_mat(i,j)*d_t_work(j) ddt2=ddt2+vfric_mat(i,j)*d_a_work(j) enddo d_t_work_new(i)=ddt1+0.5d0*ddt2 d_t_work(i)=ddt1+ddt2 enddo do j=1,3 dc(j,0)=dc_work(j) d_t(j,0)=d_t_work(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 dc(j,i)=dc_work(ind+j) d_t(j,i)=d_t_work(ind+j) enddo ind=ind+3 enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 dc(j,inres)=dc_work(ind+j) d_t(j,inres)=d_t_work(ind+j) enddo ind=ind+3 endif enddo return end c-------------------------------------------------------------------------- subroutine sd_verlet2_ciccotti c Calculating the adjusted velocities for accelerations implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6) common /stochcalc/ stochforcvec c c Compute the stochastic forces which contribute to velocity change c call stochastic_force(stochforcvecV) do i=1,dimen ddt1=0.0d0 ddt2=0.0d0 do j=1,dimen ddt1=ddt1+vfric_mat(i,j)*d_a_work(j) c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j) ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j) enddo d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2 enddo do j=1,3 d_t(j,0)=d_t_work(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_work(ind+j) enddo ind=ind+3 enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_work(ind+j) enddo ind=ind+3 endif enddo return end c------------------------------------------------------ double precision function HNose(ek,s,e,pi,Q,t_bath,dimen) implicit none double precision ek,s,e,pi,Q,t_bath,Rb integer dimen Rb=0.001986d0 HNose=ek+e+pi**2/(2*Q)+dimen*Rb*t_bath*log(s) c print '(6f15.5,i5,a2,2f15.5)',ek,s,e,pi,Q,t_bath,dimen,"--", c & pi**2/(2*Q),dimen*Rb*t_bath*log(s) return end c----------------------------------------------------------------- double precision function HNose_nh(eki,e) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.MD' HNose_nh=eki+e+dimen*Rb*t_bath*xlogs(1)+qmass(1)*vlogs(1)**2/2 do i=2,nnos HNose_nh=HNose_nh+qmass(i)*vlogs(i)**2/2+Rb*t_bath*xlogs(i) enddo c write(4,'(5e15.5)') c & vlogs(1),xlogs(1),HNose,eki,e return end c------------------------------------------------------- subroutine tnp1_step1 c Applying Nose-Poincare algorithm - step 1 to coordinates c JPSJ 70 75 (2001) S. Nose c c d_t is not updated here c implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision adt,adt2,tmp tmp=1+pi_np/(2*Q_np)*0.5*d_time s12_np=s_np*tmp**2 pistar=pi_np/tmp s12_dt=d_time/s12_np d_time_s12=d_time*0.5*s12_np do j=1,3 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt enddo do i=nnt,nct-1 do j=1,3 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt enddo endif enddo return end c--------------------------------------------------------------------- subroutine tnp1_step2 c Step 2 of the velocity Verlet algorithm: update velocities implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision d_time_s12 do i=0,2*nres do j=1,3 d_t(j,i)=d_t_new(j,i) enddo enddo call kinetic(EK) EK=EK/s12_np**2 d_time_s12=0.5d0*s12_np*d_time do j=1,3 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12 enddo do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12 enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12 enddo endif enddo cd write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np pistar=pistar+(EK-0.5*(E_old+potE) & -dimen*Rb*t_bath*log(s12_np)+H0-dimen*Rb*t_bath)*d_time tmp=1+pistar/(2*Q_np)*0.5*d_time s_np=s12_np*tmp**2 pi_np=pistar/tmp return end c--------------------------------------------------------------------- subroutine tnp_step1 c Applying Nose-Poincare algorithm - step 1 to coordinates c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird c c d_t is not updated here, it is destroyed c implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision C_np,d_time_s,tmp,d_time_ss d_time_s=d_time*0.5*s_np do j=1,3 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s enddo do i=nnt,nct-1 do j=1,3 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s enddo endif enddo do i=0,2*nres do j=1,3 d_t(j,i)=d_t_new(j,i) enddo enddo call kinetic(EK) EK=EK/s_np**2 C_np=0.5*d_time*(dimen*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0) & -pi_np pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np)) tmp=0.5*d_time*pistar/Q_np s12_np=s_np*(1.0+tmp)/(1.0-tmp) c write(iout,*) 'tnp_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np) do j=1,3 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss enddo do i=nnt,nct-1 do j=1,3 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss enddo endif enddo return end c----------------------------------------------------------------- subroutine tnp_step2 c Step 2 of the velocity Verlet algorithm: update velocities implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision d_time_s EK=EK*(s_np/s12_np)**2 HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen) pi_np=pistar+0.5*d_time*(2*EK-dimen*Rb*t_bath) & -0.5*d_time*(HNose1-H0) cd write(iout,'(a,4f)') 'mmm',EK,potE,HNose1,pi_np d_time_s=d_time*0.5*s12_np do j=1,3 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s enddo do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s enddo endif enddo s_np=s12_np return end c----------------------------------------------------------------- subroutine tnp1_respa_i_step1 c Applying Nose-Poincare algorithm - step 1 to coordinates c JPSJ 70 75 (2001) S. Nose c c d_t is not updated here c implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision adt,adt2,tmp tmp=1+pi_np/(2*Q_np)*0.5*d_time s12_np=s_np*tmp**2 pistar=pi_np/tmp s12_dt=d_time/s12_np d_time_s12=d_time*0.5*s12_np do j=1,3 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt enddo do i=nnt,nct-1 do j=1,3 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt enddo endif enddo return end c--------------------------------------------------------------------- c----------------------------------------------------------------- subroutine tnp1_respa_step1_ c Applying Nose-Poincare algorithm - step 1 to vel for RESPA c JPSJ 70 75 (2001) S. Nose c c d_t is not updated here c implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision adt,adt2,tmp tmp=1+pi_np/(2*Q_np)*0.5*d_time s12_np=s_np*tmp**2 pistar=pi_np/tmp s12_dt=d_time/s12_np d_time_s12=d_time*0.5*s12_np do j=1,3 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s12 enddo do i=nnt,nct-1 do j=1,3 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s12 enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s12 enddo endif enddo return end c--------------------------------------------------------------------- subroutine tnp1_respa_i_step2 c Step 2 of the velocity Verlet algorithm: update velocities implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision d_time_s12 do i=0,2*nres do j=1,3 d_t(j,i)=d_t_new(j,i) enddo enddo call kinetic(EK) EK=EK/s12_np**2 d_time_s12=0.5d0*s12_np*d_time do j=1,3 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12 enddo do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12 enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12 enddo endif enddo pistar=pistar+(EK-0.5*(E_old+potE) & -dimen*Rb*t_bath*log(s12_np)+Csplit-dimen*Rb*t_bath)*d_time tmp=1+pistar/(2*Q_np)*0.5*d_time s_np=s12_np*tmp**2 pi_np=pistar/tmp return end c--------------------------------------------------------------------- subroutine tnp1_respa_step2_ c Step 2 of the velocity Verlet algorithm: update velocities for RESPA implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision d_time_s12 do i=0,2*nres do j=1,3 d_t(j,i)=d_t_half(j,i) enddo enddo call kinetic(EK) EK=EK/s12_np**2 d_time_s12=0.5d0*s12_np*d_time do j=1,3 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s12 enddo do i=nnt,nct-1 do j=1,3 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s12 enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s12 enddo endif enddo cd write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np pistar=pistar+(EK-0.5*(E_old+potE) & -dimen*Rb*t_bath*log(s12_np)+H0-dimen*Rb*t_bath)*d_time tmp=1+pistar/(2*Q_np)*0.5*d_time s_np=s12_np*tmp**2 pi_np=pistar/tmp return end c----------------------------------------------------------------- subroutine tnp_respa_i_step1 c Applying Nose-Poincare algorithm - step 1 to coordinates c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird c c d_t is not updated here, it is destroyed c implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision C_np,d_time_s,tmp,d_time_ss d_time_s=d_time*0.5*s_np ct2 d_time_s=d_time*0.5*s12_np do j=1,3 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s enddo do i=nnt,nct-1 do j=1,3 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s enddo endif enddo do i=0,2*nres do j=1,3 d_t(j,i)=d_t_new(j,i) enddo enddo call kinetic(EK) EK=EK/s_np**2 C_np=0.5*d_time*(dimen*Rb*t_bath*(1.0+log(s_np))-EK+potE-Csplit) & -pi_np pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np)) tmp=0.5*d_time*pistar/Q_np s12_np=s_np*(1.0+tmp)/(1.0-tmp) d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np) ct2 d_time_ss=d_time/s12_np c d_time_ss=0.5*d_time*(1.0/sold_np+1.0/s_np) do j=1,3 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss enddo do i=nnt,nct-1 do j=1,3 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss enddo endif enddo return end c----------------------------------------------------------------- subroutine tnp_respa_step1 c Applying Nose-Poincare algorithm - step 1 to vel for RESPA c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird c c d_t is not updated here, it is destroyed c implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision C_np,d_time_s,tmp,d_time_ss double precision energia(0:n_ene) d_time_s=d_time*0.5*s_np do j=1,3 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s enddo do i=nnt,nct-1 do j=1,3 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s enddo endif enddo c C_np=0.5*d_time*(dimen*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0) c & -pi_np c c pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np)) c tmp=0.5*d_time*pistar/Q_np c s12_np=s_np*(1.0+tmp)/(1.0-tmp) c write(iout,*) 'tnp_respa_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp ct1 pi_np=pistar c sold_np=s_np c s_np=s12_np c------------------------------------- c test of reviewer's comment pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0) cr print '(a,3f)','1 pi_np,s_np',pi_np,s_np,E_long c------------------------------------- return end c--------------------------------------------------------------------- subroutine tnp_respa_i_step2 c Step 2 of the velocity Verlet algorithm: update velocities implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision d_time_s EK=EK*(s_np/s12_np)**2 HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen) pi_np=pistar+0.5*d_time*(2*EK-dimen*Rb*t_bath & -HNose1+Csplit) cr print '(a,5f)','i_step2',EK,potE,HNose1,pi_np,E_long d_time_s=d_time*0.5*s12_np c d_time_s=d_time*0.5*s_np do j=1,3 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s enddo do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s enddo endif enddo s_np=s12_np return end c--------------------------------------------------------------------- subroutine tnp_respa_step2 c Step 2 of the velocity Verlet algorithm: update velocities for RESPA implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision d_time_s ct1 s12_np=s_np ct2 pistar=pi_np ct call kinetic(EK) ct HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen) ct pi_np=pistar+0.5*d_time*(2*EK-dimen*Rb*t_bath) ct & -0.5*d_time*(HNose1-H0) c------------------------------------- c test of reviewer's comment pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0) cr print '(a,3f)','2 pi_np,s_np',pi_np,s_np,E_long c------------------------------------- d_time_s=d_time*0.5*s_np do j=1,3 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s enddo do i=nnt,nct-1 do j=1,3 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s enddo enddo do i=nnt,nct if (itype(i).ne.10) then inres=i+nres do j=1,3 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s enddo endif enddo cd s_np=s12_np return end c----------------------------------------------------------------- SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.MD' double precision akin,gnkt,dt,aa,gkt,scale double precision wdti(maxyosh),wdti2(maxyosh), & wdti4(maxyosh),wdti8(maxyosh) integer i,iresn,iyosh,inos,nnos1 dt=d_time nnos1=nnos+1 GKT = Rb*t_bath GNKT = dimen*GKT akin=akin*2 C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE C INTEGRATION FROM t=0 TO t=DT/2 C GET THE TOTAL KINETIC ENERGY SCALE = 1.D0 c CALL GETKINP(MASS,VX,VY,VZ,AKIN) C UPDATE THE FORCES GLOGS(1) = (AKIN - GNKT)/QMASS(1) C START THE MULTIPLE TIME STEP PROCEDURE DO IRESN = 1,NRESN DO IYOSH = 1,NYOSH C UPDATE THE THERMOSTAT VELOCITIES VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH) DO INOS = 1,NNOS-1 AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) ) VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA & + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA ENDDO C UPDATE THE PARTICLE VELOCITIES AA = EXP(-WDTI2(IYOSH)*VLOGS(1) ) SCALE = SCALE*AA C UPDATE THE FORCES GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1) C UPDATE THE THERMOSTAT POSITIONS DO INOS = 1,NNOS XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH) ENDDO C UPDATE THE THERMOSTAT VELOCITIES DO INOS = 1,NNOS-1 AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) ) VLOGS(INOS) = VLOGS(INOS)*AA*AA & + WDTI4(IYOSH)*GLOGS(INOS)*AA GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS) & -GKT)/QMASS(INOS+1) ENDDO VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH) ENDDO ENDDO C UPDATE THE PARTICLE VELOCITIES c outside of this subroutine c DO I = 1,N c VX(I) = VX(I)*SCALE c VY(I) = VY(I)*SCALE c VZ(I) = VZ(I)*SCALE c ENDDO RETURN END