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 C call check_ecartint
235 write (tytul,'("time",f8.2)') totT
237 call hairpin(.true.,nharp,iharp)
238 call secondary2(.true.)
239 call pdbout(potE,tytul,ipdb)
244 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
245 open(irest2,file=rest2name,status='unknown')
246 write(irest2,*) totT,EK,potE,totE,t_bath
248 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
251 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
262 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
264 & 'MD calculations setup:',t_MDsetup,
265 & 'Energy & gradient evaluation:',t_enegrad,
266 & 'Stochastic MD setup:',t_langsetup,
267 & 'Stochastic MD step setup:',t_sdsetup,
269 write (iout,'(/28(1h=),a25,27(1h=))')
270 & ' End of MD calculation '
272 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
274 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
275 & " time_fricmatmult",time_fricmatmult," time_fsample ",
280 c-------------------------------------------------------------------------------
281 subroutine velverlet_step(itime)
282 c-------------------------------------------------------------------------------
283 c Perform a single velocity Verlet step; the time step can be rescaled if
284 c increments in accelerations exceed the threshold
285 c-------------------------------------------------------------------------------
286 implicit real*8 (a-h,o-z)
290 integer ierror,ierrcode
292 include 'COMMON.SETUP'
293 include 'COMMON.CONTROL'
297 include 'COMMON.LANGEVIN'
299 include 'COMMON.LANGEVIN.lang0'
301 include 'COMMON.CHAIN'
302 include 'COMMON.DERIV'
304 include 'COMMON.LOCAL'
305 include 'COMMON.INTERACT'
306 include 'COMMON.IOUNITS'
307 include 'COMMON.NAMES'
308 include 'COMMON.TIME1'
309 include 'COMMON.MUCA'
310 double precision vcm(3),incr(3)
311 double precision cm(3),L(3)
312 integer ilen,count,rstcount
315 integer maxcount_scale /20/
317 double precision stochforcvec(MAXRES6)
318 common /stochcalc/ stochforcvec
326 else if (lang.eq.2 .or. lang.eq.3) then
328 call stochastic_force(stochforcvec)
331 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
333 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
340 icount_scale=icount_scale+1
341 if (icount_scale.gt.maxcount_scale) then
343 & "ERROR: too many attempts at scaling down the time step. ",
344 & "amax=",amax,"epdrift=",epdrift,
345 & "damax=",damax,"edriftmax=",edriftmax,
349 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
353 c First step of the velocity Verlet algorithm
358 else if (lang.eq.3) then
360 call sd_verlet1_ciccotti
362 else if (lang.eq.1) then
367 c Build the chain from the newly calculated coordinates
369 if (rattle) call rattle1
371 if (large.and. mod(itime,ntwe).eq.0) then
372 write (iout,*) "Cartesian and internal coordinates: step 1"
377 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
378 & (dc(j,i+nres),j=1,3)
380 write (iout,*) "Accelerations"
382 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
383 & (d_a(j,i+nres),j=1,3)
385 write (iout,*) "Velocities, step 1"
387 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
388 & (d_t(j,i+nres),j=1,3)
397 c Calculate energy and forces
399 call etotal(potEcomp)
400 if (large.and. mod(itime,ntwe).eq.0)
401 & call enerprint(potEcomp)
404 t_etotal=t_etotal+MPI_Wtime()-tt0
406 t_etotal=t_etotal+tcpu()-tt0
409 potE=potEcomp(0)-potEcomp(20)
411 c Get the new accelerations
414 t_enegrad=t_enegrad+MPI_Wtime()-tt0
416 t_enegrad=t_enegrad+tcpu()-tt0
418 c Determine maximum acceleration and scale down the timestep if needed
420 amax=amax/(itime_scal+1)**2
421 call predict_edrift(epdrift)
422 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
423 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
425 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
427 itime_scal=itime_scal+ifac_time
428 c fac_time=dmin1(damax/amax,0.5d0)
429 fac_time=0.5d0**ifac_time
430 d_time=d_time*fac_time
431 if (lang.eq.2 .or. lang.eq.3) then
433 c write (iout,*) "Calling sd_verlet_setup: 1"
434 c Rescale the stochastic forces and recalculate or restore
435 c the matrices of tinker integrator
436 if (itime_scal.gt.maxflag_stoch) then
437 if (large) write (iout,'(a,i5,a)')
438 & "Calculate matrices for stochastic step;",
439 & " itime_scal ",itime_scal
441 call sd_verlet_p_setup
443 call sd_verlet_ciccotti_setup
445 write (iout,'(2a,i3,a,i3,1h.)')
446 & "Warning: cannot store matrices for stochastic",
447 & " integration because the index",itime_scal,
448 & " is greater than",maxflag_stoch
449 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
450 & " integration Langevin algorithm for better efficiency."
451 else if (flag_stoch(itime_scal)) then
452 if (large) write (iout,'(a,i5,a,l1)')
453 & "Restore matrices for stochastic step; itime_scal ",
454 & itime_scal," flag ",flag_stoch(itime_scal)
457 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
458 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
459 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
460 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
461 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
462 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
466 if (large) write (iout,'(2a,i5,a,l1)')
467 & "Calculate & store matrices for stochastic step;",
468 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
470 call sd_verlet_p_setup
472 call sd_verlet_ciccotti_setup
474 flag_stoch(ifac_time)=.true.
477 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
478 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
479 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
480 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
481 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
482 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
486 fac_time=1.0d0/dsqrt(fac_time)
488 stochforcvec(i)=fac_time*stochforcvec(i)
491 else if (lang.eq.1) then
492 c Rescale the accelerations due to stochastic forces
493 fac_time=1.0d0/dsqrt(fac_time)
495 d_as_work(i)=d_as_work(i)*fac_time
498 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
499 & "itime",itime," Timestep scaled down to ",
500 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
502 c Second step of the velocity Verlet algorithm
507 else if (lang.eq.3) then
509 call sd_verlet2_ciccotti
511 else if (lang.eq.1) then
516 if (rattle) call rattle2
518 if (d_time.ne.d_time0) then
521 if (lang.eq.2 .or. lang.eq.3) then
522 if (large) write (iout,'(a)')
523 & "Restore original matrices for stochastic step"
524 c write (iout,*) "Calling sd_verlet_setup: 2"
525 c Restore the matrices of tinker integrator if the time step has been restored
528 pfric_mat(i,j)=pfric0_mat(i,j,0)
529 afric_mat(i,j)=afric0_mat(i,j,0)
530 vfric_mat(i,j)=vfric0_mat(i,j,0)
531 prand_mat(i,j)=prand0_mat(i,j,0)
532 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
533 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
542 c Calculate the kinetic and the total energy and the kinetic temperature
547 c write (iout,*) "step",itime," EK",EK," EK1",EK1
549 c Couple the system to Berendsen bath if needed
550 if (tbf .and. lang.eq.0) then
553 kinetic_T=2.0d0/(dimen3*Rb)*EK
554 c Backup the coordinates, velocities, and accelerations
558 d_t_old(j,i)=d_t(j,i)
559 d_a_old(j,i)=d_a(j,i)
563 if (mod(itime,ntwe).eq.0 .and. large) then
564 write (iout,*) "Velocities, step 2"
566 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
567 & (d_t(j,i+nres),j=1,3)
573 c-------------------------------------------------------------------------------
574 subroutine RESPA_step(itime)
575 c-------------------------------------------------------------------------------
576 c Perform a single RESPA step.
577 c-------------------------------------------------------------------------------
578 implicit real*8 (a-h,o-z)
582 integer IERROR,ERRCODE
584 include 'COMMON.SETUP'
585 include 'COMMON.CONTROL'
589 include 'COMMON.LANGEVIN'
591 include 'COMMON.LANGEVIN.lang0'
593 include 'COMMON.CHAIN'
594 include 'COMMON.DERIV'
596 include 'COMMON.LOCAL'
597 include 'COMMON.INTERACT'
598 include 'COMMON.IOUNITS'
599 include 'COMMON.NAMES'
600 include 'COMMON.TIME1'
601 double precision energia_short(0:n_ene),
602 & energia_long(0:n_ene)
603 double precision cm(3),L(3),vcm(3),incr(3)
604 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
605 & d_a_old0(3,0:maxres2)
606 logical PRINT_AMTS_MSG /.false./
607 integer ilen,count,rstcount
610 integer maxcount_scale /10/
612 double precision stochforcvec(MAXRES6)
613 common /stochcalc/ stochforcvec
616 common /cipiszcze/ itt
619 if (large.and. mod(itime,ntwe).eq.0) then
620 write (iout,*) "***************** RESPA itime",itime
621 write (iout,*) "Cartesian and internal coordinates: step 0"
623 call pdbout(0.0d0,"cipiszcze",iout)
625 write (iout,*) "Accelerations from long-range forces"
627 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
628 & (d_a(j,i+nres),j=1,3)
630 write (iout,*) "Velocities, step 0"
632 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
633 & (d_t(j,i+nres),j=1,3)
638 c Perform the initial RESPA step (increment velocities)
639 c write (iout,*) "*********************** RESPA ini"
642 if (mod(itime,ntwe).eq.0 .and. large) then
643 write (iout,*) "Velocities, end"
645 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
646 & (d_t(j,i+nres),j=1,3)
650 c Compute the short-range forces
656 C 7/2/2009 commented out
658 c call etotal_short(energia_short)
661 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
664 d_a(j,i)=d_a_short(j,i)
668 if (large.and. mod(itime,ntwe).eq.0) then
669 write (iout,*) "energia_short",energia_short(0)
670 write (iout,*) "Accelerations from short-range forces"
672 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
673 & (d_a(j,i+nres),j=1,3)
678 t_enegrad=t_enegrad+MPI_Wtime()-tt0
680 t_enegrad=t_enegrad+tcpu()-tt0
685 d_t_old(j,i)=d_t(j,i)
686 d_a_old(j,i)=d_a(j,i)
689 c 6/30/08 A-MTS: attempt at increasing the split number
692 dc_old0(j,i)=dc_old(j,i)
693 d_t_old0(j,i)=d_t_old(j,i)
694 d_a_old0(j,i)=d_a_old(j,i)
697 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
698 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
705 c write (iout,*) "itime",itime," ntime_split",ntime_split
706 c Split the time step
707 d_time=d_time0/ntime_split
708 c Perform the short-range RESPA steps (velocity Verlet increments of
709 c positions and velocities using short-range forces)
710 c write (iout,*) "*********************** RESPA split"
711 do itsplit=1,ntime_split
714 else if (lang.eq.2 .or. lang.eq.3) then
716 call stochastic_force(stochforcvec)
719 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
721 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
726 c First step of the velocity Verlet algorithm
731 else if (lang.eq.3) then
733 call sd_verlet1_ciccotti
735 else if (lang.eq.1) then
740 c Build the chain from the newly calculated coordinates
742 if (rattle) call rattle1
744 if (large.and. mod(itime,ntwe).eq.0) then
745 write (iout,*) "***** ITSPLIT",itsplit
746 write (iout,*) "Cartesian and internal coordinates: step 1"
747 call pdbout(0.0d0,"cipiszcze",iout)
750 write (iout,*) "Velocities, step 1"
752 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
753 & (d_t(j,i+nres),j=1,3)
762 c Calculate energy and forces
764 call etotal_short(energia_short)
765 if (large.and. mod(itime,ntwe).eq.0)
766 & call enerprint(energia_short)
769 t_eshort=t_eshort+MPI_Wtime()-tt0
771 t_eshort=t_eshort+tcpu()-tt0
775 c Get the new accelerations
777 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
780 d_a_short(j,i)=d_a(j,i)
784 if (large.and. mod(itime,ntwe).eq.0) then
785 write (iout,*)"energia_short",energia_short(0)
786 write (iout,*) "Accelerations from short-range forces"
788 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
789 & (d_a(j,i+nres),j=1,3)
794 c Determine maximum acceleration and scale down the timestep if needed
796 amax=amax/ntime_split**2
797 call predict_edrift(epdrift)
798 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
799 & write (iout,*) "amax",amax," damax",damax,
800 & " epdrift",epdrift," epdriftmax",epdriftmax
801 c Exit loop and try with increased split number if the change of
802 c acceleration is too big
803 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
804 if (ntime_split.lt.maxtime_split) then
806 ntime_split=ntime_split*2
809 dc_old(j,i)=dc_old0(j,i)
810 d_t_old(j,i)=d_t_old0(j,i)
811 d_a_old(j,i)=d_a_old0(j,i)
814 if (PRINT_AMTS_MSG) then
815 write (iout,*) "acceleration/energy drift too large",amax,
816 & epdrift," split increased to ",ntime_split," itime",itime,
822 & "Uh-hu. Bumpy landscape. Maximum splitting number",
824 & " already reached!!! Trying to carry on!"
828 t_enegrad=t_enegrad+MPI_Wtime()-tt0
830 t_enegrad=t_enegrad+tcpu()-tt0
832 c Second step of the velocity Verlet algorithm
837 else if (lang.eq.3) then
839 call sd_verlet2_ciccotti
841 else if (lang.eq.1) then
846 if (rattle) call rattle2
847 c Backup the coordinates, velocities, and accelerations
851 d_t_old(j,i)=d_t(j,i)
852 d_a_old(j,i)=d_a(j,i)
859 c Restore the time step
861 c Compute long-range forces
868 call etotal_long(energia_long)
869 if (large.and. mod(itime,ntwe).eq.0)
870 & call enerprint(energia_long)
873 t_elong=t_elong+MPI_Wtime()-tt0
875 t_elong=t_elong+tcpu()-tt0
881 t_enegrad=t_enegrad+MPI_Wtime()-tt0
883 t_enegrad=t_enegrad+tcpu()-tt0
885 c Compute accelerations from long-range forces
887 if (large.and. mod(itime,ntwe).eq.0) then
888 write (iout,*) "energia_long",energia_long(0)
889 write (iout,*) "Cartesian and internal coordinates: step 2"
891 call pdbout(0.0d0,"cipiszcze",iout)
893 write (iout,*) "Accelerations from long-range forces"
895 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
896 & (d_a(j,i+nres),j=1,3)
898 write (iout,*) "Velocities, step 2"
900 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
901 & (d_t(j,i+nres),j=1,3)
905 c Compute the final RESPA step (increment velocities)
906 c write (iout,*) "*********************** RESPA fin"
908 c Compute the complete potential energy
910 potEcomp(i)=energia_short(i)+energia_long(i)
912 potE=potEcomp(0)-potEcomp(20)
913 c potE=energia_short(0)+energia_long(0)
915 c Calculate the kinetic and the total energy and the kinetic temperature
918 c Couple the system to Berendsen bath if needed
919 if (tbf .and. lang.eq.0) then
922 kinetic_T=2.0d0/(dimen3*Rb)*EK
923 c Backup the coordinates, velocities, and accelerations
925 if (mod(itime,ntwe).eq.0 .and. large) then
926 write (iout,*) "Velocities, end"
928 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
929 & (d_t(j,i+nres),j=1,3)
935 c---------------------------------------------------------------------
937 c First and last RESPA step (incrementing velocities using long-range
939 implicit real*8 (a-h,o-z)
941 include 'COMMON.CONTROL'
944 include 'COMMON.CHAIN'
945 include 'COMMON.DERIV'
947 include 'COMMON.LOCAL'
948 include 'COMMON.INTERACT'
949 include 'COMMON.IOUNITS'
950 include 'COMMON.NAMES'
952 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
956 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
960 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
963 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
969 c-----------------------------------------------------------------
971 c Applying velocity Verlet algorithm - step 1 to coordinates
972 implicit real*8 (a-h,o-z)
974 include 'COMMON.CONTROL'
977 include 'COMMON.CHAIN'
978 include 'COMMON.DERIV'
980 include 'COMMON.LOCAL'
981 include 'COMMON.INTERACT'
982 include 'COMMON.IOUNITS'
983 include 'COMMON.NAMES'
984 double precision adt,adt2
987 write (iout,*) "VELVERLET1 START: DC"
989 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
990 & (dc(j,i+nres),j=1,3)
994 adt=d_a_old(j,0)*d_time
996 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
997 d_t_new(j,0)=d_t_old(j,0)+adt2
998 d_t(j,0)=d_t_old(j,0)+adt
1004 adt=d_a_old(j,i)*d_time
1006 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1007 d_t_new(j,i)=d_t_old(j,i)+adt2
1008 d_t(j,i)=d_t_old(j,i)+adt
1013 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1016 adt=d_a_old(j,inres)*d_time
1018 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1019 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1020 d_t(j,inres)=d_t_old(j,inres)+adt
1025 write (iout,*) "VELVERLET1 END: DC"
1027 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1028 & (dc(j,i+nres),j=1,3)
1033 c---------------------------------------------------------------------
1035 c Step 2 of the velocity Verlet algorithm: update velocities
1036 implicit real*8 (a-h,o-z)
1037 include 'DIMENSIONS'
1038 include 'COMMON.CONTROL'
1039 include 'COMMON.VAR'
1041 include 'COMMON.CHAIN'
1042 include 'COMMON.DERIV'
1043 include 'COMMON.GEO'
1044 include 'COMMON.LOCAL'
1045 include 'COMMON.INTERACT'
1046 include 'COMMON.IOUNITS'
1047 include 'COMMON.NAMES'
1049 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1053 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1057 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1060 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1066 c-----------------------------------------------------------------
1067 subroutine sddir_precalc
1068 c Applying velocity Verlet algorithm - step 1 to coordinates
1069 implicit real*8 (a-h,o-z)
1070 include 'DIMENSIONS'
1074 include 'COMMON.CONTROL'
1075 include 'COMMON.VAR'
1078 include 'COMMON.LANGEVIN'
1080 include 'COMMON.LANGEVIN.lang0'
1082 include 'COMMON.CHAIN'
1083 include 'COMMON.DERIV'
1084 include 'COMMON.GEO'
1085 include 'COMMON.LOCAL'
1086 include 'COMMON.INTERACT'
1087 include 'COMMON.IOUNITS'
1088 include 'COMMON.NAMES'
1089 include 'COMMON.TIME1'
1090 double precision stochforcvec(MAXRES6)
1091 common /stochcalc/ stochforcvec
1093 c Compute friction and stochastic forces
1102 time_fric=time_fric+MPI_Wtime()-time00
1105 time_fric=time_fric+tcpu()-time00
1108 call stochastic_force(stochforcvec)
1110 time_stoch=time_stoch+MPI_Wtime()-time00
1112 time_stoch=time_stoch+tcpu()-time00
1115 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1116 c forces (d_as_work)
1118 call ginv_mult(fric_work, d_af_work)
1119 call ginv_mult(stochforcvec, d_as_work)
1122 c---------------------------------------------------------------------
1123 subroutine sddir_verlet1
1124 c Applying velocity Verlet algorithm - step 1 to velocities
1125 implicit real*8 (a-h,o-z)
1126 include 'DIMENSIONS'
1127 include 'COMMON.CONTROL'
1128 include 'COMMON.VAR'
1131 include 'COMMON.LANGEVIN'
1133 include 'COMMON.LANGEVIN.lang0'
1135 include 'COMMON.CHAIN'
1136 include 'COMMON.DERIV'
1137 include 'COMMON.GEO'
1138 include 'COMMON.LOCAL'
1139 include 'COMMON.INTERACT'
1140 include 'COMMON.IOUNITS'
1141 include 'COMMON.NAMES'
1142 c Revised 3/31/05 AL: correlation between random contributions to
1143 c position and velocity increments included.
1144 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1145 double precision adt,adt2
1147 c Add the contribution from BOTH friction and stochastic force to the
1148 c coordinates, but ONLY the contribution from the friction forces to velocities
1151 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1152 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1153 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1154 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1155 d_t(j,0)=d_t_old(j,0)+adt
1160 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1161 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1162 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1163 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1164 d_t(j,i)=d_t_old(j,i)+adt
1169 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1172 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1173 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1174 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1175 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1176 d_t(j,inres)=d_t_old(j,inres)+adt
1183 c---------------------------------------------------------------------
1184 subroutine sddir_verlet2
1185 c Calculating the adjusted velocities for accelerations
1186 implicit real*8 (a-h,o-z)
1187 include 'DIMENSIONS'
1188 include 'COMMON.CONTROL'
1189 include 'COMMON.VAR'
1192 include 'COMMON.LANGEVIN'
1194 include 'COMMON.LANGEVIN.lang0'
1196 include 'COMMON.CHAIN'
1197 include 'COMMON.DERIV'
1198 include 'COMMON.GEO'
1199 include 'COMMON.LOCAL'
1200 include 'COMMON.INTERACT'
1201 include 'COMMON.IOUNITS'
1202 include 'COMMON.NAMES'
1203 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1204 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1205 c Revised 3/31/05 AL: correlation between random contributions to
1206 c position and velocity increments included.
1207 c The correlation coefficients are calculated at low-friction limit.
1208 c Also, friction forces are now not calculated with new velocities.
1210 c call friction_force
1211 call stochastic_force(stochforcvec)
1213 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1214 c forces (d_as_work)
1216 call ginv_mult(stochforcvec, d_as_work1)
1222 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1223 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1228 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1229 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1234 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1237 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1238 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1239 & +cos60*d_as_work1(ind+j))*d_time
1246 c---------------------------------------------------------------------
1247 subroutine max_accel
1249 c Find the maximum difference in the accelerations of the the sites
1250 c at the beginning and the end of the time step.
1252 implicit real*8 (a-h,o-z)
1253 include 'DIMENSIONS'
1254 include 'COMMON.CONTROL'
1255 include 'COMMON.VAR'
1257 include 'COMMON.CHAIN'
1258 include 'COMMON.DERIV'
1259 include 'COMMON.GEO'
1260 include 'COMMON.LOCAL'
1261 include 'COMMON.INTERACT'
1262 include 'COMMON.IOUNITS'
1263 double precision aux(3),accel(3),accel_old(3),dacc
1265 c aux(j)=d_a(j,0)-d_a_old(j,0)
1266 accel_old(j)=d_a_old(j,0)
1273 c 7/3/08 changed to asymmetric difference
1275 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1276 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1277 accel(j)=accel(j)+0.5d0*d_a(j,i)
1278 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1279 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1280 dacc=dabs(accel(j)-accel_old(j))
1281 c write (iout,*) i,dacc
1282 if (dacc.gt.amax) amax=dacc
1290 accel_old(j)=d_a_old(j,0)
1295 accel_old(j)=accel_old(j)+d_a_old(j,1)
1296 accel(j)=accel(j)+d_a(j,1)
1300 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1302 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1303 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1304 accel(j)=accel(j)+d_a(j,i+nres)
1308 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1309 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1310 dacc=dabs(accel(j)-accel_old(j))
1311 c write (iout,*) "side-chain",i,dacc
1312 if (dacc.gt.amax) amax=dacc
1316 accel_old(j)=accel_old(j)+d_a_old(j,i)
1317 accel(j)=accel(j)+d_a(j,i)
1318 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1323 c---------------------------------------------------------------------
1324 subroutine predict_edrift(epdrift)
1326 c Predict the drift of the potential energy
1328 implicit real*8 (a-h,o-z)
1329 include 'DIMENSIONS'
1330 include 'COMMON.CONTROL'
1331 include 'COMMON.VAR'
1333 include 'COMMON.CHAIN'
1334 include 'COMMON.DERIV'
1335 include 'COMMON.GEO'
1336 include 'COMMON.LOCAL'
1337 include 'COMMON.INTERACT'
1338 include 'COMMON.IOUNITS'
1339 include 'COMMON.MUCA'
1340 double precision epdrift,epdriftij
1341 c Drift of the potential energy
1347 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1348 if (lmuca) epdriftij=epdriftij*factor
1349 c write (iout,*) "back",i,j,epdriftij
1350 if (epdriftij.gt.epdrift) epdrift=epdriftij
1354 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1357 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1358 if (lmuca) epdriftij=epdriftij*factor
1359 c write (iout,*) "side",i,j,epdriftij
1360 if (epdriftij.gt.epdrift) epdrift=epdriftij
1364 epdrift=0.5d0*epdrift*d_time*d_time
1365 c write (iout,*) "epdrift",epdrift
1368 c-----------------------------------------------------------------------
1369 subroutine verlet_bath
1371 c Coupling to the thermostat by using the Berendsen algorithm
1373 implicit real*8 (a-h,o-z)
1374 include 'DIMENSIONS'
1375 include 'COMMON.CONTROL'
1376 include 'COMMON.VAR'
1378 include 'COMMON.CHAIN'
1379 include 'COMMON.DERIV'
1380 include 'COMMON.GEO'
1381 include 'COMMON.LOCAL'
1382 include 'COMMON.INTERACT'
1383 include 'COMMON.IOUNITS'
1384 include 'COMMON.NAMES'
1385 double precision T_half,fact
1387 T_half=2.0d0/(dimen3*Rb)*EK
1388 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1389 c write(iout,*) "T_half", T_half
1390 c write(iout,*) "EK", EK
1391 c write(iout,*) "fact", fact
1393 d_t(j,0)=fact*d_t(j,0)
1397 d_t(j,i)=fact*d_t(j,i)
1401 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1404 d_t(j,inres)=fact*d_t(j,inres)
1410 c---------------------------------------------------------
1412 c Set up the initial conditions of a MD simulation
1413 implicit real*8 (a-h,o-z)
1414 include 'DIMENSIONS'
1418 integer IERROR,ERRCODE
1420 include 'COMMON.SETUP'
1421 include 'COMMON.CONTROL'
1422 include 'COMMON.VAR'
1425 include 'COMMON.LANGEVIN'
1427 include 'COMMON.LANGEVIN.lang0'
1429 include 'COMMON.CHAIN'
1430 include 'COMMON.DERIV'
1431 include 'COMMON.GEO'
1432 include 'COMMON.LOCAL'
1433 include 'COMMON.INTERACT'
1434 include 'COMMON.IOUNITS'
1435 include 'COMMON.NAMES'
1436 include 'COMMON.REMD'
1437 real*8 energia_long(0:n_ene),
1438 & energia_short(0:n_ene),vcm(3),incr(3)
1439 double precision cm(3),L(3),xv,sigv,lowb,highb
1440 double precision varia(maxvar)
1448 c write(iout,*) "d_time", d_time
1449 c Compute the standard deviations of stochastic forces for Langevin dynamics
1450 c if the friction coefficients do not depend on surface area
1451 if (lang.gt.0 .and. .not.surfarea) then
1453 stdforcp(i)=stdfp*dsqrt(gamp)
1456 stdforcsc(i)=stdfsc(iabs(itype(i)))
1457 & *dsqrt(gamsc(iabs(itype(i))))
1460 c Open the pdb file for snapshotshots
1463 if (ilen(tmpdir).gt.0)
1464 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1465 & liczba(:ilen(liczba))//".pdb")
1467 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1471 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1472 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1473 & liczba(:ilen(liczba))//".x")
1474 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1477 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1478 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1479 & liczba(:ilen(liczba))//".cx")
1480 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1486 if (ilen(tmpdir).gt.0)
1487 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1488 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1490 if (ilen(tmpdir).gt.0)
1491 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1492 cartname=prefix(:ilen(prefix))//"_MD.cx"
1496 write (qstr,'(256(1h ))')
1499 iq = qinfrag(i,iset)*10
1500 iw = wfrag(i,iset)/100
1502 if(me.eq.king.or..not.out1file)
1503 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1504 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1509 iq = qinpair(i,iset)*10
1510 iw = wpair(i,iset)/100
1512 if(me.eq.king.or..not.out1file)
1513 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1514 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1518 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1520 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1522 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1524 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1528 if (restart1file) then
1530 & inquire(file=mremd_rst_name,exist=file_exist)
1532 write (*,*) me," Before broadcast: file_exist",file_exist
1533 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1535 write (*,*) me," After broadcast: file_exist",file_exist
1536 c inquire(file=mremd_rst_name,exist=file_exist)
1538 if(me.eq.king.or..not.out1file)
1539 & write(iout,*) "Initial state read by master and distributed"
1541 if (ilen(tmpdir).gt.0)
1542 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1543 & //liczba(:ilen(liczba))//'.rst')
1544 inquire(file=rest2name,exist=file_exist)
1547 if(.not.restart1file) then
1548 if(me.eq.king.or..not.out1file)
1549 & write(iout,*) "Initial state will be read from file ",
1550 & rest2name(:ilen(rest2name))
1553 call rescale_weights(t_bath)
1555 if(me.eq.king.or..not.out1file)then
1556 if (restart1file) then
1557 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1560 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1563 write(iout,*) "Initial velocities randomly generated"
1569 c Generate initial velocities
1570 if(me.eq.king.or..not.out1file)
1571 & write(iout,*) "Initial velocities randomly generated"
1575 c rest2name = prefix(:ilen(prefix))//'.rst'
1576 if(me.eq.king.or..not.out1file)then
1577 write (iout,*) "Initial velocities"
1579 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1580 & (d_t(j,i+nres),j=1,3)
1582 c Zeroing the total angular momentum of the system
1583 write(iout,*) "Calling the zero-angular
1584 & momentum subroutine"
1587 c Getting the potential energy and forces and velocities and accelerations
1589 c write (iout,*) "velocity of the center of the mass:"
1590 c write (iout,*) (vcm(j),j=1,3)
1592 d_t(j,0)=d_t(j,0)-vcm(j)
1594 c Removing the velocity of the center of mass
1596 if(me.eq.king.or..not.out1file)then
1597 write (iout,*) "vcm right after adjustment:"
1598 write (iout,*) (vcm(j),j=1,3)
1602 if(iranconf.ne.0) then
1604 print *, 'Calling OVERLAP_SC'
1605 call overlap_sc(fail)
1609 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1610 print *,'SC_move',nft_sc,etot
1611 if(me.eq.king.or..not.out1file)
1612 & write(iout,*) 'SC_move',nft_sc,etot
1616 print *, 'Calling MINIM_DC'
1617 call minim_dc(etot,iretcode,nfun)
1619 call geom_to_var(nvar,varia)
1620 print *,'Calling MINIMIZE.'
1621 call minimize(etot,varia,iretcode,nfun)
1622 call var_to_geom(nvar,varia)
1624 if(me.eq.king.or..not.out1file)
1625 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1628 call chainbuild_cart
1633 kinetic_T=2.0d0/(dimen3*Rb)*EK
1634 if(me.eq.king.or..not.out1file)then
1644 call etotal(potEcomp)
1645 if (large) call enerprint(potEcomp)
1648 t_etotal=t_etotal+MPI_Wtime()-tt0
1650 t_etotal=t_etotal+tcpu()-tt0
1657 if (amax*d_time .gt. dvmax) then
1658 d_time=d_time*dvmax/amax
1659 if(me.eq.king.or..not.out1file) write (iout,*)
1660 & "Time step reduced to",d_time,
1661 & " because of too large initial acceleration."
1663 if(me.eq.king.or..not.out1file)then
1664 write(iout,*) "Potential energy and its components"
1665 call enerprint(potEcomp)
1666 c write(iout,*) (potEcomp(i),i=0,n_ene)
1668 potE=potEcomp(0)-potEcomp(20)
1671 if (ntwe.ne.0) call statout(itime)
1672 if(me.eq.king.or..not.out1file)
1673 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1674 & " Kinetic energy",EK," Potential energy",potE,
1675 & " Total energy",totE," Maximum acceleration ",
1678 write (iout,*) "Initial coordinates"
1680 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1681 & (c(j,i+nres),j=1,3)
1683 write (iout,*) "Initial dC"
1685 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1686 & (dc(j,i+nres),j=1,3)
1688 write (iout,*) "Initial velocities"
1689 write (iout,"(13x,' backbone ',23x,' side chain')")
1691 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1692 & (d_t(j,i+nres),j=1,3)
1694 write (iout,*) "Initial accelerations"
1696 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1697 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1698 & (d_a(j,i+nres),j=1,3)
1704 d_t_old(j,i)=d_t(j,i)
1705 d_a_old(j,i)=d_a(j,i)
1707 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1716 call etotal_short(energia_short)
1717 if (large) call enerprint(potEcomp)
1720 t_eshort=t_eshort+MPI_Wtime()-tt0
1722 t_eshort=t_eshort+tcpu()-tt0
1727 if(.not.out1file .and. large) then
1728 write (iout,*) "energia_long",energia_long(0),
1729 & " energia_short",energia_short(0),
1730 & " total",energia_long(0)+energia_short(0)
1731 write (iout,*) "Initial fast-force accelerations"
1733 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1734 & (d_a(j,i+nres),j=1,3)
1737 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1740 d_a_short(j,i)=d_a(j,i)
1749 call etotal_long(energia_long)
1750 if (large) call enerprint(potEcomp)
1753 t_elong=t_elong+MPI_Wtime()-tt0
1755 t_elong=t_elong+tcpu()-tt0
1760 if(.not.out1file .and. large) then
1761 write (iout,*) "energia_long",energia_long(0)
1762 write (iout,*) "Initial slow-force accelerations"
1764 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1765 & (d_a(j,i+nres),j=1,3)
1769 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1771 t_enegrad=t_enegrad+tcpu()-tt0
1776 c-----------------------------------------------------------
1777 subroutine random_vel
1778 implicit real*8 (a-h,o-z)
1779 include 'DIMENSIONS'
1780 include 'COMMON.CONTROL'
1781 include 'COMMON.VAR'
1784 include 'COMMON.LANGEVIN'
1786 include 'COMMON.LANGEVIN.lang0'
1788 include 'COMMON.CHAIN'
1789 include 'COMMON.DERIV'
1790 include 'COMMON.GEO'
1791 include 'COMMON.LOCAL'
1792 include 'COMMON.INTERACT'
1793 include 'COMMON.IOUNITS'
1794 include 'COMMON.NAMES'
1795 include 'COMMON.TIME1'
1796 double precision xv,sigv,lowb,highb
1797 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1798 c First generate velocities in the eigenspace of the G matrix
1799 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1806 sigv=dsqrt((Rb*t_bath)/geigen(i))
1809 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1810 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1811 c & " d_t_work_new",d_t_work_new(ii)
1820 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1823 c write (iout,*) "Ek from eigenvectors",Ek1
1825 c Transform velocities to UNRES coordinate space
1831 d_t_work(ind)=d_t_work(ind)
1832 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1834 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1838 c Transfer to the d_t vector
1840 d_t(j,0)=d_t_work(j)
1846 d_t(j,i)=d_t_work(ind)
1850 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1853 d_t(j,i+nres)=d_t_work(ind)
1858 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1859 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1864 c-----------------------------------------------------------
1865 subroutine sd_verlet_p_setup
1866 c Sets up the parameters of stochastic Verlet algorithm
1867 implicit real*8 (a-h,o-z)
1868 include 'DIMENSIONS'
1872 include 'COMMON.CONTROL'
1873 include 'COMMON.VAR'
1876 include 'COMMON.LANGEVIN'
1878 include 'COMMON.LANGEVIN.lang0'
1880 include 'COMMON.CHAIN'
1881 include 'COMMON.DERIV'
1882 include 'COMMON.GEO'
1883 include 'COMMON.LOCAL'
1884 include 'COMMON.INTERACT'
1885 include 'COMMON.IOUNITS'
1886 include 'COMMON.NAMES'
1887 include 'COMMON.TIME1'
1888 double precision emgdt(MAXRES6),
1889 & pterm,vterm,rho,rhoc,vsig,
1890 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1891 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1892 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1893 logical lprn /.false./
1894 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1895 double precision ktm
1902 c AL 8/17/04 Code adapted from tinker
1904 c Get the frictional and random terms for stochastic dynamics in the
1905 c eigenspace of mass-scaled UNRES friction matrix
1908 gdt = fricgam(i) * d_time
1910 c Stochastic dynamics reduces to simple MD for zero friction
1912 if (gdt .le. zero) then
1913 pfric_vec(i) = 1.0d0
1914 vfric_vec(i) = d_time
1915 afric_vec(i) = 0.5d0 * d_time * d_time
1916 prand_vec(i) = 0.0d0
1917 vrand_vec1(i) = 0.0d0
1918 vrand_vec2(i) = 0.0d0
1920 c Analytical expressions when friction coefficient is large
1923 if (gdt .ge. gdt_radius) then
1926 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1927 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1928 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1929 vterm = 1.0d0 - egdt**2
1930 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1932 c Use series expansions when friction coefficient is small
1943 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1944 & - gdt5/120.0d0 + gdt6/720.0d0
1945 & - gdt7/5040.0d0 + gdt8/40320.0d0
1946 & - gdt9/362880.0d0) / fricgam(i)**2
1947 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1948 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1949 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1950 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1951 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1952 & + 127.0d0*gdt9/90720.0d0
1953 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1954 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1955 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1956 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1957 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1958 & - 17.0d0*gdt2/1280.0d0
1959 & + 17.0d0*gdt3/6144.0d0
1960 & + 40967.0d0*gdt4/34406400.0d0
1961 & - 57203.0d0*gdt5/275251200.0d0
1962 & - 1429487.0d0*gdt6/13212057600.0d0)
1965 c Compute the scaling factors of random terms for the nonzero friction case
1967 ktm = 0.5d0*d_time/fricgam(i)
1968 psig = dsqrt(ktm*pterm) / fricgam(i)
1969 vsig = dsqrt(ktm*vterm)
1970 rhoc = dsqrt(1.0d0 - rho*rho)
1972 vrand_vec1(i) = vsig * rho
1973 vrand_vec2(i) = vsig * rhoc
1978 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1981 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1982 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1986 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1989 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1990 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1991 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1992 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1993 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1994 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1997 t_sdsetup=t_sdsetup+MPI_Wtime()
1999 t_sdsetup=t_sdsetup+tcpu()-tt0
2003 c-------------------------------------------------------------
2004 subroutine eigtransf1(n,ndim,ab,d,c)
2007 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2013 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2019 c-------------------------------------------------------------
2020 subroutine eigtransf(n,ndim,a,b,d,c)
2023 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2029 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2035 c-------------------------------------------------------------
2036 subroutine sd_verlet1
2037 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2038 implicit real*8 (a-h,o-z)
2039 include 'DIMENSIONS'
2040 include 'COMMON.CONTROL'
2041 include 'COMMON.VAR'
2044 include 'COMMON.LANGEVIN'
2046 include 'COMMON.LANGEVIN.lang0'
2048 include 'COMMON.CHAIN'
2049 include 'COMMON.DERIV'
2050 include 'COMMON.GEO'
2051 include 'COMMON.LOCAL'
2052 include 'COMMON.INTERACT'
2053 include 'COMMON.IOUNITS'
2054 include 'COMMON.NAMES'
2055 double precision stochforcvec(MAXRES6)
2056 common /stochcalc/ stochforcvec
2057 logical lprn /.false./
2059 c write (iout,*) "dc_old"
2061 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2062 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2065 dc_work(j)=dc_old(j,0)
2066 d_t_work(j)=d_t_old(j,0)
2067 d_a_work(j)=d_a_old(j,0)
2072 dc_work(ind+j)=dc_old(j,i)
2073 d_t_work(ind+j)=d_t_old(j,i)
2074 d_a_work(ind+j)=d_a_old(j,i)
2079 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2081 dc_work(ind+j)=dc_old(j,i+nres)
2082 d_t_work(ind+j)=d_t_old(j,i+nres)
2083 d_a_work(ind+j)=d_a_old(j,i+nres)
2091 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2095 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2096 & vfric_mat(i,j),afric_mat(i,j),
2097 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2105 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2106 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2107 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2108 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2110 d_t_work_new(i)=ddt1+0.5d0*ddt2
2111 d_t_work(i)=ddt1+ddt2
2116 d_t(j,0)=d_t_work(j)
2121 dc(j,i)=dc_work(ind+j)
2122 d_t(j,i)=d_t_work(ind+j)
2127 if (itype(i).ne.10) then
2130 dc(j,inres)=dc_work(ind+j)
2131 d_t(j,inres)=d_t_work(ind+j)
2138 c--------------------------------------------------------------------------
2139 subroutine sd_verlet2
2140 c Calculating the adjusted velocities for accelerations
2141 implicit real*8 (a-h,o-z)
2142 include 'DIMENSIONS'
2143 include 'COMMON.CONTROL'
2144 include 'COMMON.VAR'
2147 include 'COMMON.LANGEVIN'
2149 include 'COMMON.LANGEVIN.lang0'
2151 include 'COMMON.CHAIN'
2152 include 'COMMON.DERIV'
2153 include 'COMMON.GEO'
2154 include 'COMMON.LOCAL'
2155 include 'COMMON.INTERACT'
2156 include 'COMMON.IOUNITS'
2157 include 'COMMON.NAMES'
2158 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2159 common /stochcalc/ stochforcvec
2161 c Compute the stochastic forces which contribute to velocity change
2163 call stochastic_force(stochforcvecV)
2170 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2171 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2172 & vrand_mat2(i,j)*stochforcvecV(j)
2174 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2178 d_t(j,0)=d_t_work(j)
2183 d_t(j,i)=d_t_work(ind+j)
2188 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2191 d_t(j,inres)=d_t_work(ind+j)
2198 c-----------------------------------------------------------
2199 subroutine sd_verlet_ciccotti_setup
2200 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2202 implicit real*8 (a-h,o-z)
2203 include 'DIMENSIONS'
2207 include 'COMMON.CONTROL'
2208 include 'COMMON.VAR'
2211 include 'COMMON.LANGEVIN'
2213 include 'COMMON.LANGEVIN.lang0'
2215 include 'COMMON.CHAIN'
2216 include 'COMMON.DERIV'
2217 include 'COMMON.GEO'
2218 include 'COMMON.LOCAL'
2219 include 'COMMON.INTERACT'
2220 include 'COMMON.IOUNITS'
2221 include 'COMMON.NAMES'
2222 include 'COMMON.TIME1'
2223 double precision emgdt(MAXRES6),
2224 & pterm,vterm,rho,rhoc,vsig,
2225 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2226 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2227 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2228 logical lprn /.false./
2229 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2230 double precision ktm
2237 c AL 8/17/04 Code adapted from tinker
2239 c Get the frictional and random terms for stochastic dynamics in the
2240 c eigenspace of mass-scaled UNRES friction matrix
2243 write (iout,*) "i",i," fricgam",fricgam(i)
2244 gdt = fricgam(i) * d_time
2246 c Stochastic dynamics reduces to simple MD for zero friction
2248 if (gdt .le. zero) then
2249 pfric_vec(i) = 1.0d0
2250 vfric_vec(i) = d_time
2251 afric_vec(i) = 0.5d0*d_time*d_time
2252 prand_vec(i) = afric_vec(i)
2253 vrand_vec2(i) = vfric_vec(i)
2255 c Analytical expressions when friction coefficient is large
2260 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2261 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2262 prand_vec(i) = afric_vec(i)
2263 vrand_vec2(i) = vfric_vec(i)
2265 c Compute the scaling factors of random terms for the nonzero friction case
2267 c ktm = 0.5d0*d_time/fricgam(i)
2268 c psig = dsqrt(ktm*pterm) / fricgam(i)
2269 c vsig = dsqrt(ktm*vterm)
2270 c prand_vec(i) = psig*afric_vec(i)
2271 c vrand_vec2(i) = vsig*vfric_vec(i)
2276 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2279 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2280 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2284 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2286 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2287 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2288 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2289 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2290 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2292 t_sdsetup=t_sdsetup+MPI_Wtime()
2294 t_sdsetup=t_sdsetup+tcpu()-tt0
2298 c-------------------------------------------------------------
2299 subroutine sd_verlet1_ciccotti
2300 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2301 implicit real*8 (a-h,o-z)
2302 include 'DIMENSIONS'
2306 include 'COMMON.CONTROL'
2307 include 'COMMON.VAR'
2310 include 'COMMON.LANGEVIN'
2312 include 'COMMON.LANGEVIN.lang0'
2314 include 'COMMON.CHAIN'
2315 include 'COMMON.DERIV'
2316 include 'COMMON.GEO'
2317 include 'COMMON.LOCAL'
2318 include 'COMMON.INTERACT'
2319 include 'COMMON.IOUNITS'
2320 include 'COMMON.NAMES'
2321 double precision stochforcvec(MAXRES6)
2322 common /stochcalc/ stochforcvec
2323 logical lprn /.false./
2325 c write (iout,*) "dc_old"
2327 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2328 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2331 dc_work(j)=dc_old(j,0)
2332 d_t_work(j)=d_t_old(j,0)
2333 d_a_work(j)=d_a_old(j,0)
2338 dc_work(ind+j)=dc_old(j,i)
2339 d_t_work(ind+j)=d_t_old(j,i)
2340 d_a_work(ind+j)=d_a_old(j,i)
2345 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2347 dc_work(ind+j)=dc_old(j,i+nres)
2348 d_t_work(ind+j)=d_t_old(j,i+nres)
2349 d_a_work(ind+j)=d_a_old(j,i+nres)
2358 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2362 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2363 & vfric_mat(i,j),afric_mat(i,j),
2364 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2372 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2373 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2374 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2375 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2377 d_t_work_new(i)=ddt1+0.5d0*ddt2
2378 d_t_work(i)=ddt1+ddt2
2383 d_t(j,0)=d_t_work(j)
2388 dc(j,i)=dc_work(ind+j)
2389 d_t(j,i)=d_t_work(ind+j)
2394 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2397 dc(j,inres)=dc_work(ind+j)
2398 d_t(j,inres)=d_t_work(ind+j)
2405 c--------------------------------------------------------------------------
2406 subroutine sd_verlet2_ciccotti
2407 c Calculating the adjusted velocities for accelerations
2408 implicit real*8 (a-h,o-z)
2409 include 'DIMENSIONS'
2410 include 'COMMON.CONTROL'
2411 include 'COMMON.VAR'
2414 include 'COMMON.LANGEVIN'
2416 include 'COMMON.LANGEVIN.lang0'
2418 include 'COMMON.CHAIN'
2419 include 'COMMON.DERIV'
2420 include 'COMMON.GEO'
2421 include 'COMMON.LOCAL'
2422 include 'COMMON.INTERACT'
2423 include 'COMMON.IOUNITS'
2424 include 'COMMON.NAMES'
2425 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2426 common /stochcalc/ stochforcvec
2428 c Compute the stochastic forces which contribute to velocity change
2430 call stochastic_force(stochforcvecV)
2437 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2438 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2439 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2441 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2445 d_t(j,0)=d_t_work(j)
2450 d_t(j,i)=d_t_work(ind+j)
2455 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2458 d_t(j,inres)=d_t_work(ind+j)