X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fmd-diff%2FMD_NP.F;fp=source%2Funres%2Fsrc_MD%2Fmd-diff%2FMD_NP.F;h=0000000000000000000000000000000000000000;hb=0a11a2c4ccee14ed99ae44f2565b270ba8d4bbb6;hp=06e89ea18b5e019060f7a8e867f4b9ea3cbce73b;hpb=5eb407964903815242c59de10960f42761139e10;p=unres.git diff --git a/source/unres/src_MD/md-diff/MD_NP.F b/source/unres/src_MD/md-diff/MD_NP.F deleted file mode 100644 index 06e89ea..0000000 --- a/source/unres/src_MD/md-diff/MD_NP.F +++ /dev/null @@ -1,3290 +0,0 @@ - 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