2 c------------------------------------------------
3 c The driver for molecular dynamics subroutines
4 c------------------------------------------------
5 implicit real*8 (a-h,o-z)
11 include 'COMMON.SETUP'
12 include 'COMMON.CONTROL'
16 include 'COMMON.LANGEVIN'
18 include 'COMMON.LANGEVIN.lang0'
20 include 'COMMON.CHAIN'
21 include 'COMMON.DERIV'
23 include 'COMMON.LOCAL'
24 include 'COMMON.INTERACT'
25 include 'COMMON.IOUNITS'
26 include 'COMMON.NAMES'
27 include 'COMMON.TIME1'
28 double precision cm(3),L(3),vcm(3)
30 double precision v_work(maxres6),v_transf(maxres6)
39 if (ilen(tmpdir).gt.0)
40 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
41 & //liczba(:ilen(liczba))//'.rst')
43 if (ilen(tmpdir).gt.0)
44 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
51 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
57 c Determine the inverse of the inertia matrix.
58 call setup_MD_matrices
62 t_MDsetup = MPI_Wtime()-tt0
64 t_MDsetup = tcpu()-tt0
67 c Entering the MD loop
73 if (lang.eq.2 .or. lang.eq.3) then
77 call sd_verlet_p_setup
79 call sd_verlet_ciccotti_setup
83 pfric0_mat(i,j,0)=pfric_mat(i,j)
84 afric0_mat(i,j,0)=afric_mat(i,j)
85 vfric0_mat(i,j,0)=vfric_mat(i,j)
86 prand0_mat(i,j,0)=prand_mat(i,j)
87 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
88 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
97 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
99 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
103 else if (lang.eq.1 .or. lang.eq.4) then
107 t_langsetup=MPI_Wtime()-tt0
110 t_langsetup=tcpu()-tt0
113 do itime=1,n_timestep
115 if (lang.gt.0 .and. surfarea .and.
116 & mod(itime,reset_fricmat).eq.0) then
117 if (lang.eq.2 .or. lang.eq.3) then
121 call sd_verlet_p_setup
123 call sd_verlet_ciccotti_setup
127 pfric0_mat(i,j,0)=pfric_mat(i,j)
128 afric0_mat(i,j,0)=afric_mat(i,j)
129 vfric0_mat(i,j,0)=vfric_mat(i,j)
130 prand0_mat(i,j,0)=prand_mat(i,j)
131 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
132 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
137 flag_stoch(i)=.false.
140 else if (lang.eq.1 .or. lang.eq.4) then
143 write (iout,'(a,i10)')
144 & "Friction matrix reset based on surface area, itime",itime
146 if (reset_vel .and. tbf .and. lang.eq.0
147 & .and. mod(itime,count_reset_vel).eq.0) then
149 write(iout,'(a,f20.2)')
150 & "Velocities reset to random values, time",totT
153 d_t_old(j,i)=d_t(j,i)
157 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
161 d_t(j,0)=d_t(j,0)-vcm(j)
164 kinetic_T=2.0d0/(dimen3*Rb)*EK
165 scalfac=dsqrt(T_bath/kinetic_T)
166 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
169 d_t_old(j,i)=scalfac*d_t(j,i)
175 c Time-reversible RESPA algorithm
176 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
177 call RESPA_step(itime)
179 c Variable time step algorithm.
180 call velverlet_step(itime)
184 call brown_step(itime)
186 print *,"Brown dynamics not here!"
188 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
194 if (mod(itime,ntwe).eq.0) call statout(itime)
207 if (itype(i).ne.10) then
210 v_work(ind)=d_t(j,i+nres)
215 write (66,'(80f10.5)')
216 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
220 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
222 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
224 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
227 if (mod(itime,ntwx).eq.0) then
228 write (tytul,'("time",f8.2)') totT
230 call pdbout(potE,tytul,ipdb)
235 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
236 open(irest2,file=rest2name,status='unknown')
237 write(irest2,*) totT,EK,potE,totE,t_bath
239 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
242 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
253 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
255 & 'MD calculations setup:',t_MDsetup,
256 & 'Energy & gradient evaluation:',t_enegrad,
257 & 'Stochastic MD setup:',t_langsetup,
258 & 'Stochastic MD step setup:',t_sdsetup,
260 write (iout,'(/28(1h=),a25,27(1h=))')
261 & ' End of MD calculation '
263 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
265 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
266 & " time_fricmatmult",time_fricmatmult," time_fsample ",
271 c-------------------------------------------------------------------------------
272 subroutine velverlet_step(itime)
273 c-------------------------------------------------------------------------------
274 c Perform a single velocity Verlet step; the time step can be rescaled if
275 c increments in accelerations exceed the threshold
276 c-------------------------------------------------------------------------------
277 implicit real*8 (a-h,o-z)
281 integer ierror,ierrcode
283 include 'COMMON.SETUP'
284 include 'COMMON.CONTROL'
288 include 'COMMON.LANGEVIN'
290 include 'COMMON.LANGEVIN.lang0'
292 include 'COMMON.CHAIN'
293 include 'COMMON.DERIV'
295 include 'COMMON.LOCAL'
296 include 'COMMON.INTERACT'
297 include 'COMMON.IOUNITS'
298 include 'COMMON.NAMES'
299 include 'COMMON.TIME1'
300 include 'COMMON.MUCA'
301 double precision vcm(3),incr(3)
302 double precision cm(3),L(3)
303 integer ilen,count,rstcount
306 integer maxcount_scale /20/
308 double precision stochforcvec(MAXRES6)
309 common /stochcalc/ stochforcvec
317 else if (lang.eq.2 .or. lang.eq.3) then
319 call stochastic_force(stochforcvec)
322 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
324 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
331 icount_scale=icount_scale+1
332 if (icount_scale.gt.maxcount_scale) then
334 & "ERROR: too many attempts at scaling down the time step. ",
335 & "amax=",amax,"epdrift=",epdrift,
336 & "damax=",damax,"edriftmax=",edriftmax,
340 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
344 c First step of the velocity Verlet algorithm
349 else if (lang.eq.3) then
351 call sd_verlet1_ciccotti
353 else if (lang.eq.1) then
358 c Build the chain from the newly calculated coordinates
360 if (rattle) call rattle1
362 if (large.and. mod(itime,ntwe).eq.0) then
363 write (iout,*) "Cartesian and internal coordinates: step 1"
368 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
369 & (dc(j,i+nres),j=1,3)
371 write (iout,*) "Accelerations"
373 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
374 & (d_a(j,i+nres),j=1,3)
376 write (iout,*) "Velocities, step 1"
378 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
379 & (d_t(j,i+nres),j=1,3)
388 c Calculate energy and forces
390 call etotal(potEcomp)
393 t_etotal=t_etotal+MPI_Wtime()-tt0
395 t_etotal=t_etotal+tcpu()-tt0
398 potE=potEcomp(0)-potEcomp(20)
400 c Get the new accelerations
403 t_enegrad=t_enegrad+MPI_Wtime()-tt0
405 t_enegrad=t_enegrad+tcpu()-tt0
407 c Determine maximum acceleration and scale down the timestep if needed
409 amax=amax/(itime_scal+1)**2
410 call predict_edrift(epdrift)
411 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
412 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
414 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
416 itime_scal=itime_scal+ifac_time
417 c fac_time=dmin1(damax/amax,0.5d0)
418 fac_time=0.5d0**ifac_time
419 d_time=d_time*fac_time
420 if (lang.eq.2 .or. lang.eq.3) then
422 c write (iout,*) "Calling sd_verlet_setup: 1"
423 c Rescale the stochastic forces and recalculate or restore
424 c the matrices of tinker integrator
425 if (itime_scal.gt.maxflag_stoch) then
426 if (large) write (iout,'(a,i5,a)')
427 & "Calculate matrices for stochastic step;",
428 & " itime_scal ",itime_scal
430 call sd_verlet_p_setup
432 call sd_verlet_ciccotti_setup
434 write (iout,'(2a,i3,a,i3,1h.)')
435 & "Warning: cannot store matrices for stochastic",
436 & " integration because the index",itime_scal,
437 & " is greater than",maxflag_stoch
438 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
439 & " integration Langevin algorithm for better efficiency."
440 else if (flag_stoch(itime_scal)) then
441 if (large) write (iout,'(a,i5,a,l1)')
442 & "Restore matrices for stochastic step; itime_scal ",
443 & itime_scal," flag ",flag_stoch(itime_scal)
446 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
447 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
448 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
449 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
450 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
451 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
455 if (large) write (iout,'(2a,i5,a,l1)')
456 & "Calculate & store matrices for stochastic step;",
457 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
459 call sd_verlet_p_setup
461 call sd_verlet_ciccotti_setup
463 flag_stoch(ifac_time)=.true.
466 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
467 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
468 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
469 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
470 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
471 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
475 fac_time=1.0d0/dsqrt(fac_time)
477 stochforcvec(i)=fac_time*stochforcvec(i)
480 else if (lang.eq.1) then
481 c Rescale the accelerations due to stochastic forces
482 fac_time=1.0d0/dsqrt(fac_time)
484 d_as_work(i)=d_as_work(i)*fac_time
487 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
488 & "itime",itime," Timestep scaled down to ",
489 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
491 c Second step of the velocity Verlet algorithm
496 else if (lang.eq.3) then
498 call sd_verlet2_ciccotti
500 else if (lang.eq.1) then
505 if (rattle) call rattle2
507 if (d_time.ne.d_time0) then
510 if (lang.eq.2 .or. lang.eq.3) then
511 if (large) write (iout,'(a)')
512 & "Restore original matrices for stochastic step"
513 c write (iout,*) "Calling sd_verlet_setup: 2"
514 c Restore the matrices of tinker integrator if the time step has been restored
517 pfric_mat(i,j)=pfric0_mat(i,j,0)
518 afric_mat(i,j)=afric0_mat(i,j,0)
519 vfric_mat(i,j)=vfric0_mat(i,j,0)
520 prand_mat(i,j)=prand0_mat(i,j,0)
521 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
522 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
531 c Calculate the kinetic and the total energy and the kinetic temperature
536 c write (iout,*) "step",itime," EK",EK," EK1",EK1
538 c Couple the system to Berendsen bath if needed
539 if (tbf .and. lang.eq.0) then
542 kinetic_T=2.0d0/(dimen3*Rb)*EK
543 c Backup the coordinates, velocities, and accelerations
547 d_t_old(j,i)=d_t(j,i)
548 d_a_old(j,i)=d_a(j,i)
552 if (mod(itime,ntwe).eq.0 .and. large) then
553 write (iout,*) "Velocities, step 2"
555 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
556 & (d_t(j,i+nres),j=1,3)
562 c-------------------------------------------------------------------------------
563 subroutine RESPA_step(itime)
564 c-------------------------------------------------------------------------------
565 c Perform a single RESPA step.
566 c-------------------------------------------------------------------------------
567 implicit real*8 (a-h,o-z)
571 integer IERROR,ERRCODE
573 include 'COMMON.SETUP'
574 include 'COMMON.CONTROL'
578 include 'COMMON.LANGEVIN'
580 include 'COMMON.LANGEVIN.lang0'
582 include 'COMMON.CHAIN'
583 include 'COMMON.DERIV'
585 include 'COMMON.LOCAL'
586 include 'COMMON.INTERACT'
587 include 'COMMON.IOUNITS'
588 include 'COMMON.NAMES'
589 include 'COMMON.TIME1'
590 double precision energia_short(0:n_ene),
591 & energia_long(0:n_ene)
592 double precision cm(3),L(3),vcm(3),incr(3)
593 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
594 & d_a_old0(3,0:maxres2)
595 integer ilen,count,rstcount
598 integer maxcount_scale /10/
600 double precision stochforcvec(MAXRES6)
601 common /stochcalc/ stochforcvec
604 common /cipiszcze/ itt
607 if (large.and. mod(itime,ntwe).eq.0) then
608 write (iout,*) "***************** RESPA itime",itime
609 write (iout,*) "Cartesian and internal coordinates: step 0"
611 call pdbout(0.0d0,"cipiszcze",iout)
613 write (iout,*) "Accelerations from long-range forces"
615 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
616 & (d_a(j,i+nres),j=1,3)
618 write (iout,*) "Velocities, step 0"
620 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
621 & (d_t(j,i+nres),j=1,3)
626 c Perform the initial RESPA step (increment velocities)
627 c write (iout,*) "*********************** RESPA ini"
630 if (mod(itime,ntwe).eq.0 .and. large) then
631 write (iout,*) "Velocities, end"
633 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
634 & (d_t(j,i+nres),j=1,3)
638 c Compute the short-range forces
644 C 7/2/2009 commented out
646 c call etotal_short(energia_short)
649 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
652 d_a(j,i)=d_a_short(j,i)
656 if (large.and. mod(itime,ntwe).eq.0) then
657 write (iout,*) "energia_short",energia_short(0)
658 write (iout,*) "Accelerations from short-range forces"
660 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
661 & (d_a(j,i+nres),j=1,3)
666 t_enegrad=t_enegrad+MPI_Wtime()-tt0
668 t_enegrad=t_enegrad+tcpu()-tt0
673 d_t_old(j,i)=d_t(j,i)
674 d_a_old(j,i)=d_a(j,i)
677 c 6/30/08 A-MTS: attempt at increasing the split number
680 dc_old0(j,i)=dc_old(j,i)
681 d_t_old0(j,i)=d_t_old(j,i)
682 d_a_old0(j,i)=d_a_old(j,i)
685 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
686 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
693 c write (iout,*) "itime",itime," ntime_split",ntime_split
694 c Split the time step
695 d_time=d_time0/ntime_split
696 c Perform the short-range RESPA steps (velocity Verlet increments of
697 c positions and velocities using short-range forces)
698 c write (iout,*) "*********************** RESPA split"
699 do itsplit=1,ntime_split
702 else if (lang.eq.2 .or. lang.eq.3) then
704 call stochastic_force(stochforcvec)
707 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
709 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
714 c First step of the velocity Verlet algorithm
719 else if (lang.eq.3) then
721 call sd_verlet1_ciccotti
723 else if (lang.eq.1) then
728 c Build the chain from the newly calculated coordinates
730 if (rattle) call rattle1
732 if (large.and. mod(itime,ntwe).eq.0) then
733 write (iout,*) "***** ITSPLIT",itsplit
734 write (iout,*) "Cartesian and internal coordinates: step 1"
735 call pdbout(0.0d0,"cipiszcze",iout)
738 write (iout,*) "Velocities, step 1"
740 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
741 & (d_t(j,i+nres),j=1,3)
750 c Calculate energy and forces
752 call etotal_short(energia_short)
755 t_eshort=t_eshort+MPI_Wtime()-tt0
757 t_eshort=t_eshort+tcpu()-tt0
761 c Get the new accelerations
763 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
766 d_a_short(j,i)=d_a(j,i)
770 if (large.and. mod(itime,ntwe).eq.0) then
771 write (iout,*)"energia_short",energia_short(0)
772 write (iout,*) "Accelerations from short-range forces"
774 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
775 & (d_a(j,i+nres),j=1,3)
780 c Determine maximum acceleration and scale down the timestep if needed
782 amax=amax/ntime_split**2
783 call predict_edrift(epdrift)
784 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
785 & write (iout,*) "amax",amax," damax",damax,
786 & " epdrift",epdrift," epdriftmax",epdriftmax
787 c Exit loop and try with increased split number if the change of
788 c acceleration is too big
789 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
790 if (ntime_split.lt.maxtime_split) then
792 ntime_split=ntime_split*2
795 dc_old(j,i)=dc_old0(j,i)
796 d_t_old(j,i)=d_t_old0(j,i)
797 d_a_old(j,i)=d_a_old0(j,i)
800 write (iout,*) "acceleration/energy drift too large",amax,
801 & epdrift," split increased to ",ntime_split," itime",itime,
806 & "Uh-hu. Bumpy landscape. Maximum splitting number",
808 & " already reached!!! Trying to carry on!"
812 t_enegrad=t_enegrad+MPI_Wtime()-tt0
814 t_enegrad=t_enegrad+tcpu()-tt0
816 c Second step of the velocity Verlet algorithm
821 else if (lang.eq.3) then
823 call sd_verlet2_ciccotti
825 else if (lang.eq.1) then
830 if (rattle) call rattle2
831 c Backup the coordinates, velocities, and accelerations
835 d_t_old(j,i)=d_t(j,i)
836 d_a_old(j,i)=d_a(j,i)
843 c Restore the time step
845 c Compute long-range forces
852 call etotal_long(energia_long)
855 t_elong=t_elong+MPI_Wtime()-tt0
857 t_elong=t_elong+tcpu()-tt0
863 t_enegrad=t_enegrad+MPI_Wtime()-tt0
865 t_enegrad=t_enegrad+tcpu()-tt0
867 c Compute accelerations from long-range forces
869 if (large.and. mod(itime,ntwe).eq.0) then
870 write (iout,*) "energia_long",energia_long(0)
871 write (iout,*) "Cartesian and internal coordinates: step 2"
873 call pdbout(0.0d0,"cipiszcze",iout)
875 write (iout,*) "Accelerations from long-range forces"
877 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
878 & (d_a(j,i+nres),j=1,3)
880 write (iout,*) "Velocities, step 2"
882 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
883 & (d_t(j,i+nres),j=1,3)
887 c Compute the final RESPA step (increment velocities)
888 c write (iout,*) "*********************** RESPA fin"
890 c Compute the complete potential energy
892 potEcomp(i)=energia_short(i)+energia_long(i)
894 potE=potEcomp(0)-potEcomp(20)
895 c potE=energia_short(0)+energia_long(0)
897 c Calculate the kinetic and the total energy and the kinetic temperature
900 c Couple the system to Berendsen bath if needed
901 if (tbf .and. lang.eq.0) then
904 kinetic_T=2.0d0/(dimen3*Rb)*EK
905 c Backup the coordinates, velocities, and accelerations
907 if (mod(itime,ntwe).eq.0 .and. large) then
908 write (iout,*) "Velocities, end"
910 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
911 & (d_t(j,i+nres),j=1,3)
917 c---------------------------------------------------------------------
919 c First and last RESPA step (incrementing velocities using long-range
921 implicit real*8 (a-h,o-z)
923 include 'COMMON.CONTROL'
926 include 'COMMON.CHAIN'
927 include 'COMMON.DERIV'
929 include 'COMMON.LOCAL'
930 include 'COMMON.INTERACT'
931 include 'COMMON.IOUNITS'
932 include 'COMMON.NAMES'
934 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
938 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
942 if (itype(i).ne.10) then
945 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
951 c-----------------------------------------------------------------
953 c Applying velocity Verlet algorithm - step 1 to coordinates
954 implicit real*8 (a-h,o-z)
956 include 'COMMON.CONTROL'
959 include 'COMMON.CHAIN'
960 include 'COMMON.DERIV'
962 include 'COMMON.LOCAL'
963 include 'COMMON.INTERACT'
964 include 'COMMON.IOUNITS'
965 include 'COMMON.NAMES'
966 double precision adt,adt2
969 write (iout,*) "VELVERLET1 START: DC"
971 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
972 & (dc(j,i+nres),j=1,3)
976 adt=d_a_old(j,0)*d_time
978 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
979 d_t_new(j,0)=d_t_old(j,0)+adt2
980 d_t(j,0)=d_t_old(j,0)+adt
984 adt=d_a_old(j,i)*d_time
986 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
987 d_t_new(j,i)=d_t_old(j,i)+adt2
988 d_t(j,i)=d_t_old(j,i)+adt
992 if (itype(i).ne.10) then
995 adt=d_a_old(j,inres)*d_time
997 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
998 d_t_new(j,inres)=d_t_old(j,inres)+adt2
999 d_t(j,inres)=d_t_old(j,inres)+adt
1004 write (iout,*) "VELVERLET1 END: DC"
1006 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1007 & (dc(j,i+nres),j=1,3)
1012 c---------------------------------------------------------------------
1014 c Step 2 of the velocity Verlet algorithm: update velocities
1015 implicit real*8 (a-h,o-z)
1016 include 'DIMENSIONS'
1017 include 'COMMON.CONTROL'
1018 include 'COMMON.VAR'
1020 include 'COMMON.CHAIN'
1021 include 'COMMON.DERIV'
1022 include 'COMMON.GEO'
1023 include 'COMMON.LOCAL'
1024 include 'COMMON.INTERACT'
1025 include 'COMMON.IOUNITS'
1026 include 'COMMON.NAMES'
1028 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1032 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1036 if (itype(i).ne.10) then
1039 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1045 c-----------------------------------------------------------------
1046 subroutine sddir_precalc
1047 c Applying velocity Verlet algorithm - step 1 to coordinates
1048 implicit real*8 (a-h,o-z)
1049 include 'DIMENSIONS'
1053 include 'COMMON.CONTROL'
1054 include 'COMMON.VAR'
1057 include 'COMMON.LANGEVIN'
1059 include 'COMMON.LANGEVIN.lang0'
1061 include 'COMMON.CHAIN'
1062 include 'COMMON.DERIV'
1063 include 'COMMON.GEO'
1064 include 'COMMON.LOCAL'
1065 include 'COMMON.INTERACT'
1066 include 'COMMON.IOUNITS'
1067 include 'COMMON.NAMES'
1068 include 'COMMON.TIME1'
1069 double precision stochforcvec(MAXRES6)
1070 common /stochcalc/ stochforcvec
1072 c Compute friction and stochastic forces
1076 time_fric=time_fric+MPI_Wtime()-time00
1078 call stochastic_force(stochforcvec)
1079 time_stoch=time_stoch+MPI_Wtime()-time00
1081 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1082 c forces (d_as_work)
1084 call ginv_mult(fric_work, d_af_work)
1085 call ginv_mult(stochforcvec, d_as_work)
1088 c---------------------------------------------------------------------
1089 subroutine sddir_verlet1
1090 c Applying velocity Verlet algorithm - step 1 to velocities
1091 implicit real*8 (a-h,o-z)
1092 include 'DIMENSIONS'
1093 include 'COMMON.CONTROL'
1094 include 'COMMON.VAR'
1097 include 'COMMON.LANGEVIN'
1099 include 'COMMON.LANGEVIN.lang0'
1101 include 'COMMON.CHAIN'
1102 include 'COMMON.DERIV'
1103 include 'COMMON.GEO'
1104 include 'COMMON.LOCAL'
1105 include 'COMMON.INTERACT'
1106 include 'COMMON.IOUNITS'
1107 include 'COMMON.NAMES'
1108 c Revised 3/31/05 AL: correlation between random contributions to
1109 c position and velocity increments included.
1110 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1111 double precision adt,adt2
1113 c Add the contribution from BOTH friction and stochastic force to the
1114 c coordinates, but ONLY the contribution from the friction forces to velocities
1117 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1118 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1119 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1120 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1121 d_t(j,0)=d_t_old(j,0)+adt
1126 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1127 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1128 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1129 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1130 d_t(j,i)=d_t_old(j,i)+adt
1135 if (itype(i).ne.10) then
1138 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1139 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1140 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1141 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1142 d_t(j,inres)=d_t_old(j,inres)+adt
1149 c---------------------------------------------------------------------
1150 subroutine sddir_verlet2
1151 c Calculating the adjusted velocities for accelerations
1152 implicit real*8 (a-h,o-z)
1153 include 'DIMENSIONS'
1154 include 'COMMON.CONTROL'
1155 include 'COMMON.VAR'
1158 include 'COMMON.LANGEVIN'
1160 include 'COMMON.LANGEVIN.lang0'
1162 include 'COMMON.CHAIN'
1163 include 'COMMON.DERIV'
1164 include 'COMMON.GEO'
1165 include 'COMMON.LOCAL'
1166 include 'COMMON.INTERACT'
1167 include 'COMMON.IOUNITS'
1168 include 'COMMON.NAMES'
1169 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1170 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1171 c Revised 3/31/05 AL: correlation between random contributions to
1172 c position and velocity increments included.
1173 c The correlation coefficients are calculated at low-friction limit.
1174 c Also, friction forces are now not calculated with new velocities.
1176 c call friction_force
1177 call stochastic_force(stochforcvec)
1179 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1180 c forces (d_as_work)
1182 call ginv_mult(stochforcvec, d_as_work1)
1188 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1189 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1194 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1195 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1200 if (itype(i).ne.10) then
1203 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1204 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1205 & +cos60*d_as_work1(ind+j))*d_time
1212 c---------------------------------------------------------------------
1213 subroutine max_accel
1215 c Find the maximum difference in the accelerations of the the sites
1216 c at the beginning and the end of the time step.
1218 implicit real*8 (a-h,o-z)
1219 include 'DIMENSIONS'
1220 include 'COMMON.CONTROL'
1221 include 'COMMON.VAR'
1223 include 'COMMON.CHAIN'
1224 include 'COMMON.DERIV'
1225 include 'COMMON.GEO'
1226 include 'COMMON.LOCAL'
1227 include 'COMMON.INTERACT'
1228 include 'COMMON.IOUNITS'
1229 double precision aux(3),accel(3),accel_old(3),dacc
1231 c aux(j)=d_a(j,0)-d_a_old(j,0)
1232 accel_old(j)=d_a_old(j,0)
1239 c 7/3/08 changed to asymmetric difference
1241 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1242 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1243 accel(j)=accel(j)+0.5d0*d_a(j,i)
1244 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1245 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1246 dacc=dabs(accel(j)-accel_old(j))
1247 if (dacc.gt.amax) amax=dacc
1255 accel_old(j)=d_a_old(j,0)
1260 accel_old(j)=accel_old(j)+d_a_old(j,1)
1261 accel(j)=accel(j)+d_a(j,1)
1265 if (itype(i).ne.10) then
1267 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1268 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1269 accel(j)=accel(j)+d_a(j,i+nres)
1273 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1274 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1275 dacc=dabs(accel(j)-accel_old(j))
1276 if (dacc.gt.amax) amax=dacc
1280 accel_old(j)=accel_old(j)+d_a_old(j,i)
1281 accel(j)=accel(j)+d_a(j,i)
1282 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1287 c---------------------------------------------------------------------
1288 subroutine predict_edrift(epdrift)
1290 c Predict the drift of the potential energy
1292 implicit real*8 (a-h,o-z)
1293 include 'DIMENSIONS'
1294 include 'COMMON.CONTROL'
1295 include 'COMMON.VAR'
1297 include 'COMMON.CHAIN'
1298 include 'COMMON.DERIV'
1299 include 'COMMON.GEO'
1300 include 'COMMON.LOCAL'
1301 include 'COMMON.INTERACT'
1302 include 'COMMON.IOUNITS'
1303 include 'COMMON.MUCA'
1304 double precision epdrift,epdriftij
1305 c Drift of the potential energy
1311 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1312 if (lmuca) epdriftij=epdriftij*factor
1313 c write (iout,*) "back",i,j,epdriftij
1314 if (epdriftij.gt.epdrift) epdrift=epdriftij
1318 if (itype(i).ne.10) then
1321 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1322 if (lmuca) epdriftij=epdriftij*factor
1323 c write (iout,*) "side",i,j,epdriftij
1324 if (epdriftij.gt.epdrift) epdrift=epdriftij
1328 epdrift=0.5d0*epdrift*d_time*d_time
1329 c write (iout,*) "epdrift",epdrift
1332 c-----------------------------------------------------------------------
1333 subroutine verlet_bath
1335 c Coupling to the thermostat by using the Berendsen algorithm
1337 implicit real*8 (a-h,o-z)
1338 include 'DIMENSIONS'
1339 include 'COMMON.CONTROL'
1340 include 'COMMON.VAR'
1342 include 'COMMON.CHAIN'
1343 include 'COMMON.DERIV'
1344 include 'COMMON.GEO'
1345 include 'COMMON.LOCAL'
1346 include 'COMMON.INTERACT'
1347 include 'COMMON.IOUNITS'
1348 include 'COMMON.NAMES'
1349 double precision T_half,fact
1351 T_half=2.0d0/(dimen3*Rb)*EK
1352 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1353 c write(iout,*) "T_half", T_half
1354 c write(iout,*) "EK", EK
1355 c write(iout,*) "fact", fact
1357 d_t(j,0)=fact*d_t(j,0)
1361 d_t(j,i)=fact*d_t(j,i)
1365 if (itype(i).ne.10) then
1368 d_t(j,inres)=fact*d_t(j,inres)
1374 c---------------------------------------------------------
1376 c Set up the initial conditions of a MD simulation
1377 implicit real*8 (a-h,o-z)
1378 include 'DIMENSIONS'
1382 integer IERROR,ERRCODE
1384 include 'COMMON.SETUP'
1385 include 'COMMON.CONTROL'
1386 include 'COMMON.VAR'
1389 include 'COMMON.LANGEVIN'
1391 include 'COMMON.LANGEVIN.lang0'
1393 include 'COMMON.CHAIN'
1394 include 'COMMON.DERIV'
1395 include 'COMMON.GEO'
1396 include 'COMMON.LOCAL'
1397 include 'COMMON.INTERACT'
1398 include 'COMMON.IOUNITS'
1399 include 'COMMON.NAMES'
1400 include 'COMMON.REMD'
1401 real*8 energia_long(0:n_ene),
1402 & energia_short(0:n_ene),vcm(3),incr(3)
1403 double precision cm(3),L(3),xv,sigv,lowb,highb
1404 double precision varia(maxvar)
1412 c write(iout,*) "d_time", d_time
1413 c Compute the standard deviations of stochastic forces for Langevin dynamics
1414 c if the friction coefficients do not depend on surface area
1415 if (lang.gt.0 .and. .not.surfarea) then
1417 stdforcp(i)=stdfp*dsqrt(gamp)
1420 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
1423 c Open the pdb file for snapshotshots
1426 if (ilen(tmpdir).gt.0)
1427 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1428 & liczba(:ilen(liczba))//".pdb")
1430 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1434 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1435 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1436 & liczba(:ilen(liczba))//".x")
1437 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1440 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1441 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1442 & liczba(:ilen(liczba))//".cx")
1443 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1449 if (ilen(tmpdir).gt.0)
1450 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1451 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1453 if (ilen(tmpdir).gt.0)
1454 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1455 cartname=prefix(:ilen(prefix))//"_MD.cx"
1459 write (qstr,'(256(1h ))')
1462 iq = qinfrag(i,iset)*10
1463 iw = wfrag(i,iset)/100
1465 if(me.eq.king.or..not.out1file)
1466 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1467 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1472 iq = qinpair(i,iset)*10
1473 iw = wpair(i,iset)/100
1475 if(me.eq.king.or..not.out1file)
1476 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1477 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1481 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1483 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1485 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1487 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1491 if (restart1file) then
1493 & inquire(file=mremd_rst_name,exist=file_exist)
1494 write (*,*) me," Before broadcast: file_exist",file_exist
1495 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1497 write (*,*) me," After broadcast: file_exist",file_exist
1498 c inquire(file=mremd_rst_name,exist=file_exist)
1499 if(me.eq.king.or..not.out1file)
1500 & write(iout,*) "Initial state read by master and distributed"
1502 if (ilen(tmpdir).gt.0)
1503 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1504 & //liczba(:ilen(liczba))//'.rst')
1505 inquire(file=rest2name,exist=file_exist)
1508 if(.not.restart1file) then
1509 if(me.eq.king.or..not.out1file)
1510 & write(iout,*) "Initial state will be read from file ",
1511 & rest2name(:ilen(rest2name))
1514 call rescale_weights(t_bath)
1516 if(me.eq.king.or..not.out1file)then
1517 if (restart1file) then
1518 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1521 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1524 write(iout,*) "Initial velocities randomly generated"
1530 c Generate initial velocities
1531 if(me.eq.king.or..not.out1file)
1532 & write(iout,*) "Initial velocities randomly generated"
1536 c rest2name = prefix(:ilen(prefix))//'.rst'
1537 if(me.eq.king.or..not.out1file)then
1538 write (iout,*) "Initial velocities"
1540 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1541 & (d_t(j,i+nres),j=1,3)
1543 c Zeroing the total angular momentum of the system
1544 write(iout,*) "Calling the zero-angular
1545 & momentum subroutine"
1548 c Getting the potential energy and forces and velocities and accelerations
1550 c write (iout,*) "velocity of the center of the mass:"
1551 c write (iout,*) (vcm(j),j=1,3)
1553 d_t(j,0)=d_t(j,0)-vcm(j)
1555 c Removing the velocity of the center of mass
1557 if(me.eq.king.or..not.out1file)then
1558 write (iout,*) "vcm right after adjustment:"
1559 write (iout,*) (vcm(j),j=1,3)
1563 if(iranconf.ne.0) then
1565 print *, 'Calling OVERLAP_SC'
1566 call overlap_sc(fail)
1570 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1571 print *,'SC_move',nft_sc,etot
1572 if(me.eq.king.or..not.out1file)
1573 & write(iout,*) 'SC_move',nft_sc,etot
1577 print *, 'Calling MINIM_DC'
1578 call minim_dc(etot,iretcode,nfun)
1580 call geom_to_var(nvar,varia)
1581 print *,'Calling MINIMIZE.'
1582 call minimize(etot,varia,iretcode,nfun)
1583 call var_to_geom(nvar,varia)
1585 if(me.eq.king.or..not.out1file)
1586 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1589 call chainbuild_cart
1592 call verlet_bath(EK)
1594 kinetic_T=2.0d0/(dimen3*Rb)*EK
1595 if(me.eq.king.or..not.out1file)then
1605 call etotal(potEcomp)
1608 t_etotal=t_etotal+MPI_Wtime()-tt0
1610 t_etotal=t_etotal+tcpu()-tt0
1617 if (amax*d_time .gt. dvmax) then
1618 d_time=d_time*dvmax/amax
1619 if(me.eq.king.or..not.out1file) write (iout,*)
1620 & "Time step reduced to",d_time,
1621 & " because of too large initial acceleration."
1623 if(me.eq.king.or..not.out1file)then
1624 write(iout,*) "Potential energy and its components"
1625 call enerprint(potEcomp)
1626 c write(iout,*) (potEcomp(i),i=0,n_ene)
1628 potE=potEcomp(0)-potEcomp(20)
1631 if (ntwe.ne.0) call statout(itime)
1632 if(me.eq.king.or..not.out1file)
1633 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1634 & " Kinetic energy",EK," potential energy",potE,
1635 & " total energy",totE," maximum acceleration ",
1638 write (iout,*) "Initial coordinates"
1640 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1641 & (c(j,i+nres),j=1,3)
1643 write (iout,*) "Initial dC"
1645 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1646 & (dc(j,i+nres),j=1,3)
1648 write (iout,*) "Initial velocities"
1649 write (iout,"(13x,' backbone ',23x,' side chain')")
1651 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1652 & (d_t(j,i+nres),j=1,3)
1654 write (iout,*) "Initial accelerations"
1656 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1657 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1658 & (d_a(j,i+nres),j=1,3)
1664 d_t_old(j,i)=d_t(j,i)
1665 d_a_old(j,i)=d_a(j,i)
1667 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1676 call etotal_short(energia_short)
1679 t_eshort=t_eshort+MPI_Wtime()-tt0
1681 t_eshort=t_eshort+tcpu()-tt0
1686 if(.not.out1file .and. large) then
1687 write (iout,*) "energia_long",energia_long(0),
1688 & " energia_short",energia_short(0),
1689 & " total",energia_long(0)+energia_short(0)
1690 write (iout,*) "Initial fast-force accelerations"
1692 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1693 & (d_a(j,i+nres),j=1,3)
1696 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1699 d_a_short(j,i)=d_a(j,i)
1708 call etotal_long(energia_long)
1711 t_elong=t_elong+MPI_Wtime()-tt0
1713 t_elong=t_elong+tcpu()-tt0
1718 if(.not.out1file .and. large) then
1719 write (iout,*) "energia_long",energia_long(0)
1720 write (iout,*) "Initial slow-force accelerations"
1722 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1723 & (d_a(j,i+nres),j=1,3)
1727 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1729 t_enegrad=t_enegrad+tcpu()-tt0
1734 c-----------------------------------------------------------
1735 subroutine random_vel
1736 implicit real*8 (a-h,o-z)
1737 include 'DIMENSIONS'
1738 include 'COMMON.CONTROL'
1739 include 'COMMON.VAR'
1742 include 'COMMON.LANGEVIN'
1744 include 'COMMON.LANGEVIN.lang0'
1746 include 'COMMON.CHAIN'
1747 include 'COMMON.DERIV'
1748 include 'COMMON.GEO'
1749 include 'COMMON.LOCAL'
1750 include 'COMMON.INTERACT'
1751 include 'COMMON.IOUNITS'
1752 include 'COMMON.NAMES'
1753 include 'COMMON.TIME1'
1754 double precision xv,sigv,lowb,highb
1755 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1756 c First generate velocities in the eigenspace of the G matrix
1757 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1764 sigv=dsqrt((Rb*t_bath)/geigen(i))
1767 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1768 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1769 c & " d_t_work_new",d_t_work_new(ii)
1778 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1781 c write (iout,*) "Ek from eigenvectors",Ek1
1783 c Transform velocities to UNRES coordinate space
1789 d_t_work(ind)=d_t_work(ind)
1790 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1792 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1796 c Transfer to the d_t vector
1798 d_t(j,0)=d_t_work(j)
1804 d_t(j,i)=d_t_work(ind)
1808 if (itype(i).ne.10) then
1811 d_t(j,i+nres)=d_t_work(ind)
1816 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1817 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1822 c-----------------------------------------------------------
1823 subroutine sd_verlet_p_setup
1824 c Sets up the parameters of stochastic Verlet algorithm
1825 implicit real*8 (a-h,o-z)
1826 include 'DIMENSIONS'
1830 include 'COMMON.CONTROL'
1831 include 'COMMON.VAR'
1834 include 'COMMON.LANGEVIN'
1836 include 'COMMON.LANGEVIN.lang0'
1838 include 'COMMON.CHAIN'
1839 include 'COMMON.DERIV'
1840 include 'COMMON.GEO'
1841 include 'COMMON.LOCAL'
1842 include 'COMMON.INTERACT'
1843 include 'COMMON.IOUNITS'
1844 include 'COMMON.NAMES'
1845 include 'COMMON.TIME1'
1846 double precision emgdt(MAXRES6),
1847 & pterm,vterm,rho,rhoc,vsig,
1848 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1849 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1850 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1851 logical lprn /.false./
1852 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1853 double precision ktm
1860 c AL 8/17/04 Code adapted from tinker
1862 c Get the frictional and random terms for stochastic dynamics in the
1863 c eigenspace of mass-scaled UNRES friction matrix
1866 gdt = fricgam(i) * d_time
1868 c Stochastic dynamics reduces to simple MD for zero friction
1870 if (gdt .le. zero) then
1871 pfric_vec(i) = 1.0d0
1872 vfric_vec(i) = d_time
1873 afric_vec(i) = 0.5d0 * d_time * d_time
1874 prand_vec(i) = 0.0d0
1875 vrand_vec1(i) = 0.0d0
1876 vrand_vec2(i) = 0.0d0
1878 c Analytical expressions when friction coefficient is large
1881 if (gdt .ge. gdt_radius) then
1884 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1885 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1886 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1887 vterm = 1.0d0 - egdt**2
1888 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1890 c Use series expansions when friction coefficient is small
1901 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1902 & - gdt5/120.0d0 + gdt6/720.0d0
1903 & - gdt7/5040.0d0 + gdt8/40320.0d0
1904 & - gdt9/362880.0d0) / fricgam(i)**2
1905 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1906 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1907 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1908 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1909 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1910 & + 127.0d0*gdt9/90720.0d0
1911 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1912 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1913 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1914 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1915 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1916 & - 17.0d0*gdt2/1280.0d0
1917 & + 17.0d0*gdt3/6144.0d0
1918 & + 40967.0d0*gdt4/34406400.0d0
1919 & - 57203.0d0*gdt5/275251200.0d0
1920 & - 1429487.0d0*gdt6/13212057600.0d0)
1923 c Compute the scaling factors of random terms for the nonzero friction case
1925 ktm = 0.5d0*d_time/fricgam(i)
1926 psig = dsqrt(ktm*pterm) / fricgam(i)
1927 vsig = dsqrt(ktm*vterm)
1928 rhoc = dsqrt(1.0d0 - rho*rho)
1930 vrand_vec1(i) = vsig * rho
1931 vrand_vec2(i) = vsig * rhoc
1936 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1939 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1940 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1944 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1947 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1948 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1949 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1950 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1951 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1952 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1955 t_sdsetup=t_sdsetup+MPI_Wtime()
1957 t_sdsetup=t_sdsetup+tcpu()-tt0
1961 c-------------------------------------------------------------
1962 subroutine eigtransf1(n,ndim,ab,d,c)
1965 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
1971 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
1977 c-------------------------------------------------------------
1978 subroutine eigtransf(n,ndim,a,b,d,c)
1981 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
1987 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
1993 c-------------------------------------------------------------
1994 subroutine sd_verlet1
1995 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
1996 implicit real*8 (a-h,o-z)
1997 include 'DIMENSIONS'
1998 include 'COMMON.CONTROL'
1999 include 'COMMON.VAR'
2002 include 'COMMON.LANGEVIN'
2004 include 'COMMON.LANGEVIN.lang0'
2006 include 'COMMON.CHAIN'
2007 include 'COMMON.DERIV'
2008 include 'COMMON.GEO'
2009 include 'COMMON.LOCAL'
2010 include 'COMMON.INTERACT'
2011 include 'COMMON.IOUNITS'
2012 include 'COMMON.NAMES'
2013 double precision stochforcvec(MAXRES6)
2014 common /stochcalc/ stochforcvec
2015 logical lprn /.false./
2017 c write (iout,*) "dc_old"
2019 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2020 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2023 dc_work(j)=dc_old(j,0)
2024 d_t_work(j)=d_t_old(j,0)
2025 d_a_work(j)=d_a_old(j,0)
2030 dc_work(ind+j)=dc_old(j,i)
2031 d_t_work(ind+j)=d_t_old(j,i)
2032 d_a_work(ind+j)=d_a_old(j,i)
2037 if (itype(i).ne.10) then
2039 dc_work(ind+j)=dc_old(j,i+nres)
2040 d_t_work(ind+j)=d_t_old(j,i+nres)
2041 d_a_work(ind+j)=d_a_old(j,i+nres)
2049 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2053 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2054 & vfric_mat(i,j),afric_mat(i,j),
2055 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2063 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2064 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2065 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2066 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2068 d_t_work_new(i)=ddt1+0.5d0*ddt2
2069 d_t_work(i)=ddt1+ddt2
2074 d_t(j,0)=d_t_work(j)
2079 dc(j,i)=dc_work(ind+j)
2080 d_t(j,i)=d_t_work(ind+j)
2085 if (itype(i).ne.10) then
2088 dc(j,inres)=dc_work(ind+j)
2089 d_t(j,inres)=d_t_work(ind+j)
2096 c--------------------------------------------------------------------------
2097 subroutine sd_verlet2
2098 c Calculating the adjusted velocities for accelerations
2099 implicit real*8 (a-h,o-z)
2100 include 'DIMENSIONS'
2101 include 'COMMON.CONTROL'
2102 include 'COMMON.VAR'
2105 include 'COMMON.LANGEVIN'
2107 include 'COMMON.LANGEVIN.lang0'
2109 include 'COMMON.CHAIN'
2110 include 'COMMON.DERIV'
2111 include 'COMMON.GEO'
2112 include 'COMMON.LOCAL'
2113 include 'COMMON.INTERACT'
2114 include 'COMMON.IOUNITS'
2115 include 'COMMON.NAMES'
2116 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2117 common /stochcalc/ stochforcvec
2119 c Compute the stochastic forces which contribute to velocity change
2121 call stochastic_force(stochforcvecV)
2128 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2129 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2130 & vrand_mat2(i,j)*stochforcvecV(j)
2132 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2136 d_t(j,0)=d_t_work(j)
2141 d_t(j,i)=d_t_work(ind+j)
2146 if (itype(i).ne.10) then
2149 d_t(j,inres)=d_t_work(ind+j)
2156 c-----------------------------------------------------------
2157 subroutine sd_verlet_ciccotti_setup
2158 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2160 implicit real*8 (a-h,o-z)
2161 include 'DIMENSIONS'
2165 include 'COMMON.CONTROL'
2166 include 'COMMON.VAR'
2169 include 'COMMON.LANGEVIN'
2171 include 'COMMON.LANGEVIN.lang0'
2173 include 'COMMON.CHAIN'
2174 include 'COMMON.DERIV'
2175 include 'COMMON.GEO'
2176 include 'COMMON.LOCAL'
2177 include 'COMMON.INTERACT'
2178 include 'COMMON.IOUNITS'
2179 include 'COMMON.NAMES'
2180 include 'COMMON.TIME1'
2181 double precision emgdt(MAXRES6),
2182 & pterm,vterm,rho,rhoc,vsig,
2183 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2184 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2185 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2186 logical lprn /.false./
2187 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2188 double precision ktm
2195 c AL 8/17/04 Code adapted from tinker
2197 c Get the frictional and random terms for stochastic dynamics in the
2198 c eigenspace of mass-scaled UNRES friction matrix
2201 write (iout,*) "i",i," fricgam",fricgam(i)
2202 gdt = fricgam(i) * d_time
2204 c Stochastic dynamics reduces to simple MD for zero friction
2206 if (gdt .le. zero) then
2207 pfric_vec(i) = 1.0d0
2208 vfric_vec(i) = d_time
2209 afric_vec(i) = 0.5d0*d_time*d_time
2210 prand_vec(i) = afric_vec(i)
2211 vrand_vec2(i) = vfric_vec(i)
2213 c Analytical expressions when friction coefficient is large
2218 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2219 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2220 prand_vec(i) = afric_vec(i)
2221 vrand_vec2(i) = vfric_vec(i)
2223 c Compute the scaling factors of random terms for the nonzero friction case
2225 c ktm = 0.5d0*d_time/fricgam(i)
2226 c psig = dsqrt(ktm*pterm) / fricgam(i)
2227 c vsig = dsqrt(ktm*vterm)
2228 c prand_vec(i) = psig*afric_vec(i)
2229 c vrand_vec2(i) = vsig*vfric_vec(i)
2234 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2237 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2238 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2242 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2244 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2245 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2246 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2247 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2248 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2250 t_sdsetup=t_sdsetup+MPI_Wtime()
2252 t_sdsetup=t_sdsetup+tcpu()-tt0
2256 c-------------------------------------------------------------
2257 subroutine sd_verlet1_ciccotti
2258 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2259 implicit real*8 (a-h,o-z)
2260 include 'DIMENSIONS'
2264 include 'COMMON.CONTROL'
2265 include 'COMMON.VAR'
2268 include 'COMMON.LANGEVIN'
2270 include 'COMMON.LANGEVIN.lang0'
2272 include 'COMMON.CHAIN'
2273 include 'COMMON.DERIV'
2274 include 'COMMON.GEO'
2275 include 'COMMON.LOCAL'
2276 include 'COMMON.INTERACT'
2277 include 'COMMON.IOUNITS'
2278 include 'COMMON.NAMES'
2279 double precision stochforcvec(MAXRES6)
2280 common /stochcalc/ stochforcvec
2281 logical lprn /.false./
2283 c write (iout,*) "dc_old"
2285 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2286 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2289 dc_work(j)=dc_old(j,0)
2290 d_t_work(j)=d_t_old(j,0)
2291 d_a_work(j)=d_a_old(j,0)
2296 dc_work(ind+j)=dc_old(j,i)
2297 d_t_work(ind+j)=d_t_old(j,i)
2298 d_a_work(ind+j)=d_a_old(j,i)
2303 if (itype(i).ne.10) then
2305 dc_work(ind+j)=dc_old(j,i+nres)
2306 d_t_work(ind+j)=d_t_old(j,i+nres)
2307 d_a_work(ind+j)=d_a_old(j,i+nres)
2316 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2320 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2321 & vfric_mat(i,j),afric_mat(i,j),
2322 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2330 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2331 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2332 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2333 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2335 d_t_work_new(i)=ddt1+0.5d0*ddt2
2336 d_t_work(i)=ddt1+ddt2
2341 d_t(j,0)=d_t_work(j)
2346 dc(j,i)=dc_work(ind+j)
2347 d_t(j,i)=d_t_work(ind+j)
2352 if (itype(i).ne.10) then
2355 dc(j,inres)=dc_work(ind+j)
2356 d_t(j,inres)=d_t_work(ind+j)
2363 c--------------------------------------------------------------------------
2364 subroutine sd_verlet2_ciccotti
2365 c Calculating the adjusted velocities for accelerations
2366 implicit real*8 (a-h,o-z)
2367 include 'DIMENSIONS'
2368 include 'COMMON.CONTROL'
2369 include 'COMMON.VAR'
2372 include 'COMMON.LANGEVIN'
2374 include 'COMMON.LANGEVIN.lang0'
2376 include 'COMMON.CHAIN'
2377 include 'COMMON.DERIV'
2378 include 'COMMON.GEO'
2379 include 'COMMON.LOCAL'
2380 include 'COMMON.INTERACT'
2381 include 'COMMON.IOUNITS'
2382 include 'COMMON.NAMES'
2383 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2384 common /stochcalc/ stochforcvec
2386 c Compute the stochastic forces which contribute to velocity change
2388 call stochastic_force(stochforcvecV)
2395 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2396 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2397 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2399 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2403 d_t(j,0)=d_t_work(j)
2408 d_t(j,i)=d_t_work(ind+j)
2413 if (itype(i).ne.10) then
2416 d_t(j,inres)=d_t_work(ind+j)