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 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 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 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 call verlet2 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 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) d_t_old(j,i)=d_t(j,i) d_a_old(j,i)=d_a(j,i) enddo enddo if (mod(itime,ntwe).eq.0 .and. large) then 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 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) common /stochcalc/ stochforcvec integer itime logical scale 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" call RESPA_vel 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 Compute the short-range forces call zerograd call etotal_short(energia_short) call cartgrad call lagrangian 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 enddo d_time0=d_time c Split the time step d_time=d_time/ntime_split c Perform the short-range RESPSA steps (velocity Verlet increments of c positions and velocities using short-range forces) c write (iout,*) "*********************** RESPA split" 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 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 call etotal_short(energia_short) 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 call verlet2 endif if (rattle) call rattle2 c Backup the coordinates, velocities, and accelerations 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 enddo enddo c Restore the time step d_time=d_time0 c Compute long-range forces call zerograd call etotal_long(energia_long) 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" call RESPA_vel c Compute the complete potential energy 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 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 c----------------------------------------------------------------- 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) 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) 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 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) 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