2 c------------------------------------------------
3 c The driver for molecular dynamics subroutines
4 c------------------------------------------------
5 implicit real*8 (a-h,o-z)
11 include 'COMMON.SETUP'
12 include 'COMMON.CONTROL'
16 include 'COMMON.LANGEVIN'
18 include 'COMMON.LANGEVIN.lang0'
20 include 'COMMON.CHAIN'
21 include 'COMMON.DERIV'
23 include 'COMMON.LOCAL'
24 include 'COMMON.INTERACT'
25 include 'COMMON.IOUNITS'
26 include 'COMMON.NAMES'
27 include 'COMMON.TIME1'
28 include 'COMMON.HAIRPIN'
29 double precision cm(3),L(3),vcm(3)
31 double precision v_work(maxres6),v_transf(maxres6)
41 if (ilen(tmpdir).gt.0)
42 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
43 & //liczba(:ilen(liczba))//'.rst')
45 if (ilen(tmpdir).gt.0)
46 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
53 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
59 c Determine the inverse of the inertia matrix.
60 call setup_MD_matrices
64 t_MDsetup = MPI_Wtime()-tt0
66 t_MDsetup = tcpu()-tt0
69 c Entering the MD loop
75 if (lang.eq.2 .or. lang.eq.3) then
79 call sd_verlet_p_setup
81 call sd_verlet_ciccotti_setup
85 pfric0_mat(i,j,0)=pfric_mat(i,j)
86 afric0_mat(i,j,0)=afric_mat(i,j)
87 vfric0_mat(i,j,0)=vfric_mat(i,j)
88 prand0_mat(i,j,0)=prand_mat(i,j)
89 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
90 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
99 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
101 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
105 else if (lang.eq.1 .or. lang.eq.4) then
109 t_langsetup=MPI_Wtime()-tt0
112 t_langsetup=tcpu()-tt0
115 do itime=1,n_timestep
117 if (large.and. mod(itime,ntwe).eq.0)
118 & write (iout,*) "itime",itime
120 if (lang.gt.0 .and. surfarea .and.
121 & mod(itime,reset_fricmat).eq.0) then
122 if (lang.eq.2 .or. lang.eq.3) then
126 call sd_verlet_p_setup
128 call sd_verlet_ciccotti_setup
132 pfric0_mat(i,j,0)=pfric_mat(i,j)
133 afric0_mat(i,j,0)=afric_mat(i,j)
134 vfric0_mat(i,j,0)=vfric_mat(i,j)
135 prand0_mat(i,j,0)=prand_mat(i,j)
136 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
137 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
142 flag_stoch(i)=.false.
145 else if (lang.eq.1 .or. lang.eq.4) then
148 write (iout,'(a,i10)')
149 & "Friction matrix reset based on surface area, itime",itime
151 if (reset_vel .and. tbf .and. lang.eq.0
152 & .and. mod(itime,count_reset_vel).eq.0) then
154 write(iout,'(a,f20.2)')
155 & "Velocities reset to random values, time",totT
158 d_t_old(j,i)=d_t(j,i)
162 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
166 d_t(j,0)=d_t(j,0)-vcm(j)
169 kinetic_T=2.0d0/(dimen3*Rb)*EK
170 scalfac=dsqrt(T_bath/kinetic_T)
171 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
174 d_t_old(j,i)=scalfac*d_t(j,i)
180 c Time-reversible RESPA algorithm
181 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
182 call RESPA_step(itime)
184 c Variable time step algorithm.
185 call velverlet_step(itime)
189 call brown_step(itime)
191 print *,"Brown dynamics not here!"
193 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
199 if (mod(itime,ntwe).eq.0) call statout(itime)
212 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
215 v_work(ind)=d_t(j,i+nres)
220 write (66,'(80f10.5)')
221 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
225 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
227 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
229 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
232 if (mod(itime,ntwx).eq.0) then
233 write(iout,*) 'time=',itime
234 C call check_ecartint
236 write (tytul,'("time",f8.2)') totT
238 call hairpin(.true.,nharp,iharp)
239 call secondary2(.true.)
240 call pdbout(potE,tytul,ipdb)
245 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
246 open(irest2,file=rest2name,status='unknown')
247 write(irest2,*) totT,EK,potE,totE,t_bath
249 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
252 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
263 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
265 & 'MD calculations setup:',t_MDsetup,
266 & 'Energy & gradient evaluation:',t_enegrad,
267 & 'Stochastic MD setup:',t_langsetup,
268 & 'Stochastic MD step setup:',t_sdsetup,
270 write (iout,'(/28(1h=),a25,27(1h=))')
271 & ' End of MD calculation '
273 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
275 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
276 & " time_fricmatmult",time_fricmatmult," time_fsample ",
281 c-------------------------------------------------------------------------------
282 subroutine velverlet_step(itime)
283 c-------------------------------------------------------------------------------
284 c Perform a single velocity Verlet step; the time step can be rescaled if
285 c increments in accelerations exceed the threshold
286 c-------------------------------------------------------------------------------
287 implicit real*8 (a-h,o-z)
291 integer ierror,ierrcode
293 include 'COMMON.SETUP'
294 include 'COMMON.CONTROL'
298 include 'COMMON.LANGEVIN'
300 include 'COMMON.LANGEVIN.lang0'
302 include 'COMMON.CHAIN'
303 include 'COMMON.DERIV'
305 include 'COMMON.LOCAL'
306 include 'COMMON.INTERACT'
307 include 'COMMON.IOUNITS'
308 include 'COMMON.NAMES'
309 include 'COMMON.TIME1'
310 include 'COMMON.MUCA'
311 double precision vcm(3),incr(3)
312 double precision cm(3),L(3)
313 integer ilen,count,rstcount
316 integer maxcount_scale /20/
318 double precision stochforcvec(MAXRES6)
319 common /stochcalc/ stochforcvec
327 else if (lang.eq.2 .or. lang.eq.3) then
329 call stochastic_force(stochforcvec)
332 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
334 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
341 icount_scale=icount_scale+1
342 if (icount_scale.gt.maxcount_scale) then
344 & "ERROR: too many attempts at scaling down the time step. ",
345 & "amax=",amax,"epdrift=",epdrift,
346 & "damax=",damax,"edriftmax=",edriftmax,
350 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
354 c First step of the velocity Verlet algorithm
359 else if (lang.eq.3) then
361 call sd_verlet1_ciccotti
363 else if (lang.eq.1) then
368 c Build the chain from the newly calculated coordinates
370 if (rattle) call rattle1
372 if (large.and. mod(itime,ntwe).eq.0) then
373 write (iout,*) "Cartesian and internal coordinates: step 1"
378 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
379 & (dc(j,i+nres),j=1,3)
381 write (iout,*) "Accelerations"
383 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
384 & (d_a(j,i+nres),j=1,3)
386 write (iout,*) "Velocities, step 1"
388 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
389 & (d_t(j,i+nres),j=1,3)
398 c Calculate energy and forces
400 call etotal(potEcomp)
401 if (large.and. mod(itime,ntwe).eq.0)
402 & call enerprint(potEcomp)
405 t_etotal=t_etotal+MPI_Wtime()-tt0
407 t_etotal=t_etotal+tcpu()-tt0
410 potE=potEcomp(0)-potEcomp(20)
412 c Get the new accelerations
415 t_enegrad=t_enegrad+MPI_Wtime()-tt0
417 t_enegrad=t_enegrad+tcpu()-tt0
419 c Determine maximum acceleration and scale down the timestep if needed
421 amax=amax/(itime_scal+1)**2
422 call predict_edrift(epdrift)
423 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
424 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
426 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
428 itime_scal=itime_scal+ifac_time
429 c fac_time=dmin1(damax/amax,0.5d0)
430 fac_time=0.5d0**ifac_time
431 d_time=d_time*fac_time
432 if (lang.eq.2 .or. lang.eq.3) then
434 c write (iout,*) "Calling sd_verlet_setup: 1"
435 c Rescale the stochastic forces and recalculate or restore
436 c the matrices of tinker integrator
437 if (itime_scal.gt.maxflag_stoch) then
438 if (large) write (iout,'(a,i5,a)')
439 & "Calculate matrices for stochastic step;",
440 & " itime_scal ",itime_scal
442 call sd_verlet_p_setup
444 call sd_verlet_ciccotti_setup
446 write (iout,'(2a,i3,a,i3,1h.)')
447 & "Warning: cannot store matrices for stochastic",
448 & " integration because the index",itime_scal,
449 & " is greater than",maxflag_stoch
450 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
451 & " integration Langevin algorithm for better efficiency."
452 else if (flag_stoch(itime_scal)) then
453 if (large) write (iout,'(a,i5,a,l1)')
454 & "Restore matrices for stochastic step; itime_scal ",
455 & itime_scal," flag ",flag_stoch(itime_scal)
458 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
459 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
460 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
461 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
462 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
463 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
467 if (large) write (iout,'(2a,i5,a,l1)')
468 & "Calculate & store matrices for stochastic step;",
469 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
471 call sd_verlet_p_setup
473 call sd_verlet_ciccotti_setup
475 flag_stoch(ifac_time)=.true.
478 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
479 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
480 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
481 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
482 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
483 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
487 fac_time=1.0d0/dsqrt(fac_time)
489 stochforcvec(i)=fac_time*stochforcvec(i)
492 else if (lang.eq.1) then
493 c Rescale the accelerations due to stochastic forces
494 fac_time=1.0d0/dsqrt(fac_time)
496 d_as_work(i)=d_as_work(i)*fac_time
499 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
500 & "itime",itime," Timestep scaled down to ",
501 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
503 c Second step of the velocity Verlet algorithm
508 else if (lang.eq.3) then
510 call sd_verlet2_ciccotti
512 else if (lang.eq.1) then
517 if (rattle) call rattle2
519 if (d_time.ne.d_time0) then
522 if (lang.eq.2 .or. lang.eq.3) then
523 if (large) write (iout,'(a)')
524 & "Restore original matrices for stochastic step"
525 c write (iout,*) "Calling sd_verlet_setup: 2"
526 c Restore the matrices of tinker integrator if the time step has been restored
529 pfric_mat(i,j)=pfric0_mat(i,j,0)
530 afric_mat(i,j)=afric0_mat(i,j,0)
531 vfric_mat(i,j)=vfric0_mat(i,j,0)
532 prand_mat(i,j)=prand0_mat(i,j,0)
533 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
534 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
543 c Calculate the kinetic and the total energy and the kinetic temperature
548 c write (iout,*) "step",itime," EK",EK," EK1",EK1
550 c Couple the system to Berendsen bath if needed
551 if (tbf .and. lang.eq.0) then
554 kinetic_T=2.0d0/(dimen3*Rb)*EK
555 c Backup the coordinates, velocities, and accelerations
559 d_t_old(j,i)=d_t(j,i)
560 d_a_old(j,i)=d_a(j,i)
564 if (mod(itime,ntwe).eq.0 .and. large) then
565 write (iout,*) "Velocities, step 2"
567 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
568 & (d_t(j,i+nres),j=1,3)
574 c-------------------------------------------------------------------------------
575 subroutine RESPA_step(itime)
576 c-------------------------------------------------------------------------------
577 c Perform a single RESPA step.
578 c-------------------------------------------------------------------------------
579 implicit real*8 (a-h,o-z)
583 integer IERROR,ERRCODE
585 include 'COMMON.SETUP'
586 include 'COMMON.CONTROL'
590 include 'COMMON.LANGEVIN'
592 include 'COMMON.LANGEVIN.lang0'
594 include 'COMMON.CHAIN'
595 include 'COMMON.DERIV'
597 include 'COMMON.LOCAL'
598 include 'COMMON.INTERACT'
599 include 'COMMON.IOUNITS'
600 include 'COMMON.NAMES'
601 include 'COMMON.TIME1'
602 double precision energia_short(0:n_ene),
603 & energia_long(0:n_ene)
604 double precision cm(3),L(3),vcm(3),incr(3)
605 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
606 & d_a_old0(3,0:maxres2)
607 logical PRINT_AMTS_MSG /.false./
608 integer ilen,count,rstcount
611 integer maxcount_scale /10/
613 double precision stochforcvec(MAXRES6)
614 common /stochcalc/ stochforcvec
617 common /cipiszcze/ itt
620 if (large.and. mod(itime,ntwe).eq.0) then
621 write (iout,*) "***************** RESPA itime",itime
622 write (iout,*) "Cartesian and internal coordinates: step 0"
624 call pdbout(0.0d0,"cipiszcze",iout)
626 write (iout,*) "Accelerations from long-range forces"
628 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
629 & (d_a(j,i+nres),j=1,3)
631 write (iout,*) "Velocities, step 0"
633 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
634 & (d_t(j,i+nres),j=1,3)
639 c Perform the initial RESPA step (increment velocities)
640 c write (iout,*) "*********************** RESPA ini"
643 if (mod(itime,ntwe).eq.0 .and. large) then
644 write (iout,*) "Velocities, end"
646 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
647 & (d_t(j,i+nres),j=1,3)
651 c Compute the short-range forces
657 C 7/2/2009 commented out
659 c call etotal_short(energia_short)
662 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
665 d_a(j,i)=d_a_short(j,i)
669 if (large.and. mod(itime,ntwe).eq.0) then
670 write (iout,*) "energia_short",energia_short(0)
671 write (iout,*) "Accelerations from short-range forces"
673 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
674 & (d_a(j,i+nres),j=1,3)
679 t_enegrad=t_enegrad+MPI_Wtime()-tt0
681 t_enegrad=t_enegrad+tcpu()-tt0
686 d_t_old(j,i)=d_t(j,i)
687 d_a_old(j,i)=d_a(j,i)
690 c 6/30/08 A-MTS: attempt at increasing the split number
693 dc_old0(j,i)=dc_old(j,i)
694 d_t_old0(j,i)=d_t_old(j,i)
695 d_a_old0(j,i)=d_a_old(j,i)
698 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
699 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
706 c write (iout,*) "itime",itime," ntime_split",ntime_split
707 c Split the time step
708 d_time=d_time0/ntime_split
709 c Perform the short-range RESPA steps (velocity Verlet increments of
710 c positions and velocities using short-range forces)
711 c write (iout,*) "*********************** RESPA split"
712 do itsplit=1,ntime_split
715 else if (lang.eq.2 .or. lang.eq.3) then
717 call stochastic_force(stochforcvec)
720 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
722 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
727 c First step of the velocity Verlet algorithm
732 else if (lang.eq.3) then
734 call sd_verlet1_ciccotti
736 else if (lang.eq.1) then
741 c Build the chain from the newly calculated coordinates
743 if (rattle) call rattle1
745 if (large.and. mod(itime,ntwe).eq.0) then
746 write (iout,*) "***** ITSPLIT",itsplit
747 write (iout,*) "Cartesian and internal coordinates: step 1"
748 call pdbout(0.0d0,"cipiszcze",iout)
751 write (iout,*) "Velocities, step 1"
753 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
754 & (d_t(j,i+nres),j=1,3)
763 c Calculate energy and forces
765 call etotal_short(energia_short)
766 if (large.and. mod(itime,ntwe).eq.0)
767 & call enerprint(energia_short)
770 t_eshort=t_eshort+MPI_Wtime()-tt0
772 t_eshort=t_eshort+tcpu()-tt0
776 c Get the new accelerations
778 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
781 d_a_short(j,i)=d_a(j,i)
785 if (large.and. mod(itime,ntwe).eq.0) then
786 write (iout,*)"energia_short",energia_short(0)
787 write (iout,*) "Accelerations from short-range forces"
789 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
790 & (d_a(j,i+nres),j=1,3)
795 c Determine maximum acceleration and scale down the timestep if needed
797 amax=amax/ntime_split**2
798 call predict_edrift(epdrift)
799 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
800 & write (iout,*) "amax",amax," damax",damax,
801 & " epdrift",epdrift," epdriftmax",epdriftmax
802 c Exit loop and try with increased split number if the change of
803 c acceleration is too big
804 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
805 if (ntime_split.lt.maxtime_split) then
807 ntime_split=ntime_split*2
810 dc_old(j,i)=dc_old0(j,i)
811 d_t_old(j,i)=d_t_old0(j,i)
812 d_a_old(j,i)=d_a_old0(j,i)
815 if (PRINT_AMTS_MSG) then
816 write (iout,*) "acceleration/energy drift too large",amax,
817 & epdrift," split increased to ",ntime_split," itime",itime,
823 & "Uh-hu. Bumpy landscape. Maximum splitting number",
825 & " already reached!!! Trying to carry on!"
829 t_enegrad=t_enegrad+MPI_Wtime()-tt0
831 t_enegrad=t_enegrad+tcpu()-tt0
833 c Second step of the velocity Verlet algorithm
838 else if (lang.eq.3) then
840 call sd_verlet2_ciccotti
842 else if (lang.eq.1) then
847 if (rattle) call rattle2
848 c Backup the coordinates, velocities, and accelerations
852 d_t_old(j,i)=d_t(j,i)
853 d_a_old(j,i)=d_a(j,i)
860 c Restore the time step
862 c Compute long-range forces
869 call etotal_long(energia_long)
870 if (large.and. mod(itime,ntwe).eq.0)
871 & call enerprint(energia_long)
874 t_elong=t_elong+MPI_Wtime()-tt0
876 t_elong=t_elong+tcpu()-tt0
882 t_enegrad=t_enegrad+MPI_Wtime()-tt0
884 t_enegrad=t_enegrad+tcpu()-tt0
886 c Compute accelerations from long-range forces
888 if (large.and. mod(itime,ntwe).eq.0) then
889 write (iout,*) "energia_long",energia_long(0)
890 write (iout,*) "Cartesian and internal coordinates: step 2"
892 call pdbout(0.0d0,"cipiszcze",iout)
894 write (iout,*) "Accelerations from long-range forces"
896 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
897 & (d_a(j,i+nres),j=1,3)
899 write (iout,*) "Velocities, step 2"
901 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
902 & (d_t(j,i+nres),j=1,3)
906 c Compute the final RESPA step (increment velocities)
907 c write (iout,*) "*********************** RESPA fin"
909 c Compute the complete potential energy
911 potEcomp(i)=energia_short(i)+energia_long(i)
913 potE=potEcomp(0)-potEcomp(20)
914 c potE=energia_short(0)+energia_long(0)
916 c Calculate the kinetic and the total energy and the kinetic temperature
919 c Couple the system to Berendsen bath if needed
920 if (tbf .and. lang.eq.0) then
923 kinetic_T=2.0d0/(dimen3*Rb)*EK
924 c Backup the coordinates, velocities, and accelerations
926 if (mod(itime,ntwe).eq.0 .and. large) then
927 write (iout,*) "Velocities, end"
929 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
930 & (d_t(j,i+nres),j=1,3)
936 c---------------------------------------------------------------------
938 c First and last RESPA step (incrementing velocities using long-range
940 implicit real*8 (a-h,o-z)
942 include 'COMMON.CONTROL'
945 include 'COMMON.CHAIN'
946 include 'COMMON.DERIV'
948 include 'COMMON.LOCAL'
949 include 'COMMON.INTERACT'
950 include 'COMMON.IOUNITS'
951 include 'COMMON.NAMES'
953 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
957 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
961 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
964 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
970 c-----------------------------------------------------------------
972 c Applying velocity Verlet algorithm - step 1 to coordinates
973 implicit real*8 (a-h,o-z)
975 include 'COMMON.CONTROL'
978 include 'COMMON.CHAIN'
979 include 'COMMON.DERIV'
981 include 'COMMON.LOCAL'
982 include 'COMMON.INTERACT'
983 include 'COMMON.IOUNITS'
984 include 'COMMON.NAMES'
985 double precision adt,adt2
988 write (iout,*) "VELVERLET1 START: DC"
990 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
991 & (dc(j,i+nres),j=1,3)
995 adt=d_a_old(j,0)*d_time
997 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
998 d_t_new(j,0)=d_t_old(j,0)+adt2
999 d_t(j,0)=d_t_old(j,0)+adt
1005 adt=d_a_old(j,i)*d_time
1007 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1008 d_t_new(j,i)=d_t_old(j,i)+adt2
1009 d_t(j,i)=d_t_old(j,i)+adt
1014 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1017 adt=d_a_old(j,inres)*d_time
1019 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1020 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1021 d_t(j,inres)=d_t_old(j,inres)+adt
1026 write (iout,*) "VELVERLET1 END: DC"
1028 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1029 & (dc(j,i+nres),j=1,3)
1034 c---------------------------------------------------------------------
1036 c Step 2 of the velocity Verlet algorithm: update velocities
1037 implicit real*8 (a-h,o-z)
1038 include 'DIMENSIONS'
1039 include 'COMMON.CONTROL'
1040 include 'COMMON.VAR'
1042 include 'COMMON.CHAIN'
1043 include 'COMMON.DERIV'
1044 include 'COMMON.GEO'
1045 include 'COMMON.LOCAL'
1046 include 'COMMON.INTERACT'
1047 include 'COMMON.IOUNITS'
1048 include 'COMMON.NAMES'
1050 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1054 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1058 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1061 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1067 c-----------------------------------------------------------------
1068 subroutine sddir_precalc
1069 c Applying velocity Verlet algorithm - step 1 to coordinates
1070 implicit real*8 (a-h,o-z)
1071 include 'DIMENSIONS'
1075 include 'COMMON.CONTROL'
1076 include 'COMMON.VAR'
1079 include 'COMMON.LANGEVIN'
1081 include 'COMMON.LANGEVIN.lang0'
1083 include 'COMMON.CHAIN'
1084 include 'COMMON.DERIV'
1085 include 'COMMON.GEO'
1086 include 'COMMON.LOCAL'
1087 include 'COMMON.INTERACT'
1088 include 'COMMON.IOUNITS'
1089 include 'COMMON.NAMES'
1090 include 'COMMON.TIME1'
1091 double precision stochforcvec(MAXRES6)
1092 common /stochcalc/ stochforcvec
1094 c Compute friction and stochastic forces
1103 time_fric=time_fric+MPI_Wtime()-time00
1106 time_fric=time_fric+tcpu()-time00
1109 call stochastic_force(stochforcvec)
1111 time_stoch=time_stoch+MPI_Wtime()-time00
1113 time_stoch=time_stoch+tcpu()-time00
1116 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1117 c forces (d_as_work)
1119 call ginv_mult(fric_work, d_af_work)
1120 call ginv_mult(stochforcvec, d_as_work)
1123 c---------------------------------------------------------------------
1124 subroutine sddir_verlet1
1125 c Applying velocity Verlet algorithm - step 1 to velocities
1126 implicit real*8 (a-h,o-z)
1127 include 'DIMENSIONS'
1128 include 'COMMON.CONTROL'
1129 include 'COMMON.VAR'
1132 include 'COMMON.LANGEVIN'
1134 include 'COMMON.LANGEVIN.lang0'
1136 include 'COMMON.CHAIN'
1137 include 'COMMON.DERIV'
1138 include 'COMMON.GEO'
1139 include 'COMMON.LOCAL'
1140 include 'COMMON.INTERACT'
1141 include 'COMMON.IOUNITS'
1142 include 'COMMON.NAMES'
1143 c Revised 3/31/05 AL: correlation between random contributions to
1144 c position and velocity increments included.
1145 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1146 double precision adt,adt2
1148 c Add the contribution from BOTH friction and stochastic force to the
1149 c coordinates, but ONLY the contribution from the friction forces to velocities
1152 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1153 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1154 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1155 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1156 d_t(j,0)=d_t_old(j,0)+adt
1161 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1162 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1163 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1164 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1165 d_t(j,i)=d_t_old(j,i)+adt
1170 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1173 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1174 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1175 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1176 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1177 d_t(j,inres)=d_t_old(j,inres)+adt
1184 c---------------------------------------------------------------------
1185 subroutine sddir_verlet2
1186 c Calculating the adjusted velocities for accelerations
1187 implicit real*8 (a-h,o-z)
1188 include 'DIMENSIONS'
1189 include 'COMMON.CONTROL'
1190 include 'COMMON.VAR'
1193 include 'COMMON.LANGEVIN'
1195 include 'COMMON.LANGEVIN.lang0'
1197 include 'COMMON.CHAIN'
1198 include 'COMMON.DERIV'
1199 include 'COMMON.GEO'
1200 include 'COMMON.LOCAL'
1201 include 'COMMON.INTERACT'
1202 include 'COMMON.IOUNITS'
1203 include 'COMMON.NAMES'
1204 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1205 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1206 c Revised 3/31/05 AL: correlation between random contributions to
1207 c position and velocity increments included.
1208 c The correlation coefficients are calculated at low-friction limit.
1209 c Also, friction forces are now not calculated with new velocities.
1211 c call friction_force
1212 call stochastic_force(stochforcvec)
1214 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1215 c forces (d_as_work)
1217 call ginv_mult(stochforcvec, d_as_work1)
1223 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1224 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1229 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1230 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1235 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1238 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1239 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1240 & +cos60*d_as_work1(ind+j))*d_time
1247 c---------------------------------------------------------------------
1248 subroutine max_accel
1250 c Find the maximum difference in the accelerations of the the sites
1251 c at the beginning and the end of the time step.
1253 implicit real*8 (a-h,o-z)
1254 include 'DIMENSIONS'
1255 include 'COMMON.CONTROL'
1256 include 'COMMON.VAR'
1258 include 'COMMON.CHAIN'
1259 include 'COMMON.DERIV'
1260 include 'COMMON.GEO'
1261 include 'COMMON.LOCAL'
1262 include 'COMMON.INTERACT'
1263 include 'COMMON.IOUNITS'
1264 double precision aux(3),accel(3),accel_old(3),dacc
1266 c aux(j)=d_a(j,0)-d_a_old(j,0)
1267 accel_old(j)=d_a_old(j,0)
1274 c 7/3/08 changed to asymmetric difference
1276 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1277 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1278 accel(j)=accel(j)+0.5d0*d_a(j,i)
1279 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1280 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1281 dacc=dabs(accel(j)-accel_old(j))
1282 c write (iout,*) i,dacc
1283 if (dacc.gt.amax) amax=dacc
1291 accel_old(j)=d_a_old(j,0)
1296 accel_old(j)=accel_old(j)+d_a_old(j,1)
1297 accel(j)=accel(j)+d_a(j,1)
1301 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1303 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1304 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1305 accel(j)=accel(j)+d_a(j,i+nres)
1309 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1310 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1311 dacc=dabs(accel(j)-accel_old(j))
1312 c write (iout,*) "side-chain",i,dacc
1313 if (dacc.gt.amax) amax=dacc
1317 accel_old(j)=accel_old(j)+d_a_old(j,i)
1318 accel(j)=accel(j)+d_a(j,i)
1319 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1324 c---------------------------------------------------------------------
1325 subroutine predict_edrift(epdrift)
1327 c Predict the drift of the potential energy
1329 implicit real*8 (a-h,o-z)
1330 include 'DIMENSIONS'
1331 include 'COMMON.CONTROL'
1332 include 'COMMON.VAR'
1334 include 'COMMON.CHAIN'
1335 include 'COMMON.DERIV'
1336 include 'COMMON.GEO'
1337 include 'COMMON.LOCAL'
1338 include 'COMMON.INTERACT'
1339 include 'COMMON.IOUNITS'
1340 include 'COMMON.MUCA'
1341 double precision epdrift,epdriftij
1342 c Drift of the potential energy
1348 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1349 if (lmuca) epdriftij=epdriftij*factor
1350 c write (iout,*) "back",i,j,epdriftij
1351 if (epdriftij.gt.epdrift) epdrift=epdriftij
1355 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1358 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1359 if (lmuca) epdriftij=epdriftij*factor
1360 c write (iout,*) "side",i,j,epdriftij
1361 if (epdriftij.gt.epdrift) epdrift=epdriftij
1365 epdrift=0.5d0*epdrift*d_time*d_time
1366 c write (iout,*) "epdrift",epdrift
1369 c-----------------------------------------------------------------------
1370 subroutine verlet_bath
1372 c Coupling to the thermostat by using the Berendsen algorithm
1374 implicit real*8 (a-h,o-z)
1375 include 'DIMENSIONS'
1376 include 'COMMON.CONTROL'
1377 include 'COMMON.VAR'
1379 include 'COMMON.CHAIN'
1380 include 'COMMON.DERIV'
1381 include 'COMMON.GEO'
1382 include 'COMMON.LOCAL'
1383 include 'COMMON.INTERACT'
1384 include 'COMMON.IOUNITS'
1385 include 'COMMON.NAMES'
1386 double precision T_half,fact
1388 T_half=2.0d0/(dimen3*Rb)*EK
1389 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1390 c write(iout,*) "T_half", T_half
1391 c write(iout,*) "EK", EK
1392 c write(iout,*) "fact", fact
1394 d_t(j,0)=fact*d_t(j,0)
1398 d_t(j,i)=fact*d_t(j,i)
1402 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1405 d_t(j,inres)=fact*d_t(j,inres)
1411 c---------------------------------------------------------
1413 c Set up the initial conditions of a MD simulation
1414 implicit real*8 (a-h,o-z)
1415 include 'DIMENSIONS'
1419 integer IERROR,ERRCODE
1421 include 'COMMON.SETUP'
1422 include 'COMMON.CONTROL'
1423 include 'COMMON.VAR'
1426 include 'COMMON.LANGEVIN'
1428 include 'COMMON.LANGEVIN.lang0'
1430 include 'COMMON.CHAIN'
1431 include 'COMMON.DERIV'
1432 include 'COMMON.GEO'
1433 include 'COMMON.LOCAL'
1434 include 'COMMON.INTERACT'
1435 include 'COMMON.IOUNITS'
1436 include 'COMMON.NAMES'
1437 include 'COMMON.REMD'
1438 real*8 energia_long(0:n_ene),
1439 & energia_short(0:n_ene),vcm(3),incr(3)
1440 double precision cm(3),L(3),xv,sigv,lowb,highb
1441 double precision varia(maxvar)
1449 c write(iout,*) "d_time", d_time
1450 c Compute the standard deviations of stochastic forces for Langevin dynamics
1451 c if the friction coefficients do not depend on surface area
1452 if (lang.gt.0 .and. .not.surfarea) then
1454 stdforcp(i)=stdfp*dsqrt(gamp)
1457 stdforcsc(i)=stdfsc(iabs(itype(i)))
1458 & *dsqrt(gamsc(iabs(itype(i))))
1461 c Open the pdb file for snapshotshots
1464 if (ilen(tmpdir).gt.0)
1465 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1466 & liczba(:ilen(liczba))//".pdb")
1468 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1472 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1473 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1474 & liczba(:ilen(liczba))//".x")
1475 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1478 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1479 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1480 & liczba(:ilen(liczba))//".cx")
1481 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1487 if (ilen(tmpdir).gt.0)
1488 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1489 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1491 if (ilen(tmpdir).gt.0)
1492 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1493 cartname=prefix(:ilen(prefix))//"_MD.cx"
1497 write (qstr,'(256(1h ))')
1500 iq = qinfrag(i,iset)*10
1501 iw = wfrag(i,iset)/100
1503 if(me.eq.king.or..not.out1file)
1504 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1505 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1510 iq = qinpair(i,iset)*10
1511 iw = wpair(i,iset)/100
1513 if(me.eq.king.or..not.out1file)
1514 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1515 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1519 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1521 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1523 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1525 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1529 if (restart1file) then
1531 & inquire(file=mremd_rst_name,exist=file_exist)
1533 write (*,*) me," Before broadcast: file_exist",file_exist
1534 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1536 write (*,*) me," After broadcast: file_exist",file_exist
1537 c inquire(file=mremd_rst_name,exist=file_exist)
1539 if(me.eq.king.or..not.out1file)
1540 & write(iout,*) "Initial state read by master and distributed"
1542 if (ilen(tmpdir).gt.0)
1543 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1544 & //liczba(:ilen(liczba))//'.rst')
1545 inquire(file=rest2name,exist=file_exist)
1548 if(.not.restart1file) then
1549 if(me.eq.king.or..not.out1file)
1550 & write(iout,*) "Initial state will be read from file ",
1551 & rest2name(:ilen(rest2name))
1554 call rescale_weights(t_bath)
1556 if(me.eq.king.or..not.out1file)then
1557 if (restart1file) then
1558 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1561 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1564 write(iout,*) "Initial velocities randomly generated"
1570 c Generate initial velocities
1571 if(me.eq.king.or..not.out1file)
1572 & write(iout,*) "Initial velocities randomly generated"
1576 c rest2name = prefix(:ilen(prefix))//'.rst'
1577 if(me.eq.king.or..not.out1file)then
1578 write (iout,*) "Initial velocities"
1580 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1581 & (d_t(j,i+nres),j=1,3)
1583 c Zeroing the total angular momentum of the system
1584 write(iout,*) "Calling the zero-angular
1585 & momentum subroutine"
1588 c Getting the potential energy and forces and velocities and accelerations
1590 c write (iout,*) "velocity of the center of the mass:"
1591 c write (iout,*) (vcm(j),j=1,3)
1593 d_t(j,0)=d_t(j,0)-vcm(j)
1595 c Removing the velocity of the center of mass
1597 if(me.eq.king.or..not.out1file)then
1598 write (iout,*) "vcm right after adjustment:"
1599 write (iout,*) (vcm(j),j=1,3)
1603 if(iranconf.ne.0) then
1605 print *, 'Calling OVERLAP_SC'
1606 call overlap_sc(fail)
1610 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1611 print *,'SC_move',nft_sc,etot
1612 if(me.eq.king.or..not.out1file)
1613 & write(iout,*) 'SC_move',nft_sc,etot
1617 print *, 'Calling MINIM_DC'
1618 call minim_dc(etot,iretcode,nfun)
1620 call geom_to_var(nvar,varia)
1621 print *,'Calling MINIMIZE.'
1622 call minimize(etot,varia,iretcode,nfun)
1623 call var_to_geom(nvar,varia)
1625 if(me.eq.king.or..not.out1file)
1626 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1629 call chainbuild_cart
1634 kinetic_T=2.0d0/(dimen3*Rb)*EK
1635 if(me.eq.king.or..not.out1file)then
1645 call etotal(potEcomp)
1646 if (large) call enerprint(potEcomp)
1649 t_etotal=t_etotal+MPI_Wtime()-tt0
1651 t_etotal=t_etotal+tcpu()-tt0
1658 if (amax*d_time .gt. dvmax) then
1659 d_time=d_time*dvmax/amax
1660 if(me.eq.king.or..not.out1file) write (iout,*)
1661 & "Time step reduced to",d_time,
1662 & " because of too large initial acceleration."
1664 if(me.eq.king.or..not.out1file)then
1665 write(iout,*) "Potential energy and its components"
1666 call enerprint(potEcomp)
1667 c write(iout,*) (potEcomp(i),i=0,n_ene)
1669 potE=potEcomp(0)-potEcomp(20)
1672 if (ntwe.ne.0) call statout(itime)
1673 if(me.eq.king.or..not.out1file)
1674 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1675 & " Kinetic energy",EK," Potential energy",potE,
1676 & " Total energy",totE," Maximum acceleration ",
1679 write (iout,*) "Initial coordinates"
1681 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1682 & (c(j,i+nres),j=1,3)
1684 write (iout,*) "Initial dC"
1686 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1687 & (dc(j,i+nres),j=1,3)
1689 write (iout,*) "Initial velocities"
1690 write (iout,"(13x,' backbone ',23x,' side chain')")
1692 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1693 & (d_t(j,i+nres),j=1,3)
1695 write (iout,*) "Initial accelerations"
1697 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1698 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1699 & (d_a(j,i+nres),j=1,3)
1705 d_t_old(j,i)=d_t(j,i)
1706 d_a_old(j,i)=d_a(j,i)
1708 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1717 call etotal_short(energia_short)
1718 if (large) call enerprint(potEcomp)
1721 t_eshort=t_eshort+MPI_Wtime()-tt0
1723 t_eshort=t_eshort+tcpu()-tt0
1728 if(.not.out1file .and. large) then
1729 write (iout,*) "energia_long",energia_long(0),
1730 & " energia_short",energia_short(0),
1731 & " total",energia_long(0)+energia_short(0)
1732 write (iout,*) "Initial fast-force accelerations"
1734 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1735 & (d_a(j,i+nres),j=1,3)
1738 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1741 d_a_short(j,i)=d_a(j,i)
1750 call etotal_long(energia_long)
1751 if (large) call enerprint(potEcomp)
1754 t_elong=t_elong+MPI_Wtime()-tt0
1756 t_elong=t_elong+tcpu()-tt0
1761 if(.not.out1file .and. large) then
1762 write (iout,*) "energia_long",energia_long(0)
1763 write (iout,*) "Initial slow-force accelerations"
1765 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1766 & (d_a(j,i+nres),j=1,3)
1770 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1772 t_enegrad=t_enegrad+tcpu()-tt0
1777 c-----------------------------------------------------------
1778 subroutine random_vel
1779 implicit real*8 (a-h,o-z)
1780 include 'DIMENSIONS'
1781 include 'COMMON.CONTROL'
1782 include 'COMMON.VAR'
1785 include 'COMMON.LANGEVIN'
1787 include 'COMMON.LANGEVIN.lang0'
1789 include 'COMMON.CHAIN'
1790 include 'COMMON.DERIV'
1791 include 'COMMON.GEO'
1792 include 'COMMON.LOCAL'
1793 include 'COMMON.INTERACT'
1794 include 'COMMON.IOUNITS'
1795 include 'COMMON.NAMES'
1796 include 'COMMON.TIME1'
1797 double precision xv,sigv,lowb,highb
1798 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1799 c First generate velocities in the eigenspace of the G matrix
1800 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1807 sigv=dsqrt((Rb*t_bath)/geigen(i))
1810 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1811 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1812 c & " d_t_work_new",d_t_work_new(ii)
1821 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1824 c write (iout,*) "Ek from eigenvectors",Ek1
1826 c Transform velocities to UNRES coordinate space
1832 d_t_work(ind)=d_t_work(ind)
1833 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1835 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1839 c Transfer to the d_t vector
1841 d_t(j,0)=d_t_work(j)
1847 d_t(j,i)=d_t_work(ind)
1851 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1854 d_t(j,i+nres)=d_t_work(ind)
1859 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1860 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1865 c-----------------------------------------------------------
1866 subroutine sd_verlet_p_setup
1867 c Sets up the parameters of stochastic Verlet algorithm
1868 implicit real*8 (a-h,o-z)
1869 include 'DIMENSIONS'
1873 include 'COMMON.CONTROL'
1874 include 'COMMON.VAR'
1877 include 'COMMON.LANGEVIN'
1879 include 'COMMON.LANGEVIN.lang0'
1881 include 'COMMON.CHAIN'
1882 include 'COMMON.DERIV'
1883 include 'COMMON.GEO'
1884 include 'COMMON.LOCAL'
1885 include 'COMMON.INTERACT'
1886 include 'COMMON.IOUNITS'
1887 include 'COMMON.NAMES'
1888 include 'COMMON.TIME1'
1889 double precision emgdt(MAXRES6),
1890 & pterm,vterm,rho,rhoc,vsig,
1891 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1892 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1893 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1894 logical lprn /.false./
1895 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1896 double precision ktm
1903 c AL 8/17/04 Code adapted from tinker
1905 c Get the frictional and random terms for stochastic dynamics in the
1906 c eigenspace of mass-scaled UNRES friction matrix
1909 gdt = fricgam(i) * d_time
1911 c Stochastic dynamics reduces to simple MD for zero friction
1913 if (gdt .le. zero) then
1914 pfric_vec(i) = 1.0d0
1915 vfric_vec(i) = d_time
1916 afric_vec(i) = 0.5d0 * d_time * d_time
1917 prand_vec(i) = 0.0d0
1918 vrand_vec1(i) = 0.0d0
1919 vrand_vec2(i) = 0.0d0
1921 c Analytical expressions when friction coefficient is large
1924 if (gdt .ge. gdt_radius) then
1927 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1928 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1929 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1930 vterm = 1.0d0 - egdt**2
1931 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1933 c Use series expansions when friction coefficient is small
1944 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1945 & - gdt5/120.0d0 + gdt6/720.0d0
1946 & - gdt7/5040.0d0 + gdt8/40320.0d0
1947 & - gdt9/362880.0d0) / fricgam(i)**2
1948 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1949 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1950 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1951 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1952 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1953 & + 127.0d0*gdt9/90720.0d0
1954 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1955 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1956 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1957 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1958 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1959 & - 17.0d0*gdt2/1280.0d0
1960 & + 17.0d0*gdt3/6144.0d0
1961 & + 40967.0d0*gdt4/34406400.0d0
1962 & - 57203.0d0*gdt5/275251200.0d0
1963 & - 1429487.0d0*gdt6/13212057600.0d0)
1966 c Compute the scaling factors of random terms for the nonzero friction case
1968 ktm = 0.5d0*d_time/fricgam(i)
1969 psig = dsqrt(ktm*pterm) / fricgam(i)
1970 vsig = dsqrt(ktm*vterm)
1971 rhoc = dsqrt(1.0d0 - rho*rho)
1973 vrand_vec1(i) = vsig * rho
1974 vrand_vec2(i) = vsig * rhoc
1979 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1982 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1983 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1987 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1990 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1991 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1992 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1993 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1994 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1995 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1998 t_sdsetup=t_sdsetup+MPI_Wtime()
2000 t_sdsetup=t_sdsetup+tcpu()-tt0
2004 c-------------------------------------------------------------
2005 subroutine eigtransf1(n,ndim,ab,d,c)
2008 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2014 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2020 c-------------------------------------------------------------
2021 subroutine eigtransf(n,ndim,a,b,d,c)
2024 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2030 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2036 c-------------------------------------------------------------
2037 subroutine sd_verlet1
2038 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2039 implicit real*8 (a-h,o-z)
2040 include 'DIMENSIONS'
2041 include 'COMMON.CONTROL'
2042 include 'COMMON.VAR'
2045 include 'COMMON.LANGEVIN'
2047 include 'COMMON.LANGEVIN.lang0'
2049 include 'COMMON.CHAIN'
2050 include 'COMMON.DERIV'
2051 include 'COMMON.GEO'
2052 include 'COMMON.LOCAL'
2053 include 'COMMON.INTERACT'
2054 include 'COMMON.IOUNITS'
2055 include 'COMMON.NAMES'
2056 double precision stochforcvec(MAXRES6)
2057 common /stochcalc/ stochforcvec
2058 logical lprn /.false./
2060 c write (iout,*) "dc_old"
2062 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2063 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2066 dc_work(j)=dc_old(j,0)
2067 d_t_work(j)=d_t_old(j,0)
2068 d_a_work(j)=d_a_old(j,0)
2073 dc_work(ind+j)=dc_old(j,i)
2074 d_t_work(ind+j)=d_t_old(j,i)
2075 d_a_work(ind+j)=d_a_old(j,i)
2080 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2082 dc_work(ind+j)=dc_old(j,i+nres)
2083 d_t_work(ind+j)=d_t_old(j,i+nres)
2084 d_a_work(ind+j)=d_a_old(j,i+nres)
2092 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2096 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2097 & vfric_mat(i,j),afric_mat(i,j),
2098 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2106 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2107 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2108 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2109 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2111 d_t_work_new(i)=ddt1+0.5d0*ddt2
2112 d_t_work(i)=ddt1+ddt2
2117 d_t(j,0)=d_t_work(j)
2122 dc(j,i)=dc_work(ind+j)
2123 d_t(j,i)=d_t_work(ind+j)
2128 if (itype(i).ne.10) then
2131 dc(j,inres)=dc_work(ind+j)
2132 d_t(j,inres)=d_t_work(ind+j)
2139 c--------------------------------------------------------------------------
2140 subroutine sd_verlet2
2141 c Calculating the adjusted velocities for accelerations
2142 implicit real*8 (a-h,o-z)
2143 include 'DIMENSIONS'
2144 include 'COMMON.CONTROL'
2145 include 'COMMON.VAR'
2148 include 'COMMON.LANGEVIN'
2150 include 'COMMON.LANGEVIN.lang0'
2152 include 'COMMON.CHAIN'
2153 include 'COMMON.DERIV'
2154 include 'COMMON.GEO'
2155 include 'COMMON.LOCAL'
2156 include 'COMMON.INTERACT'
2157 include 'COMMON.IOUNITS'
2158 include 'COMMON.NAMES'
2159 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2160 common /stochcalc/ stochforcvec
2162 c Compute the stochastic forces which contribute to velocity change
2164 call stochastic_force(stochforcvecV)
2171 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2172 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2173 & vrand_mat2(i,j)*stochforcvecV(j)
2175 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2179 d_t(j,0)=d_t_work(j)
2184 d_t(j,i)=d_t_work(ind+j)
2189 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2192 d_t(j,inres)=d_t_work(ind+j)
2199 c-----------------------------------------------------------
2200 subroutine sd_verlet_ciccotti_setup
2201 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2203 implicit real*8 (a-h,o-z)
2204 include 'DIMENSIONS'
2208 include 'COMMON.CONTROL'
2209 include 'COMMON.VAR'
2212 include 'COMMON.LANGEVIN'
2214 include 'COMMON.LANGEVIN.lang0'
2216 include 'COMMON.CHAIN'
2217 include 'COMMON.DERIV'
2218 include 'COMMON.GEO'
2219 include 'COMMON.LOCAL'
2220 include 'COMMON.INTERACT'
2221 include 'COMMON.IOUNITS'
2222 include 'COMMON.NAMES'
2223 include 'COMMON.TIME1'
2224 double precision emgdt(MAXRES6),
2225 & pterm,vterm,rho,rhoc,vsig,
2226 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2227 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2228 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2229 logical lprn /.false./
2230 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2231 double precision ktm
2238 c AL 8/17/04 Code adapted from tinker
2240 c Get the frictional and random terms for stochastic dynamics in the
2241 c eigenspace of mass-scaled UNRES friction matrix
2244 write (iout,*) "i",i," fricgam",fricgam(i)
2245 gdt = fricgam(i) * d_time
2247 c Stochastic dynamics reduces to simple MD for zero friction
2249 if (gdt .le. zero) then
2250 pfric_vec(i) = 1.0d0
2251 vfric_vec(i) = d_time
2252 afric_vec(i) = 0.5d0*d_time*d_time
2253 prand_vec(i) = afric_vec(i)
2254 vrand_vec2(i) = vfric_vec(i)
2256 c Analytical expressions when friction coefficient is large
2261 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2262 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2263 prand_vec(i) = afric_vec(i)
2264 vrand_vec2(i) = vfric_vec(i)
2266 c Compute the scaling factors of random terms for the nonzero friction case
2268 c ktm = 0.5d0*d_time/fricgam(i)
2269 c psig = dsqrt(ktm*pterm) / fricgam(i)
2270 c vsig = dsqrt(ktm*vterm)
2271 c prand_vec(i) = psig*afric_vec(i)
2272 c vrand_vec2(i) = vsig*vfric_vec(i)
2277 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2280 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2281 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2285 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2287 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2288 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2289 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2290 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2291 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2293 t_sdsetup=t_sdsetup+MPI_Wtime()
2295 t_sdsetup=t_sdsetup+tcpu()-tt0
2299 c-------------------------------------------------------------
2300 subroutine sd_verlet1_ciccotti
2301 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2302 implicit real*8 (a-h,o-z)
2303 include 'DIMENSIONS'
2307 include 'COMMON.CONTROL'
2308 include 'COMMON.VAR'
2311 include 'COMMON.LANGEVIN'
2313 include 'COMMON.LANGEVIN.lang0'
2315 include 'COMMON.CHAIN'
2316 include 'COMMON.DERIV'
2317 include 'COMMON.GEO'
2318 include 'COMMON.LOCAL'
2319 include 'COMMON.INTERACT'
2320 include 'COMMON.IOUNITS'
2321 include 'COMMON.NAMES'
2322 double precision stochforcvec(MAXRES6)
2323 common /stochcalc/ stochforcvec
2324 logical lprn /.false./
2326 c write (iout,*) "dc_old"
2328 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2329 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2332 dc_work(j)=dc_old(j,0)
2333 d_t_work(j)=d_t_old(j,0)
2334 d_a_work(j)=d_a_old(j,0)
2339 dc_work(ind+j)=dc_old(j,i)
2340 d_t_work(ind+j)=d_t_old(j,i)
2341 d_a_work(ind+j)=d_a_old(j,i)
2346 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2348 dc_work(ind+j)=dc_old(j,i+nres)
2349 d_t_work(ind+j)=d_t_old(j,i+nres)
2350 d_a_work(ind+j)=d_a_old(j,i+nres)
2359 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2363 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2364 & vfric_mat(i,j),afric_mat(i,j),
2365 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2373 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2374 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2375 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2376 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2378 d_t_work_new(i)=ddt1+0.5d0*ddt2
2379 d_t_work(i)=ddt1+ddt2
2384 d_t(j,0)=d_t_work(j)
2389 dc(j,i)=dc_work(ind+j)
2390 d_t(j,i)=d_t_work(ind+j)
2395 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2398 dc(j,inres)=dc_work(ind+j)
2399 d_t(j,inres)=d_t_work(ind+j)
2406 c--------------------------------------------------------------------------
2407 subroutine sd_verlet2_ciccotti
2408 c Calculating the adjusted velocities for accelerations
2409 implicit real*8 (a-h,o-z)
2410 include 'DIMENSIONS'
2411 include 'COMMON.CONTROL'
2412 include 'COMMON.VAR'
2415 include 'COMMON.LANGEVIN'
2417 include 'COMMON.LANGEVIN.lang0'
2419 include 'COMMON.CHAIN'
2420 include 'COMMON.DERIV'
2421 include 'COMMON.GEO'
2422 include 'COMMON.LOCAL'
2423 include 'COMMON.INTERACT'
2424 include 'COMMON.IOUNITS'
2425 include 'COMMON.NAMES'
2426 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2427 common /stochcalc/ stochforcvec
2429 c Compute the stochastic forces which contribute to velocity change
2431 call stochastic_force(stochforcvecV)
2438 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2439 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2440 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2442 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2446 d_t(j,0)=d_t_work(j)
2451 d_t(j,i)=d_t_work(ind+j)
2456 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2459 d_t(j,inres)=d_t_work(ind+j)