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