X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fmd-diff%2FMD_A-MTS.F;fp=source%2Funres%2Fsrc_MD%2Fmd-diff%2FMD_A-MTS.F;h=0000000000000000000000000000000000000000;hb=0a11a2c4ccee14ed99ae44f2565b270ba8d4bbb6;hp=6b523793ccf196be91e1e52e816e6341d7b06467;hpb=5eb407964903815242c59de10960f42761139e10;p=unres.git diff --git a/source/unres/src_MD/md-diff/MD_A-MTS.F b/source/unres/src_MD/md-diff/MD_A-MTS.F deleted file mode 100644 index 6b52379..0000000 --- a/source/unres/src_MD/md-diff/MD_A-MTS.F +++ /dev/null @@ -1,2423 +0,0 @@ - subroutine MD -c------------------------------------------------ -c The driver for molecular dynamics subroutines -c------------------------------------------------ - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' -#ifdef MPI - include "mpif.h" - integer IERROR,ERRCODE -#endif - include 'COMMON.SETUP' - include 'COMMON.CONTROL' - include 'COMMON.VAR' - include 'COMMON.MD' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - 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) -#ifdef VOUT - double precision v_work(maxres6),v_transf(maxres6) -#endif - integer ilen,rstcount - external ilen - character*50 tytul - common /gucio/ cm - integer itime -c -#ifdef MPI - if (ilen(tmpdir).gt.0) - & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" - & //liczba(:ilen(liczba))//'.rst') -#else - if (ilen(tmpdir).gt.0) - & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst') -#endif - 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" -#ifdef MPI - tt0=MPI_Wtime() -#else - tt0 = tcpu() -#endif -c Determine the inverse of the inertia matrix. - call setup_MD_matrices -c Initialize MD - call init_MD -#ifdef MPI - t_MDsetup = MPI_Wtime()-tt0 -#else - t_MDsetup = tcpu()-tt0 -#endif - rstcount=0 -c Entering the MD loop -#ifdef MPI - tt0 = MPI_Wtime() -#else - tt0 = tcpu() -#endif - if (lang.eq.2 .or. lang.eq.3) then -#ifndef LANG0 - 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 - write (iout,*) - & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0" -#ifdef MPI - call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) -#endif - stop -#endif - else if (lang.eq.1 .or. lang.eq.4) then - call setup_fricmat - endif -#ifdef MPI - t_langsetup=MPI_Wtime()-tt0 - tt0=MPI_Wtime() -#else - t_langsetup=tcpu()-tt0 - tt0=tcpu() -#endif - 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 -#ifndef LANG0 - 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 -#endif - 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/(dimen3*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 -#ifdef BROWN - call brown_step(itime) -#else - print *,"Brown dynamics not here!" -#ifdef MPI - call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) -#endif - stop -#endif - endif - if (ntwe.ne.0) then - if (mod(itime,ntwe).eq.0) call statout(itime) -#ifdef VOUT - do j=1,3 - v_work(j)=d_t(j,0) - enddo - ind=3 - do i=nnt,nct-1 - do j=1,3 - ind=ind+1 - v_work(ind)=d_t(j,i) - enddo - enddo - do i=nnt,nct - if (itype(i).ne.10) then - do j=1,3 - ind=ind+1 - v_work(ind)=d_t(j,i+nres) - enddo - endif - enddo - - write (66,'(80f10.5)') - & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres) - do i=1,ind - v_transf(i)=0.0d0 - do j=1,ind - v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j) - enddo - v_transf(i)= v_transf(i)*dsqrt(geigen(i)) - enddo - write (67,'(80f10.5)') (v_transf(i),i=1,ind) -#endif - endif - 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,t_bath - 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 -#ifdef MPI - t_MD=MPI_Wtime()-tt0 -#else - t_MD=tcpu()-tt0 -#endif - 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 ' -#ifdef TIMING_ENE - write (iout,*) "time for etotal",t_etotal," elong",t_elong, - & " eshort",t_eshort - write (iout,*) "time_fric",time_fric," time_stoch",time_stoch, - & " time_fricmatmult",time_fricmatmult," time_fsample ", - & time_fsample -#endif - 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' -#ifdef MPI - include 'mpif.h' - integer ierror,ierrcode -#endif - include 'COMMON.SETUP' - include 'COMMON.CONTROL' - include 'COMMON.VAR' - include 'COMMON.MD' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - 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 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 - 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 -#ifndef LANG0 - call stochastic_force(stochforcvec) -#else - write (iout,*) - & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0" -#ifdef MPI - call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) -#endif - stop -#endif - 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 - call flush(iout) -#ifdef MPI - call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE) -#endif - stop - endif -c First step of the velocity Verlet algorithm - if (lang.eq.2) then -#ifndef LANG0 - call sd_verlet1 -#endif - else if (lang.eq.3) then -#ifndef LANG0 - call sd_verlet1_ciccotti -#endif - 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 (ntwe.ne.0) then - if (large.and. mod(itime,ntwe).eq.0) then - write (iout,*) "Cartesian and internal coordinates: step 1" - call cartprint - call intout - write (iout,*) "dC" - do i=0,nres - write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3), - & (dc(j,i+nres),j=1,3) - enddo - 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 - endif -#ifdef MPI - tt0 = MPI_Wtime() -#else - tt0 = tcpu() -#endif -c Calculate energy and forces - call zerograd - call etotal(potEcomp) -#ifdef TIMING_ENE -#ifdef MPI - t_etotal=t_etotal+MPI_Wtime()-tt0 -#else - t_etotal=t_etotal+tcpu()-tt0 -#endif -#endif - potE=potEcomp(0)-potEcomp(20) - call cartgrad -c Get the new accelerations - call lagrangian -#ifdef MPI - t_enegrad=t_enegrad+MPI_Wtime()-tt0 -#else - t_enegrad=t_enegrad+tcpu()-tt0 -#endif -c Determine maximum acceleration and scale down the timestep if needed - call max_accel - amax=amax/(itime_scal+1)**2 - call predict_edrift(epdrift) - if (amax/(itime_scal+1).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 -#ifndef LANG0 -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 -#endif - 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 -#ifndef LANG0 - call sd_verlet2 -#endif - else if (lang.eq.3) then -#ifndef LANG0 - call sd_verlet2_ciccotti -#endif - 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 -#ifndef LANG0 - 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 - 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/(dimen3*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 (ntwe.ne.0) then - 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 - 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' -#ifdef MPI - include 'mpif.h' - integer IERROR,ERRCODE -#endif - include 'COMMON.SETUP' - include 'COMMON.CONTROL' - include 'COMMON.VAR' - include 'COMMON.MD' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - 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_short(0:n_ene), - & energia_long(0:n_ene) - double precision cm(3),L(3),vcm(3),incr(3) - double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2), - & d_a_old0(3,0:maxres2) - 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 - common /cipiszcze/ itt - itt=itime - if (ntwe.ne.0) then - if (large.and. mod(itime,ntwe).eq.0) then - write (iout,*) "***************** RESPA itime",itime - write (iout,*) "Cartesian and internal coordinates: step 0" -c call cartprint - call pdbout(0.0d0,"cipiszcze",iout) - call intout - write (iout,*) "Accelerations from long-range forces" - 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 - endif -c -c Perform the initial RESPA step (increment velocities) -c write (iout,*) "*********************** RESPA ini" - call RESPA_vel - if (ntwe.ne.0) then - 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 - endif -c Compute the short-range forces -#ifdef MPI - tt0 =MPI_Wtime() -#else - tt0 = tcpu() -#endif -C 7/2/2009 commented out -c call zerograd -c call etotal_short(energia_short) -c call cartgrad -c call lagrangian -C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step - do i=0,2*nres - do j=1,3 - d_a(j,i)=d_a_short(j,i) - enddo - enddo - if (ntwe.ne.0) then - if (large.and. mod(itime,ntwe).eq.0) then - write (iout,*) "energia_short",energia_short(0) - write (iout,*) "Accelerations from short-range forces" - 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 - endif -#ifdef MPI - t_enegrad=t_enegrad+MPI_Wtime()-tt0 -#else - t_enegrad=t_enegrad+tcpu()-tt0 -#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 - enddo -c 6/30/08 A-MTS: attempt at increasing the split number - do i=0,2*nres - do j=1,3 - dc_old0(j,i)=dc_old(j,i) - d_t_old0(j,i)=d_t_old(j,i) - d_a_old0(j,i)=d_a_old(j,i) - enddo - enddo - if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2 - if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0 -c - scale=.true. - d_time0=d_time - do while (scale) - - scale=.false. -c write (iout,*) "itime",itime," ntime_split",ntime_split -c Split the time step - d_time=d_time0/ntime_split -c Perform the short-range RESPA 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 -#ifndef LANG0 - call stochastic_force(stochforcvec) -#else - write (iout,*) - & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0" -#ifdef MPI - call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) -#endif - stop -#endif - endif -c First step of the velocity Verlet algorithm - if (lang.eq.2) then -#ifndef LANG0 - call sd_verlet1 -#endif - else if (lang.eq.3) then -#ifndef LANG0 - call sd_verlet1_ciccotti -#endif - 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 (ntwe.ne.0) then - if (large.and. mod(itime,ntwe).eq.0) then - write (iout,*) "***** ITSPLIT",itsplit - write (iout,*) "Cartesian and internal coordinates: step 1" - call pdbout(0.0d0,"cipiszcze",iout) -c call cartprint - call intout - 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 - endif -#ifdef MPI - tt0 = MPI_Wtime() -#else - tt0 = tcpu() -#endif -c Calculate energy and forces - call zerograd - call etotal_short(energia_short) -#ifdef TIMING_ENE -#ifdef MPI - t_eshort=t_eshort+MPI_Wtime()-tt0 -#else - t_eshort=t_eshort+tcpu()-tt0 -#endif -#endif - call cartgrad -c Get the new accelerations - call lagrangian -C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array - do i=0,2*nres - do j=1,3 - d_a_short(j,i)=d_a(j,i) - enddo - enddo - if (ntwe.ne.0) then - if (large.and. mod(itime,ntwe).eq.0) then - write (iout,*)"energia_short",energia_short(0) - write (iout,*) "Accelerations from short-range forces" - 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 - endif -c 6/30/08 A-MTS -c Determine maximum acceleration and scale down the timestep if needed - call max_accel - amax=amax/ntime_split**2 - call predict_edrift(epdrift) - if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) - & write (iout,*) "amax",amax," damax",damax, - & " epdrift",epdrift," epdriftmax",epdriftmax -c Exit loop and try with increased split number if the change of -c acceleration is too big - if (amax.gt.damax .or. epdrift.gt.edriftmax) then - if (ntime_split.lt.maxtime_split) then - scale=.true. - ntime_split=ntime_split*2 - do i=0,2*nres - do j=1,3 - dc_old(j,i)=dc_old0(j,i) - d_t_old(j,i)=d_t_old0(j,i) - d_a_old(j,i)=d_a_old0(j,i) - enddo - enddo - write (iout,*) "acceleration/energy drift too large",amax, - & epdrift," split increased to ",ntime_split," itime",itime, - & " itsplit",itsplit - exit - else - write (iout,*) - & "Uh-hu. Bumpy landscape. Maximum splitting number", - & maxtime_split, - & " already reached!!! Trying to carry on!" - endif - endif -#ifdef MPI - t_enegrad=t_enegrad+MPI_Wtime()-tt0 -#else - t_enegrad=t_enegrad+tcpu()-tt0 -#endif -c Second step of the velocity Verlet algorithm - if (lang.eq.2) then -#ifndef LANG0 - call sd_verlet2 -#endif - else if (lang.eq.3) then -#ifndef LANG0 - call sd_verlet2_ciccotti -#endif - 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 - - enddo ! while scale - -c Restore the time step - d_time=d_time0 -c Compute long-range forces -#ifdef MPI - tt0 =MPI_Wtime() -#else - tt0 = tcpu() -#endif - call zerograd - call etotal_long(energia_long) -#ifdef TIMING_ENE -#ifdef MPI - t_elong=t_elong+MPI_Wtime()-tt0 -#else - t_elong=t_elong+tcpu()-tt0 -#endif -#endif - call cartgrad - call lagrangian -#ifdef MPI - t_enegrad=t_enegrad+MPI_Wtime()-tt0 -#else - t_enegrad=t_enegrad+tcpu()-tt0 -#endif -c Compute accelerations from long-range forces - if (ntwe.ne.0) then - if (large.and. mod(itime,ntwe).eq.0) then - write (iout,*) "energia_long",energia_long(0) - write (iout,*) "Cartesian and internal coordinates: step 2" -c call cartprint - call pdbout(0.0d0,"cipiszcze",iout) - call intout - write (iout,*) "Accelerations from long-range forces" - 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 - endif -c Compute the final RESPA step (increment velocities) -c write (iout,*) "*********************** RESPA fin" - call RESPA_vel -c Compute the complete potential energy - do i=0,n_ene - potEcomp(i)=energia_short(i)+energia_long(i) - enddo - potE=potEcomp(0)-potEcomp(20) -c 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/(dimen3*Rb)*EK -c Backup the coordinates, velocities, and accelerations - if (ntwe.ne.0) then - 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 - 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 - -#ifdef DEBUG - write (iout,*) "VELVERLET1 START: DC" - do i=0,nres - write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3), - & (dc(j,i+nres),j=1,3) - enddo -#endif - 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 -#ifdef DEBUG - write (iout,*) "VELVERLET1 END: DC" - do i=0,nres - write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3), - & (dc(j,i+nres),j=1,3) - enddo -#endif - 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' -#ifdef MPI - include 'mpif.h' -#endif - include 'COMMON.CONTROL' - include 'COMMON.VAR' - include 'COMMON.MD' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - 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 stochforcvec(MAXRES6) - common /stochcalc/ stochforcvec -c -c Compute friction and stochastic forces -c - time00=MPI_Wtime() - call friction_force - time_fric=time_fric+MPI_Wtime()-time00 - time00=MPI_Wtime() - call stochastic_force(stochforcvec) - time_stoch=time_stoch+MPI_Wtime()-time00 -c -c Compute the acceleration due to friction forces (d_af_work) and stochastic -c forces (d_as_work) -c - call ginv_mult(fric_work, d_af_work) - call ginv_mult(stochforcvec, d_as_work) - 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' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - 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' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - 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 - call ginv_mult(stochforcvec, d_as_work1) - -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),accel_old(3),dacc - do j=1,3 -c aux(j)=d_a(j,0)-d_a_old(j,0) - accel_old(j)=d_a_old(j,0) - accel(j)=d_a(j,0) - enddo - amax=0.0d0 - do i=nnt,nct -c Backbone - if (i.lt.nct) then -c 7/3/08 changed to asymmetric difference - do j=1,3 -c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i)) - accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i) - accel(j)=accel(j)+0.5d0*d_a(j,i) -c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j)) - if (dabs(accel(j)).gt.dabs(accel_old(j))) then - dacc=dabs(accel(j)-accel_old(j)) - if (dacc.gt.amax) amax=dacc - endif - enddo - endif - enddo -c Side chains - do j=1,3 -c accel(j)=aux(j) - accel_old(j)=d_a_old(j,0) - accel(j)=d_a(j,0) - enddo - if (nnt.eq.2) then - do j=1,3 - accel_old(j)=accel_old(j)+d_a_old(j,1) - accel(j)=accel(j)+d_a(j,1) - enddo - endif - do i=nnt,nct - if (itype(i).ne.10) then - do j=1,3 -c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres) - accel_old(j)=accel_old(j)+d_a_old(j,i+nres) - accel(j)=accel(j)+d_a(j,i+nres) - enddo - endif - do j=1,3 -c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j)) - if (dabs(accel(j)).gt.dabs(accel_old(j))) then - dacc=dabs(accel(j)-accel_old(j)) - if (dacc.gt.amax) amax=dacc - endif - enddo - do j=1,3 - accel_old(j)=accel_old(j)+d_a_old(j,i) - accel(j)=accel(j)+d_a(j,i) -c 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/(dimen3*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' - character*16 form - integer IERROR,ERRCODE -#endif - include 'COMMON.SETUP' - include 'COMMON.CONTROL' - include 'COMMON.VAR' - include 'COMMON.MD' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - include 'COMMON.CHAIN' - include 'COMMON.DERIV' - include 'COMMON.GEO' - include 'COMMON.LOCAL' - include 'COMMON.INTERACT' - include 'COMMON.IOUNITS' - include 'COMMON.NAMES' - include 'COMMON.REMD' - real*8 energia_long(0:n_ene), - & energia_short(0:n_ene),vcm(3),incr(3) - double precision cm(3),L(3),xv,sigv,lowb,highb - double precision varia(maxvar) - character*256 qstr - integer ilen - external ilen - character*50 tytul - logical file_exist - common /gucio/ cm - d_time0=d_time -c 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(mdpdb) then - if (ilen(tmpdir).gt.0) - & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// - & liczba(:ilen(liczba))//".pdb") - open(ipdb, - & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) - & //".pdb") - else -#ifdef NOXDR - if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) - & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// - & liczba(:ilen(liczba))//".x") - cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) - & //".x" -#else - if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) - & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// - & liczba(:ilen(liczba))//".cx") - cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) - & //".cx" -#endif - endif -#else - if(mdpdb) then - if (ilen(tmpdir).gt.0) - & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb") - open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb") - else - if (ilen(tmpdir).gt.0) - & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx") - cartname=prefix(:ilen(prefix))//"_MD.cx" - endif -#endif - if (usampl) then - write (qstr,'(256(1h ))') - ipos=1 - do i=1,nfrag - iq = qinfrag(i,iset)*10 - iw = wfrag(i,iset)/100 - if (iw.gt.0) then - if(me.eq.king.or..not.out1file) - & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),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,iset)*10 - iw = wpair(i,iset)/100 - if (iw.gt.0) then - if(me.eq.king.or..not.out1file) - & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw - write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw - ipos=ipos+7 - endif - enddo -c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb' -#ifdef NOXDR -c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x' -#else -c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx' -#endif -c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat' - endif - icg=1 - if (rest) then - if (restart1file) then - if (me.eq.king) - & inquire(file=mremd_rst_name,exist=file_exist) - write (*,*) me," Before broadcast: file_exist",file_exist - call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM, - & IERR) - write (*,*) me," After broadcast: file_exist",file_exist -c inquire(file=mremd_rst_name,exist=file_exist) - if(me.eq.king.or..not.out1file) - & write(iout,*) "Initial state read by master and distributed" - else - if (ilen(tmpdir).gt.0) - & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' - & //liczba(:ilen(liczba))//'.rst') - inquire(file=rest2name,exist=file_exist) - endif - if(file_exist) then - if(.not.restart1file) then - if(me.eq.king.or..not.out1file) - & write(iout,*) "Initial state will be read from file ", - & rest2name(:ilen(rest2name)) - call readrst - endif - call rescale_weights(t_bath) - else - if(me.eq.king.or..not.out1file)then - if (restart1file) then - write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)), - & " does not exist" - else - write(iout,*) "File ",rest2name(:ilen(rest2name)), - & " does not exist" - endif - write(iout,*) "Initial velocities randomly generated" - endif - call random_vel - totT=0.0d0 - endif - else -c Generate initial velocities - if(me.eq.king.or..not.out1file) - & write(iout,*) "Initial velocities randomly generated" - call random_vel - totT=0.0d0 - endif -c rest2name = prefix(:ilen(prefix))//'.rst' - if(me.eq.king.or..not.out1file)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 -c Zeroing the total angular momentum of the system - write(iout,*) "Calling the zero-angular - & momentum subroutine" - endif - 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) - if(me.eq.king.or..not.out1file)then - write (iout,*) "vcm right after adjustment:" - write (iout,*) (vcm(j),j=1,3) - endif - if (.not.rest) then - call chainbuild - if(iranconf.ne.0) then - if (overlapsc) then - print *, 'Calling OVERLAP_SC' - call overlap_sc(fail) - endif - - if (searchsc) then - call sc_move(2,nres-1,10,1d10,nft_sc,etot) - print *,'SC_move',nft_sc,etot - if(me.eq.king.or..not.out1file) - & write(iout,*) 'SC_move',nft_sc,etot - endif - - if(dccart)then - print *, 'Calling MINIM_DC' - call minim_dc(etot,iretcode,nfun) - else - call geom_to_var(nvar,varia) - print *,'Calling MINIMIZE.' - call minimize(etot,varia,iretcode,nfun) - call var_to_geom(nvar,varia) - endif - if(me.eq.king.or..not.out1file) - & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun - endif - endif - call chainbuild_cart - call kinetic(EK) - if (tbf) then - call verlet_bath(EK) - endif - kinetic_T=2.0d0/(dimen3*Rb)*EK - if(me.eq.king.or..not.out1file)then - call cartprint - call intout - endif -#ifdef MPI - tt0=MPI_Wtime() -#else - tt0=tcpu() -#endif - call zerograd - call etotal(potEcomp) -#ifdef TIMING_ENE -#ifdef MPI - t_etotal=t_etotal+MPI_Wtime()-tt0 -#else - t_etotal=t_etotal+tcpu()-tt0 -#endif -#endif - potE=potEcomp(0) - call cartgrad - call lagrangian - call max_accel - if (amax*d_time .gt. dvmax) then - d_time=d_time*dvmax/amax - if(me.eq.king.or..not.out1file) write (iout,*) - & "Time step reduced to",d_time, - & " because of too large initial acceleration." - endif - if(me.eq.king.or..not.out1file)then - write(iout,*) "Potential energy and its components" - call enerprint(potEcomp) -c write(iout,*) (potEcomp(i),i=0,n_ene) - endif - potE=potEcomp(0)-potEcomp(20) - totE=EK+potE - itime=0 - if (ntwe.ne.0) call statout(itime) - if(me.eq.king.or..not.out1file) - & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", - & " Kinetic energy",EK," potential energy",potE, - & " total energy",totE," maximum acceleration ", - & amax - if (large) then - write (iout,*) "Initial coordinates" - do i=1,nres - write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3), - & (c(j,i+nres),j=1,3) - enddo - write (iout,*) "Initial dC" - do i=0,nres - write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3), - & (dc(j,i+nres),j=1,3) - enddo - write (iout,*) "Initial velocities" - write (iout,"(13x,' backbone ',23x,' side chain')") - 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 -c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), - write (iout,'(i3,3f15.10,3x,3f15.10)') 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 -#ifdef MPI - tt0 =MPI_Wtime() -#else - tt0 = tcpu() -#endif - call zerograd - call etotal_short(energia_short) -#ifdef TIMING_ENE -#ifdef MPI - t_eshort=t_eshort+MPI_Wtime()-tt0 -#else - t_eshort=t_eshort+tcpu()-tt0 -#endif -#endif - call cartgrad - call lagrangian - if(.not.out1file .and. large) then - write (iout,*) "energia_long",energia_long(0), - & " energia_short",energia_short(0), - & " total",energia_long(0)+energia_short(0) - write (iout,*) "Initial fast-force 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 -C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array - do i=0,2*nres - do j=1,3 - d_a_short(j,i)=d_a(j,i) - enddo - enddo -#ifdef MPI - tt0=MPI_Wtime() -#else - tt0=tcpu() -#endif - call zerograd - call etotal_long(energia_long) -#ifdef TIMING_ENE -#ifdef MPI - t_elong=t_elong+MPI_Wtime()-tt0 -#else - t_elong=t_elong+tcpu()-tt0 -#endif -#endif - call cartgrad - call lagrangian - if(.not.out1file .and. large) then - write (iout,*) "energia_long",energia_long(0) - write (iout,*) "Initial slow-force 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 -#ifdef MPI - t_enegrad=t_enegrad+MPI_Wtime()-tt0 -#else - t_enegrad=t_enegrad+tcpu()-tt0 -#endif - 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' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - 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 dimen dimen3",dimen,dimen3 -c call flush(iout) - xv=0.0d0 - ii=0 - do i=1,dimen - do k=1,3 - ii=ii+1 - sigv=dsqrt((Rb*t_bath)/geigen(i)) - lowb=-5*sigv - highb=5*sigv - d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb) -c write (iout,*) "i",i," ii",ii," geigen",geigen(i), -c & " d_t_work_new",d_t_work_new(ii) - enddo - enddo -c diagnostics -c Ek1=0.0d0 -c ii=0 -c do i=1,dimen -c do k=1,3 -c ii=ii+1 -c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2 -c enddo -c enddo -c write (iout,*) "Ek from eigenvectors",Ek1 -c end diagnostics -c Transform velocities to UNRES coordinate space - do k=0,2 - do i=1,dimen - ind=(i-1)*3+k+1 - d_t_work(ind)=0.0d0 - do j=1,dimen - d_t_work(ind)=d_t_work(ind) - & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1) - enddo -c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind) -c call flush(iout) - 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/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1 -c call flush(iout) - return - end -#ifndef LANG0 -c----------------------------------------------------------- - subroutine sd_verlet_p_setup -c Sets up the parameters of stochastic Verlet algorithm - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' -#ifdef MPI - include 'mpif.h' -#endif - include 'COMMON.CONTROL' - include 'COMMON.VAR' - include 'COMMON.MD' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - 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 -#ifdef MPI - tt0 = MPI_Wtime() -#else - tt0 = tcpu() -#endif -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 -#ifndef LANG0 - call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat) - call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat) - call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat) - call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat) - call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1) - call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2) -#endif -#ifdef MPI - t_sdsetup=t_sdsetup+MPI_Wtime() -#else - t_sdsetup=t_sdsetup+tcpu()-tt0 -#endif - 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' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - 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 -#ifndef LANG0 - 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 -#endif - 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' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - 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) - -#ifndef LANG0 - 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 -#endif - 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' -#ifdef MPI - include 'mpif.h' -#endif - include 'COMMON.CONTROL' - include 'COMMON.VAR' - include 'COMMON.MD' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - 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 -#ifdef MPI - tt0 = MPI_Wtime() -#else - tt0 = tcpu() -#endif -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,maxres2,mt3,mt2,pfric_vec,pfric_mat) - call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat) - call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat) - call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat) - call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2) -#ifdef MPI - t_sdsetup=t_sdsetup+MPI_Wtime() -#else - t_sdsetup=t_sdsetup+tcpu()-tt0 -#endif - 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' -#ifdef MPI - include 'mpif.h' -#endif - include 'COMMON.CONTROL' - include 'COMMON.VAR' - include 'COMMON.MD' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - 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 - -#ifndef LANG0 - 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 -#endif - 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' -#ifndef LANG0 - include 'COMMON.LANGEVIN' -#else - include 'COMMON.LANGEVIN.lang0' -#endif - 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) -#ifndef LANG0 - 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 -#endif - 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 -#endif