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 include 'COMMON.HAIRPIN'
29 double precision cm(3),L(3),vcm(3)
31 double precision v_work(maxres6),v_transf(maxres6)
41 if (ilen(tmpdir).gt.0)
42 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
43 & //liczba(:ilen(liczba))//'.rst')
45 if (ilen(tmpdir).gt.0)
46 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
53 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
59 c Determine the inverse of the inertia matrix.
60 call setup_MD_matrices
64 t_MDsetup = MPI_Wtime()-tt0
66 t_MDsetup = tcpu()-tt0
69 c Entering the MD loop
75 if (lang.eq.2 .or. lang.eq.3) then
79 call sd_verlet_p_setup
81 call sd_verlet_ciccotti_setup
85 pfric0_mat(i,j,0)=pfric_mat(i,j)
86 afric0_mat(i,j,0)=afric_mat(i,j)
87 vfric0_mat(i,j,0)=vfric_mat(i,j)
88 prand0_mat(i,j,0)=prand_mat(i,j)
89 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
90 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
99 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
101 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
105 else if (lang.eq.1 .or. lang.eq.4) then
109 t_langsetup=MPI_Wtime()-tt0
112 t_langsetup=tcpu()-tt0
115 do itime=1,n_timestep
117 if (large.and. mod(itime,ntwe).eq.0)
118 & write (iout,*) "itime",itime
120 if (lang.gt.0 .and. surfarea .and.
121 & mod(itime,reset_fricmat).eq.0) then
122 if (lang.eq.2 .or. lang.eq.3) then
126 call sd_verlet_p_setup
128 call sd_verlet_ciccotti_setup
132 pfric0_mat(i,j,0)=pfric_mat(i,j)
133 afric0_mat(i,j,0)=afric_mat(i,j)
134 vfric0_mat(i,j,0)=vfric_mat(i,j)
135 prand0_mat(i,j,0)=prand_mat(i,j)
136 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
137 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
142 flag_stoch(i)=.false.
145 else if (lang.eq.1 .or. lang.eq.4) then
148 write (iout,'(a,i10)')
149 & "Friction matrix reset based on surface area, itime",itime
151 if (reset_vel .and. tbf .and. lang.eq.0
152 & .and. mod(itime,count_reset_vel).eq.0) then
154 write(iout,'(a,f20.2)')
155 & "Velocities reset to random values, time",totT
158 d_t_old(j,i)=d_t(j,i)
162 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
166 d_t(j,0)=d_t(j,0)-vcm(j)
169 kinetic_T=2.0d0/(dimen3*Rb)*EK
170 scalfac=dsqrt(T_bath/kinetic_T)
171 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
174 d_t_old(j,i)=scalfac*d_t(j,i)
180 c Time-reversible RESPA algorithm
181 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
182 call RESPA_step(itime)
184 c Variable time step algorithm.
185 call velverlet_step(itime)
189 call brown_step(itime)
191 print *,"Brown dynamics not here!"
193 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
199 if (mod(itime,ntwe).eq.0) call statout(itime)
212 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
215 v_work(ind)=d_t(j,i+nres)
220 write (66,'(80f10.5)')
221 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
225 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
227 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
229 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
232 if (mod(itime,ntwx).eq.0) then
233 write (tytul,'("time",f8.2)') totT
236 call hairpin(.true.,nharp,iharp)
237 call secondary2(.true.)
238 call pdbout(potE,tytul,ipdb)
243 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
244 open(irest2,file=rest2name,status='unknown')
245 write(irest2,*) totT,EK,potE,totE,t_bath
247 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
250 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
261 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
263 & 'MD calculations setup:',t_MDsetup,
264 & 'Energy & gradient evaluation:',t_enegrad,
265 & 'Stochastic MD setup:',t_langsetup,
266 & 'Stochastic MD step setup:',t_sdsetup,
268 write (iout,'(/28(1h=),a25,27(1h=))')
269 & ' End of MD calculation '
271 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
273 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
274 & " time_fricmatmult",time_fricmatmult," time_fsample ",
279 c-------------------------------------------------------------------------------
280 subroutine velverlet_step(itime)
281 c-------------------------------------------------------------------------------
282 c Perform a single velocity Verlet step; the time step can be rescaled if
283 c increments in accelerations exceed the threshold
284 c-------------------------------------------------------------------------------
285 implicit real*8 (a-h,o-z)
289 integer ierror,ierrcode
291 include 'COMMON.SETUP'
292 include 'COMMON.CONTROL'
296 include 'COMMON.LANGEVIN'
298 include 'COMMON.LANGEVIN.lang0'
300 include 'COMMON.CHAIN'
301 include 'COMMON.DERIV'
303 include 'COMMON.LOCAL'
304 include 'COMMON.INTERACT'
305 include 'COMMON.IOUNITS'
306 include 'COMMON.NAMES'
307 include 'COMMON.TIME1'
308 include 'COMMON.MUCA'
309 double precision vcm(3),incr(3)
310 double precision cm(3),L(3)
311 integer ilen,count,rstcount
314 integer maxcount_scale /20/
316 double precision stochforcvec(MAXRES6)
317 common /stochcalc/ stochforcvec
325 else if (lang.eq.2 .or. lang.eq.3) then
327 call stochastic_force(stochforcvec)
330 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
332 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
339 icount_scale=icount_scale+1
340 if (icount_scale.gt.maxcount_scale) then
342 & "ERROR: too many attempts at scaling down the time step. ",
343 & "amax=",amax,"epdrift=",epdrift,
344 & "damax=",damax,"edriftmax=",edriftmax,
348 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
352 c First step of the velocity Verlet algorithm
357 else if (lang.eq.3) then
359 call sd_verlet1_ciccotti
361 else if (lang.eq.1) then
366 c Build the chain from the newly calculated coordinates
368 if (rattle) call rattle1
370 if (large.and. mod(itime,ntwe).eq.0) then
371 write (iout,*) "Cartesian and internal coordinates: step 1"
376 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
377 & (dc(j,i+nres),j=1,3)
379 write (iout,*) "Accelerations"
381 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
382 & (d_a(j,i+nres),j=1,3)
384 write (iout,*) "Velocities, step 1"
386 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
387 & (d_t(j,i+nres),j=1,3)
396 c Calculate energy and forces
398 call etotal(potEcomp)
399 if (large.and. mod(itime,ntwe).eq.0)
400 & call enerprint(potEcomp)
403 t_etotal=t_etotal+MPI_Wtime()-tt0
405 t_etotal=t_etotal+tcpu()-tt0
408 potE=potEcomp(0)-potEcomp(20)
410 c Get the new accelerations
413 t_enegrad=t_enegrad+MPI_Wtime()-tt0
415 t_enegrad=t_enegrad+tcpu()-tt0
417 c Determine maximum acceleration and scale down the timestep if needed
419 amax=amax/(itime_scal+1)**2
420 call predict_edrift(epdrift)
421 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
422 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
424 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
426 itime_scal=itime_scal+ifac_time
427 c fac_time=dmin1(damax/amax,0.5d0)
428 fac_time=0.5d0**ifac_time
429 d_time=d_time*fac_time
430 if (lang.eq.2 .or. lang.eq.3) then
432 c write (iout,*) "Calling sd_verlet_setup: 1"
433 c Rescale the stochastic forces and recalculate or restore
434 c the matrices of tinker integrator
435 if (itime_scal.gt.maxflag_stoch) then
436 if (large) write (iout,'(a,i5,a)')
437 & "Calculate matrices for stochastic step;",
438 & " itime_scal ",itime_scal
440 call sd_verlet_p_setup
442 call sd_verlet_ciccotti_setup
444 write (iout,'(2a,i3,a,i3,1h.)')
445 & "Warning: cannot store matrices for stochastic",
446 & " integration because the index",itime_scal,
447 & " is greater than",maxflag_stoch
448 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
449 & " integration Langevin algorithm for better efficiency."
450 else if (flag_stoch(itime_scal)) then
451 if (large) write (iout,'(a,i5,a,l1)')
452 & "Restore matrices for stochastic step; itime_scal ",
453 & itime_scal," flag ",flag_stoch(itime_scal)
456 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
457 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
458 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
459 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
460 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
461 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
465 if (large) write (iout,'(2a,i5,a,l1)')
466 & "Calculate & store matrices for stochastic step;",
467 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
469 call sd_verlet_p_setup
471 call sd_verlet_ciccotti_setup
473 flag_stoch(ifac_time)=.true.
476 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
477 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
478 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
479 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
480 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
481 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
485 fac_time=1.0d0/dsqrt(fac_time)
487 stochforcvec(i)=fac_time*stochforcvec(i)
490 else if (lang.eq.1) then
491 c Rescale the accelerations due to stochastic forces
492 fac_time=1.0d0/dsqrt(fac_time)
494 d_as_work(i)=d_as_work(i)*fac_time
497 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
498 & "itime",itime," Timestep scaled down to ",
499 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
501 c Second step of the velocity Verlet algorithm
506 else if (lang.eq.3) then
508 call sd_verlet2_ciccotti
510 else if (lang.eq.1) then
515 if (rattle) call rattle2
517 if (d_time.ne.d_time0) then
520 if (lang.eq.2 .or. lang.eq.3) then
521 if (large) write (iout,'(a)')
522 & "Restore original matrices for stochastic step"
523 c write (iout,*) "Calling sd_verlet_setup: 2"
524 c Restore the matrices of tinker integrator if the time step has been restored
527 pfric_mat(i,j)=pfric0_mat(i,j,0)
528 afric_mat(i,j)=afric0_mat(i,j,0)
529 vfric_mat(i,j)=vfric0_mat(i,j,0)
530 prand_mat(i,j)=prand0_mat(i,j,0)
531 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
532 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
541 c Calculate the kinetic and the total energy and the kinetic temperature
546 c write (iout,*) "step",itime," EK",EK," EK1",EK1
548 c Couple the system to Berendsen bath if needed
549 if (tbf .and. lang.eq.0) then
552 kinetic_T=2.0d0/(dimen3*Rb)*EK
553 c Backup the coordinates, velocities, and accelerations
557 d_t_old(j,i)=d_t(j,i)
558 d_a_old(j,i)=d_a(j,i)
562 if (mod(itime,ntwe).eq.0 .and. large) then
563 write (iout,*) "Velocities, step 2"
565 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
566 & (d_t(j,i+nres),j=1,3)
572 c-------------------------------------------------------------------------------
573 subroutine RESPA_step(itime)
574 c-------------------------------------------------------------------------------
575 c Perform a single RESPA step.
576 c-------------------------------------------------------------------------------
577 implicit real*8 (a-h,o-z)
581 integer IERROR,ERRCODE
583 include 'COMMON.SETUP'
584 include 'COMMON.CONTROL'
588 include 'COMMON.LANGEVIN'
590 include 'COMMON.LANGEVIN.lang0'
592 include 'COMMON.CHAIN'
593 include 'COMMON.DERIV'
595 include 'COMMON.LOCAL'
596 include 'COMMON.INTERACT'
597 include 'COMMON.IOUNITS'
598 include 'COMMON.NAMES'
599 include 'COMMON.TIME1'
600 double precision energia_short(0:n_ene),
601 & energia_long(0:n_ene)
602 double precision cm(3),L(3),vcm(3),incr(3)
603 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
604 & d_a_old0(3,0:maxres2)
605 logical PRINT_AMTS_MSG /.false./
606 integer ilen,count,rstcount
609 integer maxcount_scale /10/
611 double precision stochforcvec(MAXRES6)
612 common /stochcalc/ stochforcvec
615 common /cipiszcze/ itt
618 if (large.and. mod(itime,ntwe).eq.0) then
619 write (iout,*) "***************** RESPA itime",itime
620 write (iout,*) "Cartesian and internal coordinates: step 0"
622 call pdbout(0.0d0,"cipiszcze",iout)
624 write (iout,*) "Accelerations from long-range forces"
626 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
627 & (d_a(j,i+nres),j=1,3)
629 write (iout,*) "Velocities, step 0"
631 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
632 & (d_t(j,i+nres),j=1,3)
637 c Perform the initial RESPA step (increment velocities)
638 c write (iout,*) "*********************** RESPA ini"
641 if (mod(itime,ntwe).eq.0 .and. large) then
642 write (iout,*) "Velocities, end"
644 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
645 & (d_t(j,i+nres),j=1,3)
649 c Compute the short-range forces
655 C 7/2/2009 commented out
657 c call etotal_short(energia_short)
660 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
663 d_a(j,i)=d_a_short(j,i)
667 if (large.and. mod(itime,ntwe).eq.0) then
668 write (iout,*) "energia_short",energia_short(0)
669 write (iout,*) "Accelerations from short-range forces"
671 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
672 & (d_a(j,i+nres),j=1,3)
677 t_enegrad=t_enegrad+MPI_Wtime()-tt0
679 t_enegrad=t_enegrad+tcpu()-tt0
684 d_t_old(j,i)=d_t(j,i)
685 d_a_old(j,i)=d_a(j,i)
688 c 6/30/08 A-MTS: attempt at increasing the split number
691 dc_old0(j,i)=dc_old(j,i)
692 d_t_old0(j,i)=d_t_old(j,i)
693 d_a_old0(j,i)=d_a_old(j,i)
696 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
697 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
704 c write (iout,*) "itime",itime," ntime_split",ntime_split
705 c Split the time step
706 d_time=d_time0/ntime_split
707 c Perform the short-range RESPA steps (velocity Verlet increments of
708 c positions and velocities using short-range forces)
709 c write (iout,*) "*********************** RESPA split"
710 do itsplit=1,ntime_split
713 else if (lang.eq.2 .or. lang.eq.3) then
715 call stochastic_force(stochforcvec)
718 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
720 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
725 c First step of the velocity Verlet algorithm
730 else if (lang.eq.3) then
732 call sd_verlet1_ciccotti
734 else if (lang.eq.1) then
739 c Build the chain from the newly calculated coordinates
741 if (rattle) call rattle1
743 if (large.and. mod(itime,ntwe).eq.0) then
744 write (iout,*) "***** ITSPLIT",itsplit
745 write (iout,*) "Cartesian and internal coordinates: step 1"
746 call pdbout(0.0d0,"cipiszcze",iout)
749 write (iout,*) "Velocities, step 1"
751 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
752 & (d_t(j,i+nres),j=1,3)
761 c Calculate energy and forces
763 call etotal_short(energia_short)
764 if (large.and. mod(itime,ntwe).eq.0)
765 & call enerprint(energia_short)
768 t_eshort=t_eshort+MPI_Wtime()-tt0
770 t_eshort=t_eshort+tcpu()-tt0
774 c Get the new accelerations
776 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
779 d_a_short(j,i)=d_a(j,i)
783 if (large.and. mod(itime,ntwe).eq.0) then
784 write (iout,*)"energia_short",energia_short(0)
785 write (iout,*) "Accelerations from short-range forces"
787 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
788 & (d_a(j,i+nres),j=1,3)
793 c Determine maximum acceleration and scale down the timestep if needed
795 amax=amax/ntime_split**2
796 call predict_edrift(epdrift)
797 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
798 & write (iout,*) "amax",amax," damax",damax,
799 & " epdrift",epdrift," epdriftmax",epdriftmax
800 c Exit loop and try with increased split number if the change of
801 c acceleration is too big
802 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
803 if (ntime_split.lt.maxtime_split) then
805 ntime_split=ntime_split*2
808 dc_old(j,i)=dc_old0(j,i)
809 d_t_old(j,i)=d_t_old0(j,i)
810 d_a_old(j,i)=d_a_old0(j,i)
813 if (PRINT_AMTS_MSG) then
814 write (iout,*) "acceleration/energy drift too large",amax,
815 & epdrift," split increased to ",ntime_split," itime",itime,
821 & "Uh-hu. Bumpy landscape. Maximum splitting number",
823 & " already reached!!! Trying to carry on!"
827 t_enegrad=t_enegrad+MPI_Wtime()-tt0
829 t_enegrad=t_enegrad+tcpu()-tt0
831 c Second step of the velocity Verlet algorithm
836 else if (lang.eq.3) then
838 call sd_verlet2_ciccotti
840 else if (lang.eq.1) then
845 if (rattle) call rattle2
846 c Backup the coordinates, velocities, and accelerations
850 d_t_old(j,i)=d_t(j,i)
851 d_a_old(j,i)=d_a(j,i)
858 c Restore the time step
860 c Compute long-range forces
867 call etotal_long(energia_long)
868 if (large.and. mod(itime,ntwe).eq.0)
869 & call enerprint(energia_long)
872 t_elong=t_elong+MPI_Wtime()-tt0
874 t_elong=t_elong+tcpu()-tt0
880 t_enegrad=t_enegrad+MPI_Wtime()-tt0
882 t_enegrad=t_enegrad+tcpu()-tt0
884 c Compute accelerations from long-range forces
886 if (large.and. mod(itime,ntwe).eq.0) then
887 write (iout,*) "energia_long",energia_long(0)
888 write (iout,*) "Cartesian and internal coordinates: step 2"
890 call pdbout(0.0d0,"cipiszcze",iout)
892 write (iout,*) "Accelerations from long-range forces"
894 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
895 & (d_a(j,i+nres),j=1,3)
897 write (iout,*) "Velocities, step 2"
899 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
900 & (d_t(j,i+nres),j=1,3)
904 c Compute the final RESPA step (increment velocities)
905 c write (iout,*) "*********************** RESPA fin"
907 c Compute the complete potential energy
909 potEcomp(i)=energia_short(i)+energia_long(i)
911 potE=potEcomp(0)-potEcomp(20)
912 c potE=energia_short(0)+energia_long(0)
914 c Calculate the kinetic and the total energy and the kinetic temperature
917 c Couple the system to Berendsen bath if needed
918 if (tbf .and. lang.eq.0) then
921 kinetic_T=2.0d0/(dimen3*Rb)*EK
922 c Backup the coordinates, velocities, and accelerations
924 if (mod(itime,ntwe).eq.0 .and. large) then
925 write (iout,*) "Velocities, end"
927 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
928 & (d_t(j,i+nres),j=1,3)
934 c---------------------------------------------------------------------
936 c First and last RESPA step (incrementing velocities using long-range
938 implicit real*8 (a-h,o-z)
940 include 'COMMON.CONTROL'
943 include 'COMMON.CHAIN'
944 include 'COMMON.DERIV'
946 include 'COMMON.LOCAL'
947 include 'COMMON.INTERACT'
948 include 'COMMON.IOUNITS'
949 include 'COMMON.NAMES'
951 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
955 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
959 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
962 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
968 c-----------------------------------------------------------------
970 c Applying velocity Verlet algorithm - step 1 to coordinates
971 implicit real*8 (a-h,o-z)
973 include 'COMMON.CONTROL'
976 include 'COMMON.CHAIN'
977 include 'COMMON.DERIV'
979 include 'COMMON.LOCAL'
980 include 'COMMON.INTERACT'
981 include 'COMMON.IOUNITS'
982 include 'COMMON.NAMES'
983 double precision adt,adt2
986 write (iout,*) "VELVERLET1 START: DC"
988 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
989 & (dc(j,i+nres),j=1,3)
993 adt=d_a_old(j,0)*d_time
995 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
996 d_t_new(j,0)=d_t_old(j,0)+adt2
997 d_t(j,0)=d_t_old(j,0)+adt
1001 adt=d_a_old(j,i)*d_time
1003 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1004 d_t_new(j,i)=d_t_old(j,i)+adt2
1005 d_t(j,i)=d_t_old(j,i)+adt
1009 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1012 adt=d_a_old(j,inres)*d_time
1014 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1015 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1016 d_t(j,inres)=d_t_old(j,inres)+adt
1021 write (iout,*) "VELVERLET1 END: DC"
1023 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1024 & (dc(j,i+nres),j=1,3)
1029 c---------------------------------------------------------------------
1031 c Step 2 of the velocity Verlet algorithm: update velocities
1032 implicit real*8 (a-h,o-z)
1033 include 'DIMENSIONS'
1034 include 'COMMON.CONTROL'
1035 include 'COMMON.VAR'
1037 include 'COMMON.CHAIN'
1038 include 'COMMON.DERIV'
1039 include 'COMMON.GEO'
1040 include 'COMMON.LOCAL'
1041 include 'COMMON.INTERACT'
1042 include 'COMMON.IOUNITS'
1043 include 'COMMON.NAMES'
1045 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1049 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1053 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1056 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1062 c-----------------------------------------------------------------
1063 subroutine sddir_precalc
1064 c Applying velocity Verlet algorithm - step 1 to coordinates
1065 implicit real*8 (a-h,o-z)
1066 include 'DIMENSIONS'
1070 include 'COMMON.CONTROL'
1071 include 'COMMON.VAR'
1074 include 'COMMON.LANGEVIN'
1076 include 'COMMON.LANGEVIN.lang0'
1078 include 'COMMON.CHAIN'
1079 include 'COMMON.DERIV'
1080 include 'COMMON.GEO'
1081 include 'COMMON.LOCAL'
1082 include 'COMMON.INTERACT'
1083 include 'COMMON.IOUNITS'
1084 include 'COMMON.NAMES'
1085 include 'COMMON.TIME1'
1086 double precision stochforcvec(MAXRES6)
1087 common /stochcalc/ stochforcvec
1089 c Compute friction and stochastic forces
1093 time_fric=time_fric+MPI_Wtime()-time00
1095 call stochastic_force(stochforcvec)
1096 time_stoch=time_stoch+MPI_Wtime()-time00
1098 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1099 c forces (d_as_work)
1101 call ginv_mult(fric_work, d_af_work)
1102 call ginv_mult(stochforcvec, d_as_work)
1105 c---------------------------------------------------------------------
1106 subroutine sddir_verlet1
1107 c Applying velocity Verlet algorithm - step 1 to velocities
1108 implicit real*8 (a-h,o-z)
1109 include 'DIMENSIONS'
1110 include 'COMMON.CONTROL'
1111 include 'COMMON.VAR'
1114 include 'COMMON.LANGEVIN'
1116 include 'COMMON.LANGEVIN.lang0'
1118 include 'COMMON.CHAIN'
1119 include 'COMMON.DERIV'
1120 include 'COMMON.GEO'
1121 include 'COMMON.LOCAL'
1122 include 'COMMON.INTERACT'
1123 include 'COMMON.IOUNITS'
1124 include 'COMMON.NAMES'
1125 c Revised 3/31/05 AL: correlation between random contributions to
1126 c position and velocity increments included.
1127 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1128 double precision adt,adt2
1130 c Add the contribution from BOTH friction and stochastic force to the
1131 c coordinates, but ONLY the contribution from the friction forces to velocities
1134 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1135 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1136 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1137 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1138 d_t(j,0)=d_t_old(j,0)+adt
1143 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1144 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1145 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1146 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1147 d_t(j,i)=d_t_old(j,i)+adt
1152 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1155 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1156 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1157 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1158 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1159 d_t(j,inres)=d_t_old(j,inres)+adt
1166 c---------------------------------------------------------------------
1167 subroutine sddir_verlet2
1168 c Calculating the adjusted velocities for accelerations
1169 implicit real*8 (a-h,o-z)
1170 include 'DIMENSIONS'
1171 include 'COMMON.CONTROL'
1172 include 'COMMON.VAR'
1175 include 'COMMON.LANGEVIN'
1177 include 'COMMON.LANGEVIN.lang0'
1179 include 'COMMON.CHAIN'
1180 include 'COMMON.DERIV'
1181 include 'COMMON.GEO'
1182 include 'COMMON.LOCAL'
1183 include 'COMMON.INTERACT'
1184 include 'COMMON.IOUNITS'
1185 include 'COMMON.NAMES'
1186 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1187 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1188 c Revised 3/31/05 AL: correlation between random contributions to
1189 c position and velocity increments included.
1190 c The correlation coefficients are calculated at low-friction limit.
1191 c Also, friction forces are now not calculated with new velocities.
1193 c call friction_force
1194 call stochastic_force(stochforcvec)
1196 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1197 c forces (d_as_work)
1199 call ginv_mult(stochforcvec, d_as_work1)
1205 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1206 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1211 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1212 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1217 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1220 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1221 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1222 & +cos60*d_as_work1(ind+j))*d_time
1229 c---------------------------------------------------------------------
1230 subroutine max_accel
1232 c Find the maximum difference in the accelerations of the the sites
1233 c at the beginning and the end of the time step.
1235 implicit real*8 (a-h,o-z)
1236 include 'DIMENSIONS'
1237 include 'COMMON.CONTROL'
1238 include 'COMMON.VAR'
1240 include 'COMMON.CHAIN'
1241 include 'COMMON.DERIV'
1242 include 'COMMON.GEO'
1243 include 'COMMON.LOCAL'
1244 include 'COMMON.INTERACT'
1245 include 'COMMON.IOUNITS'
1246 double precision aux(3),accel(3),accel_old(3),dacc
1248 c aux(j)=d_a(j,0)-d_a_old(j,0)
1249 accel_old(j)=d_a_old(j,0)
1256 c 7/3/08 changed to asymmetric difference
1258 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1259 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1260 accel(j)=accel(j)+0.5d0*d_a(j,i)
1261 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1262 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1263 dacc=dabs(accel(j)-accel_old(j))
1264 c write (iout,*) i,dacc
1265 if (dacc.gt.amax) amax=dacc
1273 accel_old(j)=d_a_old(j,0)
1278 accel_old(j)=accel_old(j)+d_a_old(j,1)
1279 accel(j)=accel(j)+d_a(j,1)
1283 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1285 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1286 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1287 accel(j)=accel(j)+d_a(j,i+nres)
1291 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1292 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1293 dacc=dabs(accel(j)-accel_old(j))
1294 c write (iout,*) "side-chain",i,dacc
1295 if (dacc.gt.amax) amax=dacc
1299 accel_old(j)=accel_old(j)+d_a_old(j,i)
1300 accel(j)=accel(j)+d_a(j,i)
1301 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1306 c---------------------------------------------------------------------
1307 subroutine predict_edrift(epdrift)
1309 c Predict the drift of the potential energy
1311 implicit real*8 (a-h,o-z)
1312 include 'DIMENSIONS'
1313 include 'COMMON.CONTROL'
1314 include 'COMMON.VAR'
1316 include 'COMMON.CHAIN'
1317 include 'COMMON.DERIV'
1318 include 'COMMON.GEO'
1319 include 'COMMON.LOCAL'
1320 include 'COMMON.INTERACT'
1321 include 'COMMON.IOUNITS'
1322 include 'COMMON.MUCA'
1323 double precision epdrift,epdriftij
1324 c Drift of the potential energy
1330 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1331 if (lmuca) epdriftij=epdriftij*factor
1332 c write (iout,*) "back",i,j,epdriftij
1333 if (epdriftij.gt.epdrift) epdrift=epdriftij
1337 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1340 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1341 if (lmuca) epdriftij=epdriftij*factor
1342 c write (iout,*) "side",i,j,epdriftij
1343 if (epdriftij.gt.epdrift) epdrift=epdriftij
1347 epdrift=0.5d0*epdrift*d_time*d_time
1348 c write (iout,*) "epdrift",epdrift
1351 c-----------------------------------------------------------------------
1352 subroutine verlet_bath
1354 c Coupling to the thermostat by using the Berendsen algorithm
1356 implicit real*8 (a-h,o-z)
1357 include 'DIMENSIONS'
1358 include 'COMMON.CONTROL'
1359 include 'COMMON.VAR'
1361 include 'COMMON.CHAIN'
1362 include 'COMMON.DERIV'
1363 include 'COMMON.GEO'
1364 include 'COMMON.LOCAL'
1365 include 'COMMON.INTERACT'
1366 include 'COMMON.IOUNITS'
1367 include 'COMMON.NAMES'
1368 double precision T_half,fact
1370 T_half=2.0d0/(dimen3*Rb)*EK
1371 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1372 c write(iout,*) "T_half", T_half
1373 c write(iout,*) "EK", EK
1374 c write(iout,*) "fact", fact
1376 d_t(j,0)=fact*d_t(j,0)
1380 d_t(j,i)=fact*d_t(j,i)
1384 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1387 d_t(j,inres)=fact*d_t(j,inres)
1393 c---------------------------------------------------------
1395 c Set up the initial conditions of a MD simulation
1396 implicit real*8 (a-h,o-z)
1397 include 'DIMENSIONS'
1401 integer IERROR,ERRCODE
1403 include 'COMMON.SETUP'
1404 include 'COMMON.CONTROL'
1405 include 'COMMON.VAR'
1408 include 'COMMON.LANGEVIN'
1410 include 'COMMON.LANGEVIN.lang0'
1412 include 'COMMON.CHAIN'
1413 include 'COMMON.DERIV'
1414 include 'COMMON.GEO'
1415 include 'COMMON.LOCAL'
1416 include 'COMMON.INTERACT'
1417 include 'COMMON.IOUNITS'
1418 include 'COMMON.NAMES'
1419 include 'COMMON.REMD'
1420 real*8 energia_long(0:n_ene),
1421 & energia_short(0:n_ene),vcm(3),incr(3)
1422 double precision cm(3),L(3),xv,sigv,lowb,highb
1423 double precision varia(maxvar)
1431 c write(iout,*) "d_time", d_time
1432 c Compute the standard deviations of stochastic forces for Langevin dynamics
1433 c if the friction coefficients do not depend on surface area
1434 if (lang.gt.0 .and. .not.surfarea) then
1436 stdforcp(i)=stdfp*dsqrt(gamp)
1439 stdforcsc(i)=stdfsc(iabs(itype(i)))
1440 & *dsqrt(gamsc(iabs(itype(i))))
1443 c Open the pdb file for snapshotshots
1446 if (ilen(tmpdir).gt.0)
1447 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1448 & liczba(:ilen(liczba))//".pdb")
1450 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1454 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1455 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1456 & liczba(:ilen(liczba))//".x")
1457 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1460 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1461 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1462 & liczba(:ilen(liczba))//".cx")
1463 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1469 if (ilen(tmpdir).gt.0)
1470 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1471 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1473 if (ilen(tmpdir).gt.0)
1474 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1475 cartname=prefix(:ilen(prefix))//"_MD.cx"
1479 write (qstr,'(256(1h ))')
1482 iq = qinfrag(i,iset)*10
1483 iw = wfrag(i,iset)/100
1485 if(me.eq.king.or..not.out1file)
1486 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1487 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1492 iq = qinpair(i,iset)*10
1493 iw = wpair(i,iset)/100
1495 if(me.eq.king.or..not.out1file)
1496 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1497 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1501 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1503 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1505 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1507 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1511 if (restart1file) then
1513 & inquire(file=mremd_rst_name,exist=file_exist)
1514 write (*,*) me," Before broadcast: file_exist",file_exist
1515 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1517 write (*,*) me," After broadcast: file_exist",file_exist
1518 c inquire(file=mremd_rst_name,exist=file_exist)
1519 if(me.eq.king.or..not.out1file)
1520 & write(iout,*) "Initial state read by master and distributed"
1522 if (ilen(tmpdir).gt.0)
1523 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1524 & //liczba(:ilen(liczba))//'.rst')
1525 inquire(file=rest2name,exist=file_exist)
1528 if(.not.restart1file) then
1529 if(me.eq.king.or..not.out1file)
1530 & write(iout,*) "Initial state will be read from file ",
1531 & rest2name(:ilen(rest2name))
1534 call rescale_weights(t_bath)
1536 if(me.eq.king.or..not.out1file)then
1537 if (restart1file) then
1538 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1541 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1544 write(iout,*) "Initial velocities randomly generated"
1550 c Generate initial velocities
1551 if(me.eq.king.or..not.out1file)
1552 & write(iout,*) "Initial velocities randomly generated"
1556 c rest2name = prefix(:ilen(prefix))//'.rst'
1557 if(me.eq.king.or..not.out1file)then
1558 write (iout,*) "Initial velocities"
1560 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1561 & (d_t(j,i+nres),j=1,3)
1563 c Zeroing the total angular momentum of the system
1564 write(iout,*) "Calling the zero-angular
1565 & momentum subroutine"
1568 c Getting the potential energy and forces and velocities and accelerations
1570 c write (iout,*) "velocity of the center of the mass:"
1571 c write (iout,*) (vcm(j),j=1,3)
1573 d_t(j,0)=d_t(j,0)-vcm(j)
1575 c Removing the velocity of the center of mass
1577 if(me.eq.king.or..not.out1file)then
1578 write (iout,*) "vcm right after adjustment:"
1579 write (iout,*) (vcm(j),j=1,3)
1583 if(iranconf.ne.0) then
1585 print *, 'Calling OVERLAP_SC'
1586 call overlap_sc(fail)
1590 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1591 print *,'SC_move',nft_sc,etot
1592 if(me.eq.king.or..not.out1file)
1593 & write(iout,*) 'SC_move',nft_sc,etot
1597 print *, 'Calling MINIM_DC'
1598 call minim_dc(etot,iretcode,nfun)
1600 call geom_to_var(nvar,varia)
1601 print *,'Calling MINIMIZE.'
1602 call minimize(etot,varia,iretcode,nfun)
1603 call var_to_geom(nvar,varia)
1605 if(me.eq.king.or..not.out1file)
1606 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1609 call chainbuild_cart
1614 kinetic_T=2.0d0/(dimen3*Rb)*EK
1615 if(me.eq.king.or..not.out1file)then
1625 call etotal(potEcomp)
1626 if (large) call enerprint(potEcomp)
1629 t_etotal=t_etotal+MPI_Wtime()-tt0
1631 t_etotal=t_etotal+tcpu()-tt0
1638 if (amax*d_time .gt. dvmax) then
1639 d_time=d_time*dvmax/amax
1640 if(me.eq.king.or..not.out1file) write (iout,*)
1641 & "Time step reduced to",d_time,
1642 & " because of too large initial acceleration."
1644 if(me.eq.king.or..not.out1file)then
1645 write(iout,*) "Potential energy and its components"
1646 call enerprint(potEcomp)
1647 c write(iout,*) (potEcomp(i),i=0,n_ene)
1649 potE=potEcomp(0)-potEcomp(20)
1652 if (ntwe.ne.0) call statout(itime)
1653 if(me.eq.king.or..not.out1file)
1654 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1655 & " Kinetic energy",EK," Potential energy",potE,
1656 & " Total energy",totE," Maximum acceleration ",
1659 write (iout,*) "Initial coordinates"
1661 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1662 & (c(j,i+nres),j=1,3)
1664 write (iout,*) "Initial dC"
1666 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1667 & (dc(j,i+nres),j=1,3)
1669 write (iout,*) "Initial velocities"
1670 write (iout,"(13x,' backbone ',23x,' side chain')")
1672 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1673 & (d_t(j,i+nres),j=1,3)
1675 write (iout,*) "Initial accelerations"
1677 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1678 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1679 & (d_a(j,i+nres),j=1,3)
1685 d_t_old(j,i)=d_t(j,i)
1686 d_a_old(j,i)=d_a(j,i)
1688 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1697 call etotal_short(energia_short)
1698 if (large) call enerprint(potEcomp)
1701 t_eshort=t_eshort+MPI_Wtime()-tt0
1703 t_eshort=t_eshort+tcpu()-tt0
1708 if(.not.out1file .and. large) then
1709 write (iout,*) "energia_long",energia_long(0),
1710 & " energia_short",energia_short(0),
1711 & " total",energia_long(0)+energia_short(0)
1712 write (iout,*) "Initial fast-force accelerations"
1714 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1715 & (d_a(j,i+nres),j=1,3)
1718 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1721 d_a_short(j,i)=d_a(j,i)
1730 call etotal_long(energia_long)
1731 if (large) call enerprint(potEcomp)
1734 t_elong=t_elong+MPI_Wtime()-tt0
1736 t_elong=t_elong+tcpu()-tt0
1741 if(.not.out1file .and. large) then
1742 write (iout,*) "energia_long",energia_long(0)
1743 write (iout,*) "Initial slow-force accelerations"
1745 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1746 & (d_a(j,i+nres),j=1,3)
1750 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1752 t_enegrad=t_enegrad+tcpu()-tt0
1757 c-----------------------------------------------------------
1758 subroutine random_vel
1759 implicit real*8 (a-h,o-z)
1760 include 'DIMENSIONS'
1761 include 'COMMON.CONTROL'
1762 include 'COMMON.VAR'
1765 include 'COMMON.LANGEVIN'
1767 include 'COMMON.LANGEVIN.lang0'
1769 include 'COMMON.CHAIN'
1770 include 'COMMON.DERIV'
1771 include 'COMMON.GEO'
1772 include 'COMMON.LOCAL'
1773 include 'COMMON.INTERACT'
1774 include 'COMMON.IOUNITS'
1775 include 'COMMON.NAMES'
1776 include 'COMMON.TIME1'
1777 double precision xv,sigv,lowb,highb
1778 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1779 c First generate velocities in the eigenspace of the G matrix
1780 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1787 sigv=dsqrt((Rb*t_bath)/geigen(i))
1790 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1791 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1792 c & " d_t_work_new",d_t_work_new(ii)
1801 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1804 c write (iout,*) "Ek from eigenvectors",Ek1
1806 c Transform velocities to UNRES coordinate space
1812 d_t_work(ind)=d_t_work(ind)
1813 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1815 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1819 c Transfer to the d_t vector
1821 d_t(j,0)=d_t_work(j)
1827 d_t(j,i)=d_t_work(ind)
1831 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1834 d_t(j,i+nres)=d_t_work(ind)
1839 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1840 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1845 c-----------------------------------------------------------
1846 subroutine sd_verlet_p_setup
1847 c Sets up the parameters of stochastic Verlet algorithm
1848 implicit real*8 (a-h,o-z)
1849 include 'DIMENSIONS'
1853 include 'COMMON.CONTROL'
1854 include 'COMMON.VAR'
1857 include 'COMMON.LANGEVIN'
1859 include 'COMMON.LANGEVIN.lang0'
1861 include 'COMMON.CHAIN'
1862 include 'COMMON.DERIV'
1863 include 'COMMON.GEO'
1864 include 'COMMON.LOCAL'
1865 include 'COMMON.INTERACT'
1866 include 'COMMON.IOUNITS'
1867 include 'COMMON.NAMES'
1868 include 'COMMON.TIME1'
1869 double precision emgdt(MAXRES6),
1870 & pterm,vterm,rho,rhoc,vsig,
1871 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1872 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1873 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1874 logical lprn /.false./
1875 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1876 double precision ktm
1883 c AL 8/17/04 Code adapted from tinker
1885 c Get the frictional and random terms for stochastic dynamics in the
1886 c eigenspace of mass-scaled UNRES friction matrix
1889 gdt = fricgam(i) * d_time
1891 c Stochastic dynamics reduces to simple MD for zero friction
1893 if (gdt .le. zero) then
1894 pfric_vec(i) = 1.0d0
1895 vfric_vec(i) = d_time
1896 afric_vec(i) = 0.5d0 * d_time * d_time
1897 prand_vec(i) = 0.0d0
1898 vrand_vec1(i) = 0.0d0
1899 vrand_vec2(i) = 0.0d0
1901 c Analytical expressions when friction coefficient is large
1904 if (gdt .ge. gdt_radius) then
1907 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1908 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1909 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1910 vterm = 1.0d0 - egdt**2
1911 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1913 c Use series expansions when friction coefficient is small
1924 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1925 & - gdt5/120.0d0 + gdt6/720.0d0
1926 & - gdt7/5040.0d0 + gdt8/40320.0d0
1927 & - gdt9/362880.0d0) / fricgam(i)**2
1928 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1929 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1930 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1931 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1932 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1933 & + 127.0d0*gdt9/90720.0d0
1934 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1935 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1936 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1937 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1938 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1939 & - 17.0d0*gdt2/1280.0d0
1940 & + 17.0d0*gdt3/6144.0d0
1941 & + 40967.0d0*gdt4/34406400.0d0
1942 & - 57203.0d0*gdt5/275251200.0d0
1943 & - 1429487.0d0*gdt6/13212057600.0d0)
1946 c Compute the scaling factors of random terms for the nonzero friction case
1948 ktm = 0.5d0*d_time/fricgam(i)
1949 psig = dsqrt(ktm*pterm) / fricgam(i)
1950 vsig = dsqrt(ktm*vterm)
1951 rhoc = dsqrt(1.0d0 - rho*rho)
1953 vrand_vec1(i) = vsig * rho
1954 vrand_vec2(i) = vsig * rhoc
1959 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1962 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1963 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1967 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1970 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1971 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1972 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1973 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1974 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1975 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1978 t_sdsetup=t_sdsetup+MPI_Wtime()
1980 t_sdsetup=t_sdsetup+tcpu()-tt0
1984 c-------------------------------------------------------------
1985 subroutine eigtransf1(n,ndim,ab,d,c)
1988 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
1994 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2000 c-------------------------------------------------------------
2001 subroutine eigtransf(n,ndim,a,b,d,c)
2004 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2010 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2016 c-------------------------------------------------------------
2017 subroutine sd_verlet1
2018 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2019 implicit real*8 (a-h,o-z)
2020 include 'DIMENSIONS'
2021 include 'COMMON.CONTROL'
2022 include 'COMMON.VAR'
2025 include 'COMMON.LANGEVIN'
2027 include 'COMMON.LANGEVIN.lang0'
2029 include 'COMMON.CHAIN'
2030 include 'COMMON.DERIV'
2031 include 'COMMON.GEO'
2032 include 'COMMON.LOCAL'
2033 include 'COMMON.INTERACT'
2034 include 'COMMON.IOUNITS'
2035 include 'COMMON.NAMES'
2036 double precision stochforcvec(MAXRES6)
2037 common /stochcalc/ stochforcvec
2038 logical lprn /.false./
2040 c write (iout,*) "dc_old"
2042 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2043 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2046 dc_work(j)=dc_old(j,0)
2047 d_t_work(j)=d_t_old(j,0)
2048 d_a_work(j)=d_a_old(j,0)
2053 dc_work(ind+j)=dc_old(j,i)
2054 d_t_work(ind+j)=d_t_old(j,i)
2055 d_a_work(ind+j)=d_a_old(j,i)
2060 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2062 dc_work(ind+j)=dc_old(j,i+nres)
2063 d_t_work(ind+j)=d_t_old(j,i+nres)
2064 d_a_work(ind+j)=d_a_old(j,i+nres)
2072 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2076 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2077 & vfric_mat(i,j),afric_mat(i,j),
2078 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2086 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2087 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2088 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2089 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2091 d_t_work_new(i)=ddt1+0.5d0*ddt2
2092 d_t_work(i)=ddt1+ddt2
2097 d_t(j,0)=d_t_work(j)
2102 dc(j,i)=dc_work(ind+j)
2103 d_t(j,i)=d_t_work(ind+j)
2108 if (itype(i).ne.10) then
2111 dc(j,inres)=dc_work(ind+j)
2112 d_t(j,inres)=d_t_work(ind+j)
2119 c--------------------------------------------------------------------------
2120 subroutine sd_verlet2
2121 c Calculating the adjusted velocities for accelerations
2122 implicit real*8 (a-h,o-z)
2123 include 'DIMENSIONS'
2124 include 'COMMON.CONTROL'
2125 include 'COMMON.VAR'
2128 include 'COMMON.LANGEVIN'
2130 include 'COMMON.LANGEVIN.lang0'
2132 include 'COMMON.CHAIN'
2133 include 'COMMON.DERIV'
2134 include 'COMMON.GEO'
2135 include 'COMMON.LOCAL'
2136 include 'COMMON.INTERACT'
2137 include 'COMMON.IOUNITS'
2138 include 'COMMON.NAMES'
2139 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2140 common /stochcalc/ stochforcvec
2142 c Compute the stochastic forces which contribute to velocity change
2144 call stochastic_force(stochforcvecV)
2151 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2152 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2153 & vrand_mat2(i,j)*stochforcvecV(j)
2155 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2159 d_t(j,0)=d_t_work(j)
2164 d_t(j,i)=d_t_work(ind+j)
2169 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2172 d_t(j,inres)=d_t_work(ind+j)
2179 c-----------------------------------------------------------
2180 subroutine sd_verlet_ciccotti_setup
2181 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2183 implicit real*8 (a-h,o-z)
2184 include 'DIMENSIONS'
2188 include 'COMMON.CONTROL'
2189 include 'COMMON.VAR'
2192 include 'COMMON.LANGEVIN'
2194 include 'COMMON.LANGEVIN.lang0'
2196 include 'COMMON.CHAIN'
2197 include 'COMMON.DERIV'
2198 include 'COMMON.GEO'
2199 include 'COMMON.LOCAL'
2200 include 'COMMON.INTERACT'
2201 include 'COMMON.IOUNITS'
2202 include 'COMMON.NAMES'
2203 include 'COMMON.TIME1'
2204 double precision emgdt(MAXRES6),
2205 & pterm,vterm,rho,rhoc,vsig,
2206 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2207 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2208 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2209 logical lprn /.false./
2210 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2211 double precision ktm
2218 c AL 8/17/04 Code adapted from tinker
2220 c Get the frictional and random terms for stochastic dynamics in the
2221 c eigenspace of mass-scaled UNRES friction matrix
2224 write (iout,*) "i",i," fricgam",fricgam(i)
2225 gdt = fricgam(i) * d_time
2227 c Stochastic dynamics reduces to simple MD for zero friction
2229 if (gdt .le. zero) then
2230 pfric_vec(i) = 1.0d0
2231 vfric_vec(i) = d_time
2232 afric_vec(i) = 0.5d0*d_time*d_time
2233 prand_vec(i) = afric_vec(i)
2234 vrand_vec2(i) = vfric_vec(i)
2236 c Analytical expressions when friction coefficient is large
2241 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2242 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2243 prand_vec(i) = afric_vec(i)
2244 vrand_vec2(i) = vfric_vec(i)
2246 c Compute the scaling factors of random terms for the nonzero friction case
2248 c ktm = 0.5d0*d_time/fricgam(i)
2249 c psig = dsqrt(ktm*pterm) / fricgam(i)
2250 c vsig = dsqrt(ktm*vterm)
2251 c prand_vec(i) = psig*afric_vec(i)
2252 c vrand_vec2(i) = vsig*vfric_vec(i)
2257 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2260 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2261 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2265 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2267 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2268 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2269 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2270 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2271 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2273 t_sdsetup=t_sdsetup+MPI_Wtime()
2275 t_sdsetup=t_sdsetup+tcpu()-tt0
2279 c-------------------------------------------------------------
2280 subroutine sd_verlet1_ciccotti
2281 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2282 implicit real*8 (a-h,o-z)
2283 include 'DIMENSIONS'
2287 include 'COMMON.CONTROL'
2288 include 'COMMON.VAR'
2291 include 'COMMON.LANGEVIN'
2293 include 'COMMON.LANGEVIN.lang0'
2295 include 'COMMON.CHAIN'
2296 include 'COMMON.DERIV'
2297 include 'COMMON.GEO'
2298 include 'COMMON.LOCAL'
2299 include 'COMMON.INTERACT'
2300 include 'COMMON.IOUNITS'
2301 include 'COMMON.NAMES'
2302 double precision stochforcvec(MAXRES6)
2303 common /stochcalc/ stochforcvec
2304 logical lprn /.false./
2306 c write (iout,*) "dc_old"
2308 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2309 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2312 dc_work(j)=dc_old(j,0)
2313 d_t_work(j)=d_t_old(j,0)
2314 d_a_work(j)=d_a_old(j,0)
2319 dc_work(ind+j)=dc_old(j,i)
2320 d_t_work(ind+j)=d_t_old(j,i)
2321 d_a_work(ind+j)=d_a_old(j,i)
2326 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2328 dc_work(ind+j)=dc_old(j,i+nres)
2329 d_t_work(ind+j)=d_t_old(j,i+nres)
2330 d_a_work(ind+j)=d_a_old(j,i+nres)
2339 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2343 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2344 & vfric_mat(i,j),afric_mat(i,j),
2345 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2353 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2354 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2355 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2356 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2358 d_t_work_new(i)=ddt1+0.5d0*ddt2
2359 d_t_work(i)=ddt1+ddt2
2364 d_t(j,0)=d_t_work(j)
2369 dc(j,i)=dc_work(ind+j)
2370 d_t(j,i)=d_t_work(ind+j)
2375 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2378 dc(j,inres)=dc_work(ind+j)
2379 d_t(j,inres)=d_t_work(ind+j)
2386 c--------------------------------------------------------------------------
2387 subroutine sd_verlet2_ciccotti
2388 c Calculating the adjusted velocities for accelerations
2389 implicit real*8 (a-h,o-z)
2390 include 'DIMENSIONS'
2391 include 'COMMON.CONTROL'
2392 include 'COMMON.VAR'
2395 include 'COMMON.LANGEVIN'
2397 include 'COMMON.LANGEVIN.lang0'
2399 include 'COMMON.CHAIN'
2400 include 'COMMON.DERIV'
2401 include 'COMMON.GEO'
2402 include 'COMMON.LOCAL'
2403 include 'COMMON.INTERACT'
2404 include 'COMMON.IOUNITS'
2405 include 'COMMON.NAMES'
2406 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2407 common /stochcalc/ stochforcvec
2409 c Compute the stochastic forces which contribute to velocity change
2411 call stochastic_force(stochforcvecV)
2418 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2419 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2420 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2422 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2426 d_t(j,0)=d_t_work(j)
2431 d_t(j,i)=d_t_work(ind+j)
2436 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2439 d_t(j,inres)=d_t_work(ind+j)