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
234 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
1098 time_fric=time_fric+MPI_Wtime()-time00
1101 time_fric=time_fric+tcpu()-time00
1104 call stochastic_force(stochforcvec)
1106 time_stoch=time_stoch+MPI_Wtime()-time00
1108 time_stoch=time_stoch+tcpu()-time00
1111 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1112 c forces (d_as_work)
1114 call ginv_mult(fric_work, d_af_work)
1115 call ginv_mult(stochforcvec, d_as_work)
1118 c---------------------------------------------------------------------
1119 subroutine sddir_verlet1
1120 c Applying velocity Verlet algorithm - step 1 to velocities
1121 implicit real*8 (a-h,o-z)
1122 include 'DIMENSIONS'
1123 include 'COMMON.CONTROL'
1124 include 'COMMON.VAR'
1127 include 'COMMON.LANGEVIN'
1129 include 'COMMON.LANGEVIN.lang0'
1131 include 'COMMON.CHAIN'
1132 include 'COMMON.DERIV'
1133 include 'COMMON.GEO'
1134 include 'COMMON.LOCAL'
1135 include 'COMMON.INTERACT'
1136 include 'COMMON.IOUNITS'
1137 include 'COMMON.NAMES'
1138 c Revised 3/31/05 AL: correlation between random contributions to
1139 c position and velocity increments included.
1140 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1141 double precision adt,adt2
1143 c Add the contribution from BOTH friction and stochastic force to the
1144 c coordinates, but ONLY the contribution from the friction forces to velocities
1147 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1148 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1149 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1150 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1151 d_t(j,0)=d_t_old(j,0)+adt
1156 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1157 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1158 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1159 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1160 d_t(j,i)=d_t_old(j,i)+adt
1165 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1168 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1169 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1170 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1171 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1172 d_t(j,inres)=d_t_old(j,inres)+adt
1179 c---------------------------------------------------------------------
1180 subroutine sddir_verlet2
1181 c Calculating the adjusted velocities for accelerations
1182 implicit real*8 (a-h,o-z)
1183 include 'DIMENSIONS'
1184 include 'COMMON.CONTROL'
1185 include 'COMMON.VAR'
1188 include 'COMMON.LANGEVIN'
1190 include 'COMMON.LANGEVIN.lang0'
1192 include 'COMMON.CHAIN'
1193 include 'COMMON.DERIV'
1194 include 'COMMON.GEO'
1195 include 'COMMON.LOCAL'
1196 include 'COMMON.INTERACT'
1197 include 'COMMON.IOUNITS'
1198 include 'COMMON.NAMES'
1199 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1200 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1201 c Revised 3/31/05 AL: correlation between random contributions to
1202 c position and velocity increments included.
1203 c The correlation coefficients are calculated at low-friction limit.
1204 c Also, friction forces are now not calculated with new velocities.
1206 c call friction_force
1207 call stochastic_force(stochforcvec)
1209 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1210 c forces (d_as_work)
1212 call ginv_mult(stochforcvec, d_as_work1)
1218 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1219 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1224 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1225 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1230 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1233 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1234 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1235 & +cos60*d_as_work1(ind+j))*d_time
1242 c---------------------------------------------------------------------
1243 subroutine max_accel
1245 c Find the maximum difference in the accelerations of the the sites
1246 c at the beginning and the end of the time step.
1248 implicit real*8 (a-h,o-z)
1249 include 'DIMENSIONS'
1250 include 'COMMON.CONTROL'
1251 include 'COMMON.VAR'
1253 include 'COMMON.CHAIN'
1254 include 'COMMON.DERIV'
1255 include 'COMMON.GEO'
1256 include 'COMMON.LOCAL'
1257 include 'COMMON.INTERACT'
1258 include 'COMMON.IOUNITS'
1259 double precision aux(3),accel(3),accel_old(3),dacc
1261 c aux(j)=d_a(j,0)-d_a_old(j,0)
1262 accel_old(j)=d_a_old(j,0)
1269 c 7/3/08 changed to asymmetric difference
1271 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1272 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1273 accel(j)=accel(j)+0.5d0*d_a(j,i)
1274 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1275 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1276 dacc=dabs(accel(j)-accel_old(j))
1277 c write (iout,*) i,dacc
1278 if (dacc.gt.amax) amax=dacc
1286 accel_old(j)=d_a_old(j,0)
1291 accel_old(j)=accel_old(j)+d_a_old(j,1)
1292 accel(j)=accel(j)+d_a(j,1)
1296 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1298 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1299 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1300 accel(j)=accel(j)+d_a(j,i+nres)
1304 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1305 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1306 dacc=dabs(accel(j)-accel_old(j))
1307 c write (iout,*) "side-chain",i,dacc
1308 if (dacc.gt.amax) amax=dacc
1312 accel_old(j)=accel_old(j)+d_a_old(j,i)
1313 accel(j)=accel(j)+d_a(j,i)
1314 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1319 c---------------------------------------------------------------------
1320 subroutine predict_edrift(epdrift)
1322 c Predict the drift of the potential energy
1324 implicit real*8 (a-h,o-z)
1325 include 'DIMENSIONS'
1326 include 'COMMON.CONTROL'
1327 include 'COMMON.VAR'
1329 include 'COMMON.CHAIN'
1330 include 'COMMON.DERIV'
1331 include 'COMMON.GEO'
1332 include 'COMMON.LOCAL'
1333 include 'COMMON.INTERACT'
1334 include 'COMMON.IOUNITS'
1335 include 'COMMON.MUCA'
1336 double precision epdrift,epdriftij
1337 c Drift of the potential energy
1343 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1344 if (lmuca) epdriftij=epdriftij*factor
1345 c write (iout,*) "back",i,j,epdriftij
1346 if (epdriftij.gt.epdrift) epdrift=epdriftij
1350 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1353 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1354 if (lmuca) epdriftij=epdriftij*factor
1355 c write (iout,*) "side",i,j,epdriftij
1356 if (epdriftij.gt.epdrift) epdrift=epdriftij
1360 epdrift=0.5d0*epdrift*d_time*d_time
1361 c write (iout,*) "epdrift",epdrift
1364 c-----------------------------------------------------------------------
1365 subroutine verlet_bath
1367 c Coupling to the thermostat by using the Berendsen algorithm
1369 implicit real*8 (a-h,o-z)
1370 include 'DIMENSIONS'
1371 include 'COMMON.CONTROL'
1372 include 'COMMON.VAR'
1374 include 'COMMON.CHAIN'
1375 include 'COMMON.DERIV'
1376 include 'COMMON.GEO'
1377 include 'COMMON.LOCAL'
1378 include 'COMMON.INTERACT'
1379 include 'COMMON.IOUNITS'
1380 include 'COMMON.NAMES'
1381 double precision T_half,fact
1383 T_half=2.0d0/(dimen3*Rb)*EK
1384 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1385 c write(iout,*) "T_half", T_half
1386 c write(iout,*) "EK", EK
1387 c write(iout,*) "fact", fact
1389 d_t(j,0)=fact*d_t(j,0)
1393 d_t(j,i)=fact*d_t(j,i)
1397 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1400 d_t(j,inres)=fact*d_t(j,inres)
1406 c---------------------------------------------------------
1408 c Set up the initial conditions of a MD simulation
1409 implicit real*8 (a-h,o-z)
1410 include 'DIMENSIONS'
1414 integer IERROR,ERRCODE
1416 include 'COMMON.SETUP'
1417 include 'COMMON.CONTROL'
1418 include 'COMMON.VAR'
1421 include 'COMMON.LANGEVIN'
1423 include 'COMMON.LANGEVIN.lang0'
1425 include 'COMMON.CHAIN'
1426 include 'COMMON.DERIV'
1427 include 'COMMON.GEO'
1428 include 'COMMON.LOCAL'
1429 include 'COMMON.INTERACT'
1430 include 'COMMON.IOUNITS'
1431 include 'COMMON.NAMES'
1432 include 'COMMON.REMD'
1433 real*8 energia_long(0:n_ene),
1434 & energia_short(0:n_ene),vcm(3),incr(3)
1435 double precision cm(3),L(3),xv,sigv,lowb,highb
1436 double precision varia(maxvar)
1444 c write(iout,*) "d_time", d_time
1445 c Compute the standard deviations of stochastic forces for Langevin dynamics
1446 c if the friction coefficients do not depend on surface area
1447 if (lang.gt.0 .and. .not.surfarea) then
1449 stdforcp(i)=stdfp*dsqrt(gamp)
1452 stdforcsc(i)=stdfsc(iabs(itype(i)))
1453 & *dsqrt(gamsc(iabs(itype(i))))
1456 c Open the pdb file for snapshotshots
1459 if (ilen(tmpdir).gt.0)
1460 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1461 & liczba(:ilen(liczba))//".pdb")
1463 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1467 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1468 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1469 & liczba(:ilen(liczba))//".x")
1470 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1473 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1474 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1475 & liczba(:ilen(liczba))//".cx")
1476 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1482 if (ilen(tmpdir).gt.0)
1483 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1484 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1486 if (ilen(tmpdir).gt.0)
1487 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1488 cartname=prefix(:ilen(prefix))//"_MD.cx"
1492 write (qstr,'(256(1h ))')
1495 iq = qinfrag(i,iset)*10
1496 iw = wfrag(i,iset)/100
1498 if(me.eq.king.or..not.out1file)
1499 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1500 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1505 iq = qinpair(i,iset)*10
1506 iw = wpair(i,iset)/100
1508 if(me.eq.king.or..not.out1file)
1509 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1510 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1514 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1516 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1518 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1520 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1524 if (restart1file) then
1526 & inquire(file=mremd_rst_name,exist=file_exist)
1528 write (*,*) me," Before broadcast: file_exist",file_exist
1529 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1531 write (*,*) me," After broadcast: file_exist",file_exist
1532 c inquire(file=mremd_rst_name,exist=file_exist)
1534 if(me.eq.king.or..not.out1file)
1535 & write(iout,*) "Initial state read by master and distributed"
1537 if (ilen(tmpdir).gt.0)
1538 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1539 & //liczba(:ilen(liczba))//'.rst')
1540 inquire(file=rest2name,exist=file_exist)
1543 if(.not.restart1file) then
1544 if(me.eq.king.or..not.out1file)
1545 & write(iout,*) "Initial state will be read from file ",
1546 & rest2name(:ilen(rest2name))
1549 call rescale_weights(t_bath)
1551 if(me.eq.king.or..not.out1file)then
1552 if (restart1file) then
1553 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1556 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1559 write(iout,*) "Initial velocities randomly generated"
1565 c Generate initial velocities
1566 if(me.eq.king.or..not.out1file)
1567 & write(iout,*) "Initial velocities randomly generated"
1571 c rest2name = prefix(:ilen(prefix))//'.rst'
1572 if(me.eq.king.or..not.out1file)then
1573 write (iout,*) "Initial velocities"
1575 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1576 & (d_t(j,i+nres),j=1,3)
1578 c Zeroing the total angular momentum of the system
1579 write(iout,*) "Calling the zero-angular
1580 & momentum subroutine"
1583 c Getting the potential energy and forces and velocities and accelerations
1585 c write (iout,*) "velocity of the center of the mass:"
1586 c write (iout,*) (vcm(j),j=1,3)
1588 d_t(j,0)=d_t(j,0)-vcm(j)
1590 c Removing the velocity of the center of mass
1592 if(me.eq.king.or..not.out1file)then
1593 write (iout,*) "vcm right after adjustment:"
1594 write (iout,*) (vcm(j),j=1,3)
1598 if(iranconf.ne.0) then
1600 print *, 'Calling OVERLAP_SC'
1601 call overlap_sc(fail)
1605 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1606 print *,'SC_move',nft_sc,etot
1607 if(me.eq.king.or..not.out1file)
1608 & write(iout,*) 'SC_move',nft_sc,etot
1612 print *, 'Calling MINIM_DC'
1613 call minim_dc(etot,iretcode,nfun)
1615 call geom_to_var(nvar,varia)
1616 print *,'Calling MINIMIZE.'
1617 call minimize(etot,varia,iretcode,nfun)
1618 call var_to_geom(nvar,varia)
1620 if(me.eq.king.or..not.out1file)
1621 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1624 call chainbuild_cart
1629 kinetic_T=2.0d0/(dimen3*Rb)*EK
1630 if(me.eq.king.or..not.out1file)then
1640 call etotal(potEcomp)
1641 if (large) call enerprint(potEcomp)
1644 t_etotal=t_etotal+MPI_Wtime()-tt0
1646 t_etotal=t_etotal+tcpu()-tt0
1653 if (amax*d_time .gt. dvmax) then
1654 d_time=d_time*dvmax/amax
1655 if(me.eq.king.or..not.out1file) write (iout,*)
1656 & "Time step reduced to",d_time,
1657 & " because of too large initial acceleration."
1659 if(me.eq.king.or..not.out1file)then
1660 write(iout,*) "Potential energy and its components"
1661 call enerprint(potEcomp)
1662 c write(iout,*) (potEcomp(i),i=0,n_ene)
1664 potE=potEcomp(0)-potEcomp(20)
1667 if (ntwe.ne.0) call statout(itime)
1668 if(me.eq.king.or..not.out1file)
1669 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1670 & " Kinetic energy",EK," Potential energy",potE,
1671 & " Total energy",totE," Maximum acceleration ",
1674 write (iout,*) "Initial coordinates"
1676 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1677 & (c(j,i+nres),j=1,3)
1679 write (iout,*) "Initial dC"
1681 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1682 & (dc(j,i+nres),j=1,3)
1684 write (iout,*) "Initial velocities"
1685 write (iout,"(13x,' backbone ',23x,' side chain')")
1687 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1688 & (d_t(j,i+nres),j=1,3)
1690 write (iout,*) "Initial accelerations"
1692 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1693 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1694 & (d_a(j,i+nres),j=1,3)
1700 d_t_old(j,i)=d_t(j,i)
1701 d_a_old(j,i)=d_a(j,i)
1703 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1712 call etotal_short(energia_short)
1713 if (large) call enerprint(potEcomp)
1716 t_eshort=t_eshort+MPI_Wtime()-tt0
1718 t_eshort=t_eshort+tcpu()-tt0
1723 if(.not.out1file .and. large) then
1724 write (iout,*) "energia_long",energia_long(0),
1725 & " energia_short",energia_short(0),
1726 & " total",energia_long(0)+energia_short(0)
1727 write (iout,*) "Initial fast-force accelerations"
1729 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1730 & (d_a(j,i+nres),j=1,3)
1733 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1736 d_a_short(j,i)=d_a(j,i)
1745 call etotal_long(energia_long)
1746 if (large) call enerprint(potEcomp)
1749 t_elong=t_elong+MPI_Wtime()-tt0
1751 t_elong=t_elong+tcpu()-tt0
1756 if(.not.out1file .and. large) then
1757 write (iout,*) "energia_long",energia_long(0)
1758 write (iout,*) "Initial slow-force accelerations"
1760 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1761 & (d_a(j,i+nres),j=1,3)
1765 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1767 t_enegrad=t_enegrad+tcpu()-tt0
1772 c-----------------------------------------------------------
1773 subroutine random_vel
1774 implicit real*8 (a-h,o-z)
1775 include 'DIMENSIONS'
1776 include 'COMMON.CONTROL'
1777 include 'COMMON.VAR'
1780 include 'COMMON.LANGEVIN'
1782 include 'COMMON.LANGEVIN.lang0'
1784 include 'COMMON.CHAIN'
1785 include 'COMMON.DERIV'
1786 include 'COMMON.GEO'
1787 include 'COMMON.LOCAL'
1788 include 'COMMON.INTERACT'
1789 include 'COMMON.IOUNITS'
1790 include 'COMMON.NAMES'
1791 include 'COMMON.TIME1'
1792 double precision xv,sigv,lowb,highb
1793 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1794 c First generate velocities in the eigenspace of the G matrix
1795 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1802 sigv=dsqrt((Rb*t_bath)/geigen(i))
1805 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1806 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1807 c & " d_t_work_new",d_t_work_new(ii)
1816 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1819 c write (iout,*) "Ek from eigenvectors",Ek1
1821 c Transform velocities to UNRES coordinate space
1827 d_t_work(ind)=d_t_work(ind)
1828 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1830 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1834 c Transfer to the d_t vector
1836 d_t(j,0)=d_t_work(j)
1842 d_t(j,i)=d_t_work(ind)
1846 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1849 d_t(j,i+nres)=d_t_work(ind)
1854 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1855 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1860 c-----------------------------------------------------------
1861 subroutine sd_verlet_p_setup
1862 c Sets up the parameters of stochastic Verlet algorithm
1863 implicit real*8 (a-h,o-z)
1864 include 'DIMENSIONS'
1868 include 'COMMON.CONTROL'
1869 include 'COMMON.VAR'
1872 include 'COMMON.LANGEVIN'
1874 include 'COMMON.LANGEVIN.lang0'
1876 include 'COMMON.CHAIN'
1877 include 'COMMON.DERIV'
1878 include 'COMMON.GEO'
1879 include 'COMMON.LOCAL'
1880 include 'COMMON.INTERACT'
1881 include 'COMMON.IOUNITS'
1882 include 'COMMON.NAMES'
1883 include 'COMMON.TIME1'
1884 double precision emgdt(MAXRES6),
1885 & pterm,vterm,rho,rhoc,vsig,
1886 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1887 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1888 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1889 logical lprn /.false./
1890 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1891 double precision ktm
1898 c AL 8/17/04 Code adapted from tinker
1900 c Get the frictional and random terms for stochastic dynamics in the
1901 c eigenspace of mass-scaled UNRES friction matrix
1904 gdt = fricgam(i) * d_time
1906 c Stochastic dynamics reduces to simple MD for zero friction
1908 if (gdt .le. zero) then
1909 pfric_vec(i) = 1.0d0
1910 vfric_vec(i) = d_time
1911 afric_vec(i) = 0.5d0 * d_time * d_time
1912 prand_vec(i) = 0.0d0
1913 vrand_vec1(i) = 0.0d0
1914 vrand_vec2(i) = 0.0d0
1916 c Analytical expressions when friction coefficient is large
1919 if (gdt .ge. gdt_radius) then
1922 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1923 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1924 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1925 vterm = 1.0d0 - egdt**2
1926 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1928 c Use series expansions when friction coefficient is small
1939 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1940 & - gdt5/120.0d0 + gdt6/720.0d0
1941 & - gdt7/5040.0d0 + gdt8/40320.0d0
1942 & - gdt9/362880.0d0) / fricgam(i)**2
1943 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1944 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1945 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1946 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1947 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1948 & + 127.0d0*gdt9/90720.0d0
1949 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1950 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1951 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1952 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1953 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1954 & - 17.0d0*gdt2/1280.0d0
1955 & + 17.0d0*gdt3/6144.0d0
1956 & + 40967.0d0*gdt4/34406400.0d0
1957 & - 57203.0d0*gdt5/275251200.0d0
1958 & - 1429487.0d0*gdt6/13212057600.0d0)
1961 c Compute the scaling factors of random terms for the nonzero friction case
1963 ktm = 0.5d0*d_time/fricgam(i)
1964 psig = dsqrt(ktm*pterm) / fricgam(i)
1965 vsig = dsqrt(ktm*vterm)
1966 rhoc = dsqrt(1.0d0 - rho*rho)
1968 vrand_vec1(i) = vsig * rho
1969 vrand_vec2(i) = vsig * rhoc
1974 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1977 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1978 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1982 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1985 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1986 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1987 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1988 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1989 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1990 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1993 t_sdsetup=t_sdsetup+MPI_Wtime()
1995 t_sdsetup=t_sdsetup+tcpu()-tt0
1999 c-------------------------------------------------------------
2000 subroutine eigtransf1(n,ndim,ab,d,c)
2003 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2009 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2015 c-------------------------------------------------------------
2016 subroutine eigtransf(n,ndim,a,b,d,c)
2019 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2025 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2031 c-------------------------------------------------------------
2032 subroutine sd_verlet1
2033 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2034 implicit real*8 (a-h,o-z)
2035 include 'DIMENSIONS'
2036 include 'COMMON.CONTROL'
2037 include 'COMMON.VAR'
2040 include 'COMMON.LANGEVIN'
2042 include 'COMMON.LANGEVIN.lang0'
2044 include 'COMMON.CHAIN'
2045 include 'COMMON.DERIV'
2046 include 'COMMON.GEO'
2047 include 'COMMON.LOCAL'
2048 include 'COMMON.INTERACT'
2049 include 'COMMON.IOUNITS'
2050 include 'COMMON.NAMES'
2051 double precision stochforcvec(MAXRES6)
2052 common /stochcalc/ stochforcvec
2053 logical lprn /.false./
2055 c write (iout,*) "dc_old"
2057 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2058 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2061 dc_work(j)=dc_old(j,0)
2062 d_t_work(j)=d_t_old(j,0)
2063 d_a_work(j)=d_a_old(j,0)
2068 dc_work(ind+j)=dc_old(j,i)
2069 d_t_work(ind+j)=d_t_old(j,i)
2070 d_a_work(ind+j)=d_a_old(j,i)
2075 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2077 dc_work(ind+j)=dc_old(j,i+nres)
2078 d_t_work(ind+j)=d_t_old(j,i+nres)
2079 d_a_work(ind+j)=d_a_old(j,i+nres)
2087 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2091 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2092 & vfric_mat(i,j),afric_mat(i,j),
2093 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2101 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2102 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2103 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2104 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2106 d_t_work_new(i)=ddt1+0.5d0*ddt2
2107 d_t_work(i)=ddt1+ddt2
2112 d_t(j,0)=d_t_work(j)
2117 dc(j,i)=dc_work(ind+j)
2118 d_t(j,i)=d_t_work(ind+j)
2123 if (itype(i).ne.10) then
2126 dc(j,inres)=dc_work(ind+j)
2127 d_t(j,inres)=d_t_work(ind+j)
2134 c--------------------------------------------------------------------------
2135 subroutine sd_verlet2
2136 c Calculating the adjusted velocities for accelerations
2137 implicit real*8 (a-h,o-z)
2138 include 'DIMENSIONS'
2139 include 'COMMON.CONTROL'
2140 include 'COMMON.VAR'
2143 include 'COMMON.LANGEVIN'
2145 include 'COMMON.LANGEVIN.lang0'
2147 include 'COMMON.CHAIN'
2148 include 'COMMON.DERIV'
2149 include 'COMMON.GEO'
2150 include 'COMMON.LOCAL'
2151 include 'COMMON.INTERACT'
2152 include 'COMMON.IOUNITS'
2153 include 'COMMON.NAMES'
2154 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2155 common /stochcalc/ stochforcvec
2157 c Compute the stochastic forces which contribute to velocity change
2159 call stochastic_force(stochforcvecV)
2166 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2167 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2168 & vrand_mat2(i,j)*stochforcvecV(j)
2170 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2174 d_t(j,0)=d_t_work(j)
2179 d_t(j,i)=d_t_work(ind+j)
2184 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2187 d_t(j,inres)=d_t_work(ind+j)
2194 c-----------------------------------------------------------
2195 subroutine sd_verlet_ciccotti_setup
2196 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2198 implicit real*8 (a-h,o-z)
2199 include 'DIMENSIONS'
2203 include 'COMMON.CONTROL'
2204 include 'COMMON.VAR'
2207 include 'COMMON.LANGEVIN'
2209 include 'COMMON.LANGEVIN.lang0'
2211 include 'COMMON.CHAIN'
2212 include 'COMMON.DERIV'
2213 include 'COMMON.GEO'
2214 include 'COMMON.LOCAL'
2215 include 'COMMON.INTERACT'
2216 include 'COMMON.IOUNITS'
2217 include 'COMMON.NAMES'
2218 include 'COMMON.TIME1'
2219 double precision emgdt(MAXRES6),
2220 & pterm,vterm,rho,rhoc,vsig,
2221 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2222 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2223 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2224 logical lprn /.false./
2225 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2226 double precision ktm
2233 c AL 8/17/04 Code adapted from tinker
2235 c Get the frictional and random terms for stochastic dynamics in the
2236 c eigenspace of mass-scaled UNRES friction matrix
2239 write (iout,*) "i",i," fricgam",fricgam(i)
2240 gdt = fricgam(i) * d_time
2242 c Stochastic dynamics reduces to simple MD for zero friction
2244 if (gdt .le. zero) then
2245 pfric_vec(i) = 1.0d0
2246 vfric_vec(i) = d_time
2247 afric_vec(i) = 0.5d0*d_time*d_time
2248 prand_vec(i) = afric_vec(i)
2249 vrand_vec2(i) = vfric_vec(i)
2251 c Analytical expressions when friction coefficient is large
2256 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2257 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2258 prand_vec(i) = afric_vec(i)
2259 vrand_vec2(i) = vfric_vec(i)
2261 c Compute the scaling factors of random terms for the nonzero friction case
2263 c ktm = 0.5d0*d_time/fricgam(i)
2264 c psig = dsqrt(ktm*pterm) / fricgam(i)
2265 c vsig = dsqrt(ktm*vterm)
2266 c prand_vec(i) = psig*afric_vec(i)
2267 c vrand_vec2(i) = vsig*vfric_vec(i)
2272 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2275 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2276 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2280 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2282 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2283 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2284 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2285 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2286 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2288 t_sdsetup=t_sdsetup+MPI_Wtime()
2290 t_sdsetup=t_sdsetup+tcpu()-tt0
2294 c-------------------------------------------------------------
2295 subroutine sd_verlet1_ciccotti
2296 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2297 implicit real*8 (a-h,o-z)
2298 include 'DIMENSIONS'
2302 include 'COMMON.CONTROL'
2303 include 'COMMON.VAR'
2306 include 'COMMON.LANGEVIN'
2308 include 'COMMON.LANGEVIN.lang0'
2310 include 'COMMON.CHAIN'
2311 include 'COMMON.DERIV'
2312 include 'COMMON.GEO'
2313 include 'COMMON.LOCAL'
2314 include 'COMMON.INTERACT'
2315 include 'COMMON.IOUNITS'
2316 include 'COMMON.NAMES'
2317 double precision stochforcvec(MAXRES6)
2318 common /stochcalc/ stochforcvec
2319 logical lprn /.false./
2321 c write (iout,*) "dc_old"
2323 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2324 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2327 dc_work(j)=dc_old(j,0)
2328 d_t_work(j)=d_t_old(j,0)
2329 d_a_work(j)=d_a_old(j,0)
2334 dc_work(ind+j)=dc_old(j,i)
2335 d_t_work(ind+j)=d_t_old(j,i)
2336 d_a_work(ind+j)=d_a_old(j,i)
2341 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2343 dc_work(ind+j)=dc_old(j,i+nres)
2344 d_t_work(ind+j)=d_t_old(j,i+nres)
2345 d_a_work(ind+j)=d_a_old(j,i+nres)
2354 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2358 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2359 & vfric_mat(i,j),afric_mat(i,j),
2360 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2368 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2369 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2370 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2371 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2373 d_t_work_new(i)=ddt1+0.5d0*ddt2
2374 d_t_work(i)=ddt1+ddt2
2379 d_t(j,0)=d_t_work(j)
2384 dc(j,i)=dc_work(ind+j)
2385 d_t(j,i)=d_t_work(ind+j)
2390 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2393 dc(j,inres)=dc_work(ind+j)
2394 d_t(j,inres)=d_t_work(ind+j)
2401 c--------------------------------------------------------------------------
2402 subroutine sd_verlet2_ciccotti
2403 c Calculating the adjusted velocities for accelerations
2404 implicit real*8 (a-h,o-z)
2405 include 'DIMENSIONS'
2406 include 'COMMON.CONTROL'
2407 include 'COMMON.VAR'
2410 include 'COMMON.LANGEVIN'
2412 include 'COMMON.LANGEVIN.lang0'
2414 include 'COMMON.CHAIN'
2415 include 'COMMON.DERIV'
2416 include 'COMMON.GEO'
2417 include 'COMMON.LOCAL'
2418 include 'COMMON.INTERACT'
2419 include 'COMMON.IOUNITS'
2420 include 'COMMON.NAMES'
2421 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2422 common /stochcalc/ stochforcvec
2424 c Compute the stochastic forces which contribute to velocity change
2426 call stochastic_force(stochforcvecV)
2433 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2434 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2435 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2437 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2441 d_t(j,0)=d_t_work(j)
2446 d_t(j,i)=d_t_work(ind+j)
2451 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2454 d_t(j,inres)=d_t_work(ind+j)