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) then
201 C call enerprint(potEcomp)
202 C print *,itime,'AFM',Eafmforc,etot
216 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
219 v_work(ind)=d_t(j,i+nres)
224 write (66,'(80f10.5)')
225 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
229 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
231 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
233 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
236 if (mod(itime,ntwx).eq.0) then
237 write(iout,*) 'time=',itime
238 C call check_ecartint
240 write (tytul,'("time",f8.2)') totT
242 call hairpin(.true.,nharp,iharp)
243 call secondary2(.true.)
244 call pdbout(potE,tytul,ipdb)
249 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
250 open(irest2,file=rest2name,status='unknown')
251 write(irest2,*) totT,EK,potE,totE,t_bath
253 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
256 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
267 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
269 & 'MD calculations setup:',t_MDsetup,
270 & 'Energy & gradient evaluation:',t_enegrad,
271 & 'Stochastic MD setup:',t_langsetup,
272 & 'Stochastic MD step setup:',t_sdsetup,
274 write (iout,'(/28(1h=),a25,27(1h=))')
275 & ' End of MD calculation '
277 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
279 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
280 & " time_fricmatmult",time_fricmatmult," time_fsample ",
285 c-------------------------------------------------------------------------------
286 subroutine velverlet_step(itime)
287 c-------------------------------------------------------------------------------
288 c Perform a single velocity Verlet step; the time step can be rescaled if
289 c increments in accelerations exceed the threshold
290 c-------------------------------------------------------------------------------
291 implicit real*8 (a-h,o-z)
295 integer ierror,ierrcode
297 include 'COMMON.SETUP'
298 include 'COMMON.CONTROL'
302 include 'COMMON.LANGEVIN'
304 include 'COMMON.LANGEVIN.lang0'
306 include 'COMMON.CHAIN'
307 include 'COMMON.DERIV'
309 include 'COMMON.LOCAL'
310 include 'COMMON.INTERACT'
311 include 'COMMON.IOUNITS'
312 include 'COMMON.NAMES'
313 include 'COMMON.TIME1'
314 include 'COMMON.MUCA'
315 double precision vcm(3),incr(3)
316 double precision cm(3),L(3)
317 integer ilen,count,rstcount
320 integer maxcount_scale /20/
322 double precision stochforcvec(MAXRES6)
323 common /stochcalc/ stochforcvec
331 else if (lang.eq.2 .or. lang.eq.3) then
333 call stochastic_force(stochforcvec)
336 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
338 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
345 icount_scale=icount_scale+1
346 if (icount_scale.gt.maxcount_scale) then
348 & "ERROR: too many attempts at scaling down the time step. ",
349 & "amax=",amax,"epdrift=",epdrift,
350 & "damax=",damax,"edriftmax=",edriftmax,
354 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
358 c First step of the velocity Verlet algorithm
363 else if (lang.eq.3) then
365 call sd_verlet1_ciccotti
367 else if (lang.eq.1) then
372 c Build the chain from the newly calculated coordinates
374 if (rattle) call rattle1
376 if (large.and. mod(itime,ntwe).eq.0) then
377 write (iout,*) "Cartesian and internal coordinates: step 1"
382 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
383 & (dc(j,i+nres),j=1,3)
385 write (iout,*) "Accelerations"
387 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
388 & (d_a(j,i+nres),j=1,3)
390 write (iout,*) "Velocities, step 1"
392 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
393 & (d_t(j,i+nres),j=1,3)
402 c Calculate energy and forces
404 call etotal(potEcomp)
405 if (large.and. mod(itime,ntwe).eq.0)
406 & call enerprint(potEcomp)
409 t_etotal=t_etotal+MPI_Wtime()-tt0
411 t_etotal=t_etotal+tcpu()-tt0
414 potE=potEcomp(0)-potEcomp(20)
416 c Get the new accelerations
419 t_enegrad=t_enegrad+MPI_Wtime()-tt0
421 t_enegrad=t_enegrad+tcpu()-tt0
423 c Determine maximum acceleration and scale down the timestep if needed
425 amax=amax/(itime_scal+1)**2
426 call predict_edrift(epdrift)
427 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
428 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
430 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
432 itime_scal=itime_scal+ifac_time
433 c fac_time=dmin1(damax/amax,0.5d0)
434 fac_time=0.5d0**ifac_time
435 d_time=d_time*fac_time
436 if (lang.eq.2 .or. lang.eq.3) then
438 c write (iout,*) "Calling sd_verlet_setup: 1"
439 c Rescale the stochastic forces and recalculate or restore
440 c the matrices of tinker integrator
441 if (itime_scal.gt.maxflag_stoch) then
442 if (large) write (iout,'(a,i5,a)')
443 & "Calculate matrices for stochastic step;",
444 & " itime_scal ",itime_scal
446 call sd_verlet_p_setup
448 call sd_verlet_ciccotti_setup
450 write (iout,'(2a,i3,a,i3,1h.)')
451 & "Warning: cannot store matrices for stochastic",
452 & " integration because the index",itime_scal,
453 & " is greater than",maxflag_stoch
454 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
455 & " integration Langevin algorithm for better efficiency."
456 else if (flag_stoch(itime_scal)) then
457 if (large) write (iout,'(a,i5,a,l1)')
458 & "Restore matrices for stochastic step; itime_scal ",
459 & itime_scal," flag ",flag_stoch(itime_scal)
462 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
463 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
464 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
465 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
466 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
467 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
471 if (large) write (iout,'(2a,i5,a,l1)')
472 & "Calculate & store matrices for stochastic step;",
473 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
475 call sd_verlet_p_setup
477 call sd_verlet_ciccotti_setup
479 flag_stoch(ifac_time)=.true.
482 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
483 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
484 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
485 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
486 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
487 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
491 fac_time=1.0d0/dsqrt(fac_time)
493 stochforcvec(i)=fac_time*stochforcvec(i)
496 else if (lang.eq.1) then
497 c Rescale the accelerations due to stochastic forces
498 fac_time=1.0d0/dsqrt(fac_time)
500 d_as_work(i)=d_as_work(i)*fac_time
503 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
504 & "itime",itime," Timestep scaled down to ",
505 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
507 c Second step of the velocity Verlet algorithm
512 else if (lang.eq.3) then
514 call sd_verlet2_ciccotti
516 else if (lang.eq.1) then
521 if (rattle) call rattle2
524 C print *,totTafm,"TU?"
525 if (d_time.ne.d_time0) then
528 if (lang.eq.2 .or. lang.eq.3) then
529 if (large) write (iout,'(a)')
530 & "Restore original matrices for stochastic step"
531 c write (iout,*) "Calling sd_verlet_setup: 2"
532 c Restore the matrices of tinker integrator if the time step has been restored
535 pfric_mat(i,j)=pfric0_mat(i,j,0)
536 afric_mat(i,j)=afric0_mat(i,j,0)
537 vfric_mat(i,j)=vfric0_mat(i,j,0)
538 prand_mat(i,j)=prand0_mat(i,j,0)
539 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
540 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
549 c Calculate the kinetic and the total energy and the kinetic temperature
554 c write (iout,*) "step",itime," EK",EK," EK1",EK1
556 c Couple the system to Berendsen bath if needed
557 if (tbf .and. lang.eq.0) then
560 kinetic_T=2.0d0/(dimen3*Rb)*EK
561 c Backup the coordinates, velocities, and accelerations
565 d_t_old(j,i)=d_t(j,i)
566 d_a_old(j,i)=d_a(j,i)
570 if (mod(itime,ntwe).eq.0 .and. large) then
571 write (iout,*) "Velocities, step 2"
573 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
574 & (d_t(j,i+nres),j=1,3)
580 c-------------------------------------------------------------------------------
581 subroutine RESPA_step(itime)
582 c-------------------------------------------------------------------------------
583 c Perform a single RESPA step.
584 c-------------------------------------------------------------------------------
585 implicit real*8 (a-h,o-z)
589 integer IERROR,ERRCODE
591 include 'COMMON.SETUP'
592 include 'COMMON.CONTROL'
596 include 'COMMON.LANGEVIN'
598 include 'COMMON.LANGEVIN.lang0'
600 include 'COMMON.CHAIN'
601 include 'COMMON.DERIV'
603 include 'COMMON.LOCAL'
604 include 'COMMON.INTERACT'
605 include 'COMMON.IOUNITS'
606 include 'COMMON.NAMES'
607 include 'COMMON.TIME1'
608 double precision energia_short(0:n_ene),
609 & energia_long(0:n_ene)
610 double precision cm(3),L(3),vcm(3),incr(3)
611 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
612 & d_a_old0(3,0:maxres2)
613 logical PRINT_AMTS_MSG /.false./
614 integer ilen,count,rstcount
617 integer maxcount_scale /10/
619 double precision stochforcvec(MAXRES6)
620 common /stochcalc/ stochforcvec
623 common /cipiszcze/ itt
626 if (large.and. mod(itime,ntwe).eq.0) then
627 write (iout,*) "***************** RESPA itime",itime
628 write (iout,*) "Cartesian and internal coordinates: step 0"
630 call pdbout(0.0d0,"cipiszcze",iout)
632 write (iout,*) "Accelerations from long-range forces"
634 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
635 & (d_a(j,i+nres),j=1,3)
637 write (iout,*) "Velocities, step 0"
639 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
640 & (d_t(j,i+nres),j=1,3)
645 c Perform the initial RESPA step (increment velocities)
646 c write (iout,*) "*********************** RESPA ini"
649 if (mod(itime,ntwe).eq.0 .and. large) then
650 write (iout,*) "Velocities, end"
652 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
653 & (d_t(j,i+nres),j=1,3)
657 c Compute the short-range forces
663 C 7/2/2009 commented out
665 c call etotal_short(energia_short)
668 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
671 d_a(j,i)=d_a_short(j,i)
675 if (large.and. mod(itime,ntwe).eq.0) then
676 write (iout,*) "energia_short",energia_short(0)
677 write (iout,*) "Accelerations from short-range forces"
679 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
680 & (d_a(j,i+nres),j=1,3)
685 t_enegrad=t_enegrad+MPI_Wtime()-tt0
687 t_enegrad=t_enegrad+tcpu()-tt0
692 d_t_old(j,i)=d_t(j,i)
693 d_a_old(j,i)=d_a(j,i)
696 c 6/30/08 A-MTS: attempt at increasing the split number
699 dc_old0(j,i)=dc_old(j,i)
700 d_t_old0(j,i)=d_t_old(j,i)
701 d_a_old0(j,i)=d_a_old(j,i)
704 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
705 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
712 c write (iout,*) "itime",itime," ntime_split",ntime_split
713 c Split the time step
714 d_time=d_time0/ntime_split
715 c Perform the short-range RESPA steps (velocity Verlet increments of
716 c positions and velocities using short-range forces)
717 c write (iout,*) "*********************** RESPA split"
718 do itsplit=1,ntime_split
721 else if (lang.eq.2 .or. lang.eq.3) then
723 call stochastic_force(stochforcvec)
726 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
728 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
733 c First step of the velocity Verlet algorithm
738 else if (lang.eq.3) then
740 call sd_verlet1_ciccotti
742 else if (lang.eq.1) then
747 c Build the chain from the newly calculated coordinates
749 if (rattle) call rattle1
751 if (large.and. mod(itime,ntwe).eq.0) then
752 write (iout,*) "***** ITSPLIT",itsplit
753 write (iout,*) "Cartesian and internal coordinates: step 1"
754 call pdbout(0.0d0,"cipiszcze",iout)
757 write (iout,*) "Velocities, step 1"
759 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
760 & (d_t(j,i+nres),j=1,3)
769 c Calculate energy and forces
771 call etotal_short(energia_short)
772 if (large.and. mod(itime,ntwe).eq.0)
773 & call enerprint(energia_short)
776 t_eshort=t_eshort+MPI_Wtime()-tt0
778 t_eshort=t_eshort+tcpu()-tt0
782 c Get the new accelerations
784 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
787 d_a_short(j,i)=d_a(j,i)
791 if (large.and. mod(itime,ntwe).eq.0) then
792 write (iout,*)"energia_short",energia_short(0)
793 write (iout,*) "Accelerations from short-range forces"
795 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
796 & (d_a(j,i+nres),j=1,3)
801 c Determine maximum acceleration and scale down the timestep if needed
803 amax=amax/ntime_split**2
804 call predict_edrift(epdrift)
805 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
806 & write (iout,*) "amax",amax," damax",damax,
807 & " epdrift",epdrift," epdriftmax",epdriftmax
808 c Exit loop and try with increased split number if the change of
809 c acceleration is too big
810 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
811 if (ntime_split.lt.maxtime_split) then
813 ntime_split=ntime_split*2
816 dc_old(j,i)=dc_old0(j,i)
817 d_t_old(j,i)=d_t_old0(j,i)
818 d_a_old(j,i)=d_a_old0(j,i)
821 if (PRINT_AMTS_MSG) then
822 write (iout,*) "acceleration/energy drift too large",amax,
823 & epdrift," split increased to ",ntime_split," itime",itime,
829 & "Uh-hu. Bumpy landscape. Maximum splitting number",
831 & " already reached!!! Trying to carry on!"
835 t_enegrad=t_enegrad+MPI_Wtime()-tt0
837 t_enegrad=t_enegrad+tcpu()-tt0
839 c Second step of the velocity Verlet algorithm
844 else if (lang.eq.3) then
846 call sd_verlet2_ciccotti
848 else if (lang.eq.1) then
853 if (rattle) call rattle2
854 c Backup the coordinates, velocities, and accelerations
858 d_t_old(j,i)=d_t(j,i)
859 d_a_old(j,i)=d_a(j,i)
866 c Restore the time step
868 c Compute long-range forces
875 call etotal_long(energia_long)
876 if (large.and. mod(itime,ntwe).eq.0)
877 & call enerprint(energia_long)
880 t_elong=t_elong+MPI_Wtime()-tt0
882 t_elong=t_elong+tcpu()-tt0
888 t_enegrad=t_enegrad+MPI_Wtime()-tt0
890 t_enegrad=t_enegrad+tcpu()-tt0
892 c Compute accelerations from long-range forces
894 if (large.and. mod(itime,ntwe).eq.0) then
895 write (iout,*) "energia_long",energia_long(0)
896 write (iout,*) "Cartesian and internal coordinates: step 2"
898 call pdbout(0.0d0,"cipiszcze",iout)
900 write (iout,*) "Accelerations from long-range forces"
902 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
903 & (d_a(j,i+nres),j=1,3)
905 write (iout,*) "Velocities, step 2"
907 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
908 & (d_t(j,i+nres),j=1,3)
912 c Compute the final RESPA step (increment velocities)
913 c write (iout,*) "*********************** RESPA fin"
915 c Compute the complete potential energy
917 potEcomp(i)=energia_short(i)+energia_long(i)
919 potE=potEcomp(0)-potEcomp(20)
920 c potE=energia_short(0)+energia_long(0)
923 c Calculate the kinetic and the total energy and the kinetic temperature
926 c Couple the system to Berendsen bath if needed
927 if (tbf .and. lang.eq.0) then
930 kinetic_T=2.0d0/(dimen3*Rb)*EK
931 c Backup the coordinates, velocities, and accelerations
933 if (mod(itime,ntwe).eq.0 .and. large) then
934 write (iout,*) "Velocities, end"
936 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
937 & (d_t(j,i+nres),j=1,3)
943 c---------------------------------------------------------------------
945 c First and last RESPA step (incrementing velocities using long-range
947 implicit real*8 (a-h,o-z)
949 include 'COMMON.CONTROL'
952 include 'COMMON.CHAIN'
953 include 'COMMON.DERIV'
955 include 'COMMON.LOCAL'
956 include 'COMMON.INTERACT'
957 include 'COMMON.IOUNITS'
958 include 'COMMON.NAMES'
960 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
964 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
968 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
971 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
977 c-----------------------------------------------------------------
979 c Applying velocity Verlet algorithm - step 1 to coordinates
980 implicit real*8 (a-h,o-z)
982 include 'COMMON.CONTROL'
985 include 'COMMON.CHAIN'
986 include 'COMMON.DERIV'
988 include 'COMMON.LOCAL'
989 include 'COMMON.INTERACT'
990 include 'COMMON.IOUNITS'
991 include 'COMMON.NAMES'
992 double precision adt,adt2
995 write (iout,*) "VELVERLET1 START: DC"
997 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
998 & (dc(j,i+nres),j=1,3)
1002 adt=d_a_old(j,0)*d_time
1004 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1005 d_t_new(j,0)=d_t_old(j,0)+adt2
1006 d_t(j,0)=d_t_old(j,0)+adt
1012 adt=d_a_old(j,i)*d_time
1014 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1015 d_t_new(j,i)=d_t_old(j,i)+adt2
1016 d_t(j,i)=d_t_old(j,i)+adt
1021 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1024 adt=d_a_old(j,inres)*d_time
1026 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1027 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1028 d_t(j,inres)=d_t_old(j,inres)+adt
1033 write (iout,*) "VELVERLET1 END: DC"
1035 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1036 & (dc(j,i+nres),j=1,3)
1041 c---------------------------------------------------------------------
1043 c Step 2 of the velocity Verlet algorithm: update velocities
1044 implicit real*8 (a-h,o-z)
1045 include 'DIMENSIONS'
1046 include 'COMMON.CONTROL'
1047 include 'COMMON.VAR'
1049 include 'COMMON.CHAIN'
1050 include 'COMMON.DERIV'
1051 include 'COMMON.GEO'
1052 include 'COMMON.LOCAL'
1053 include 'COMMON.INTERACT'
1054 include 'COMMON.IOUNITS'
1055 include 'COMMON.NAMES'
1057 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1061 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1065 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1068 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1074 c-----------------------------------------------------------------
1075 subroutine sddir_precalc
1076 c Applying velocity Verlet algorithm - step 1 to coordinates
1077 implicit real*8 (a-h,o-z)
1078 include 'DIMENSIONS'
1082 include 'COMMON.CONTROL'
1083 include 'COMMON.VAR'
1086 include 'COMMON.LANGEVIN'
1088 include 'COMMON.LANGEVIN.lang0'
1090 include 'COMMON.CHAIN'
1091 include 'COMMON.DERIV'
1092 include 'COMMON.GEO'
1093 include 'COMMON.LOCAL'
1094 include 'COMMON.INTERACT'
1095 include 'COMMON.IOUNITS'
1096 include 'COMMON.NAMES'
1097 include 'COMMON.TIME1'
1098 double precision stochforcvec(MAXRES6)
1099 common /stochcalc/ stochforcvec
1101 c Compute friction and stochastic forces
1110 time_fric=time_fric+MPI_Wtime()-time00
1113 time_fric=time_fric+tcpu()-time00
1116 call stochastic_force(stochforcvec)
1118 time_stoch=time_stoch+MPI_Wtime()-time00
1120 time_stoch=time_stoch+tcpu()-time00
1123 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1124 c forces (d_as_work)
1126 call ginv_mult(fric_work, d_af_work)
1127 call ginv_mult(stochforcvec, d_as_work)
1130 c---------------------------------------------------------------------
1131 subroutine sddir_verlet1
1132 c Applying velocity Verlet algorithm - step 1 to velocities
1133 implicit real*8 (a-h,o-z)
1134 include 'DIMENSIONS'
1135 include 'COMMON.CONTROL'
1136 include 'COMMON.VAR'
1139 include 'COMMON.LANGEVIN'
1141 include 'COMMON.LANGEVIN.lang0'
1143 include 'COMMON.CHAIN'
1144 include 'COMMON.DERIV'
1145 include 'COMMON.GEO'
1146 include 'COMMON.LOCAL'
1147 include 'COMMON.INTERACT'
1148 include 'COMMON.IOUNITS'
1149 include 'COMMON.NAMES'
1150 c Revised 3/31/05 AL: correlation between random contributions to
1151 c position and velocity increments included.
1152 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1153 double precision adt,adt2
1155 c Add the contribution from BOTH friction and stochastic force to the
1156 c coordinates, but ONLY the contribution from the friction forces to velocities
1159 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1160 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1161 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1162 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1163 d_t(j,0)=d_t_old(j,0)+adt
1168 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1169 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1170 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1171 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1172 d_t(j,i)=d_t_old(j,i)+adt
1177 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1180 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1181 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1182 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1183 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1184 d_t(j,inres)=d_t_old(j,inres)+adt
1191 c---------------------------------------------------------------------
1192 subroutine sddir_verlet2
1193 c Calculating the adjusted velocities for accelerations
1194 implicit real*8 (a-h,o-z)
1195 include 'DIMENSIONS'
1196 include 'COMMON.CONTROL'
1197 include 'COMMON.VAR'
1200 include 'COMMON.LANGEVIN'
1202 include 'COMMON.LANGEVIN.lang0'
1204 include 'COMMON.CHAIN'
1205 include 'COMMON.DERIV'
1206 include 'COMMON.GEO'
1207 include 'COMMON.LOCAL'
1208 include 'COMMON.INTERACT'
1209 include 'COMMON.IOUNITS'
1210 include 'COMMON.NAMES'
1211 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1212 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1213 c Revised 3/31/05 AL: correlation between random contributions to
1214 c position and velocity increments included.
1215 c The correlation coefficients are calculated at low-friction limit.
1216 c Also, friction forces are now not calculated with new velocities.
1218 c call friction_force
1219 call stochastic_force(stochforcvec)
1221 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1222 c forces (d_as_work)
1224 call ginv_mult(stochforcvec, d_as_work1)
1230 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1231 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1236 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1237 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1242 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1245 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1246 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1247 & +cos60*d_as_work1(ind+j))*d_time
1254 c---------------------------------------------------------------------
1255 subroutine max_accel
1257 c Find the maximum difference in the accelerations of the the sites
1258 c at the beginning and the end of the time step.
1260 implicit real*8 (a-h,o-z)
1261 include 'DIMENSIONS'
1262 include 'COMMON.CONTROL'
1263 include 'COMMON.VAR'
1265 include 'COMMON.CHAIN'
1266 include 'COMMON.DERIV'
1267 include 'COMMON.GEO'
1268 include 'COMMON.LOCAL'
1269 include 'COMMON.INTERACT'
1270 include 'COMMON.IOUNITS'
1271 double precision aux(3),accel(3),accel_old(3),dacc
1273 c aux(j)=d_a(j,0)-d_a_old(j,0)
1274 accel_old(j)=d_a_old(j,0)
1281 c 7/3/08 changed to asymmetric difference
1283 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1284 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1285 accel(j)=accel(j)+0.5d0*d_a(j,i)
1286 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1287 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1288 dacc=dabs(accel(j)-accel_old(j))
1289 c write (iout,*) i,dacc
1290 if (dacc.gt.amax) amax=dacc
1298 accel_old(j)=d_a_old(j,0)
1303 accel_old(j)=accel_old(j)+d_a_old(j,1)
1304 accel(j)=accel(j)+d_a(j,1)
1308 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1310 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1311 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1312 accel(j)=accel(j)+d_a(j,i+nres)
1316 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1317 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1318 dacc=dabs(accel(j)-accel_old(j))
1319 c write (iout,*) "side-chain",i,dacc
1320 if (dacc.gt.amax) amax=dacc
1324 accel_old(j)=accel_old(j)+d_a_old(j,i)
1325 accel(j)=accel(j)+d_a(j,i)
1326 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1331 c---------------------------------------------------------------------
1332 subroutine predict_edrift(epdrift)
1334 c Predict the drift of the potential energy
1336 implicit real*8 (a-h,o-z)
1337 include 'DIMENSIONS'
1338 include 'COMMON.CONTROL'
1339 include 'COMMON.VAR'
1341 include 'COMMON.CHAIN'
1342 include 'COMMON.DERIV'
1343 include 'COMMON.GEO'
1344 include 'COMMON.LOCAL'
1345 include 'COMMON.INTERACT'
1346 include 'COMMON.IOUNITS'
1347 include 'COMMON.MUCA'
1348 double precision epdrift,epdriftij
1349 c Drift of the potential energy
1355 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1356 if (lmuca) epdriftij=epdriftij*factor
1357 c write (iout,*) "back",i,j,epdriftij
1358 if (epdriftij.gt.epdrift) epdrift=epdriftij
1362 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1365 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1366 if (lmuca) epdriftij=epdriftij*factor
1367 c write (iout,*) "side",i,j,epdriftij
1368 if (epdriftij.gt.epdrift) epdrift=epdriftij
1372 epdrift=0.5d0*epdrift*d_time*d_time
1373 c write (iout,*) "epdrift",epdrift
1376 c-----------------------------------------------------------------------
1377 subroutine verlet_bath
1379 c Coupling to the thermostat by using the Berendsen algorithm
1381 implicit real*8 (a-h,o-z)
1382 include 'DIMENSIONS'
1383 include 'COMMON.CONTROL'
1384 include 'COMMON.VAR'
1386 include 'COMMON.CHAIN'
1387 include 'COMMON.DERIV'
1388 include 'COMMON.GEO'
1389 include 'COMMON.LOCAL'
1390 include 'COMMON.INTERACT'
1391 include 'COMMON.IOUNITS'
1392 include 'COMMON.NAMES'
1393 double precision T_half,fact
1395 T_half=2.0d0/(dimen3*Rb)*EK
1396 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1397 c write(iout,*) "T_half", T_half
1398 c write(iout,*) "EK", EK
1399 c write(iout,*) "fact", fact
1401 d_t(j,0)=fact*d_t(j,0)
1405 d_t(j,i)=fact*d_t(j,i)
1409 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1412 d_t(j,inres)=fact*d_t(j,inres)
1418 c---------------------------------------------------------
1420 c Set up the initial conditions of a MD simulation
1421 implicit real*8 (a-h,o-z)
1422 include 'DIMENSIONS'
1426 integer IERROR,ERRCODE
1428 include 'COMMON.SETUP'
1429 include 'COMMON.CONTROL'
1430 include 'COMMON.VAR'
1433 include 'COMMON.LANGEVIN'
1435 include 'COMMON.LANGEVIN.lang0'
1437 include 'COMMON.CHAIN'
1438 include 'COMMON.DERIV'
1439 include 'COMMON.GEO'
1440 include 'COMMON.LOCAL'
1441 include 'COMMON.INTERACT'
1442 include 'COMMON.IOUNITS'
1443 include 'COMMON.NAMES'
1444 include 'COMMON.REMD'
1445 real*8 energia_long(0:n_ene),
1446 & energia_short(0:n_ene),vcm(3),incr(3)
1447 double precision cm(3),L(3),xv,sigv,lowb,highb
1448 double precision varia(maxvar)
1456 c write(iout,*) "d_time", d_time
1457 c Compute the standard deviations of stochastic forces for Langevin dynamics
1458 c if the friction coefficients do not depend on surface area
1459 if (lang.gt.0 .and. .not.surfarea) then
1461 stdforcp(i)=stdfp*dsqrt(gamp)
1464 stdforcsc(i)=stdfsc(iabs(itype(i)))
1465 & *dsqrt(gamsc(iabs(itype(i))))
1468 c Open the pdb file for snapshotshots
1471 if (ilen(tmpdir).gt.0)
1472 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1473 & liczba(:ilen(liczba))//".pdb")
1475 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1479 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1480 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1481 & liczba(:ilen(liczba))//".x")
1482 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1485 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1486 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1487 & liczba(:ilen(liczba))//".cx")
1488 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1494 if (ilen(tmpdir).gt.0)
1495 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1496 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1498 if (ilen(tmpdir).gt.0)
1499 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1500 cartname=prefix(:ilen(prefix))//"_MD.cx"
1504 write (qstr,'(256(1h ))')
1507 iq = qinfrag(i,iset)*10
1508 iw = wfrag(i,iset)/100
1510 if(me.eq.king.or..not.out1file)
1511 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1512 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1517 iq = qinpair(i,iset)*10
1518 iw = wpair(i,iset)/100
1520 if(me.eq.king.or..not.out1file)
1521 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1522 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1526 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1528 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1530 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1532 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1536 if (restart1file) then
1538 & inquire(file=mremd_rst_name,exist=file_exist)
1540 write (*,*) me," Before broadcast: file_exist",file_exist
1541 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1543 write (*,*) me," After broadcast: file_exist",file_exist
1544 c inquire(file=mremd_rst_name,exist=file_exist)
1546 if(me.eq.king.or..not.out1file)
1547 & write(iout,*) "Initial state read by master and distributed"
1549 if (ilen(tmpdir).gt.0)
1550 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1551 & //liczba(:ilen(liczba))//'.rst')
1552 inquire(file=rest2name,exist=file_exist)
1555 if(.not.restart1file) then
1556 if(me.eq.king.or..not.out1file)
1557 & write(iout,*) "Initial state will be read from file ",
1558 & rest2name(:ilen(rest2name))
1561 call rescale_weights(t_bath)
1563 if(me.eq.king.or..not.out1file)then
1564 if (restart1file) then
1565 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1568 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1571 write(iout,*) "Initial velocities randomly generated"
1578 c Generate initial velocities
1579 if(me.eq.king.or..not.out1file)
1580 & write(iout,*) "Initial velocities randomly generated"
1583 CtotTafm is the variable for AFM time which eclipsed during
1586 c rest2name = prefix(:ilen(prefix))//'.rst'
1587 if(me.eq.king.or..not.out1file)then
1588 write (iout,*) "Initial velocities"
1590 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1591 & (d_t(j,i+nres),j=1,3)
1593 c Zeroing the total angular momentum of the system
1594 write(iout,*) "Calling the zero-angular
1595 & momentum subroutine"
1598 c Getting the potential energy and forces and velocities and accelerations
1600 c write (iout,*) "velocity of the center of the mass:"
1601 c write (iout,*) (vcm(j),j=1,3)
1603 d_t(j,0)=d_t(j,0)-vcm(j)
1605 c Removing the velocity of the center of mass
1607 if(me.eq.king.or..not.out1file)then
1608 write (iout,*) "vcm right after adjustment:"
1609 write (iout,*) (vcm(j),j=1,3)
1613 if(iranconf.ne.0) then
1615 print *, 'Calling OVERLAP_SC'
1616 call overlap_sc(fail)
1620 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1621 print *,'SC_move',nft_sc,etot
1622 if(me.eq.king.or..not.out1file)
1623 & write(iout,*) 'SC_move',nft_sc,etot
1627 print *, 'Calling MINIM_DC'
1628 call minim_dc(etot,iretcode,nfun)
1630 call geom_to_var(nvar,varia)
1631 print *,'Calling MINIMIZE.'
1632 call minimize(etot,varia,iretcode,nfun)
1633 call var_to_geom(nvar,varia)
1635 if(me.eq.king.or..not.out1file)
1636 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1639 call chainbuild_cart
1644 kinetic_T=2.0d0/(dimen3*Rb)*EK
1645 if(me.eq.king.or..not.out1file)then
1655 call etotal(potEcomp)
1656 if (large) call enerprint(potEcomp)
1659 t_etotal=t_etotal+MPI_Wtime()-tt0
1661 t_etotal=t_etotal+tcpu()-tt0
1668 if (amax*d_time .gt. dvmax) then
1669 d_time=d_time*dvmax/amax
1670 if(me.eq.king.or..not.out1file) write (iout,*)
1671 & "Time step reduced to",d_time,
1672 & " because of too large initial acceleration."
1674 if(me.eq.king.or..not.out1file)then
1675 write(iout,*) "Potential energy and its components"
1676 call enerprint(potEcomp)
1677 c write(iout,*) (potEcomp(i),i=0,n_ene)
1679 potE=potEcomp(0)-potEcomp(20)
1682 if (ntwe.ne.0) call statout(itime)
1683 if(me.eq.king.or..not.out1file)
1684 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1685 & " Kinetic energy",EK," Potential energy",potE,
1686 & " Total energy",totE," Maximum acceleration ",
1689 write (iout,*) "Initial coordinates"
1691 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1692 & (c(j,i+nres),j=1,3)
1694 write (iout,*) "Initial dC"
1696 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1697 & (dc(j,i+nres),j=1,3)
1699 write (iout,*) "Initial velocities"
1700 write (iout,"(13x,' backbone ',23x,' side chain')")
1702 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1703 & (d_t(j,i+nres),j=1,3)
1705 write (iout,*) "Initial accelerations"
1707 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1708 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1709 & (d_a(j,i+nres),j=1,3)
1715 d_t_old(j,i)=d_t(j,i)
1716 d_a_old(j,i)=d_a(j,i)
1718 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1727 call etotal_short(energia_short)
1728 if (large) call enerprint(potEcomp)
1731 t_eshort=t_eshort+MPI_Wtime()-tt0
1733 t_eshort=t_eshort+tcpu()-tt0
1738 if(.not.out1file .and. large) then
1739 write (iout,*) "energia_long",energia_long(0),
1740 & " energia_short",energia_short(0),
1741 & " total",energia_long(0)+energia_short(0)
1742 write (iout,*) "Initial fast-force accelerations"
1744 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1745 & (d_a(j,i+nres),j=1,3)
1748 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1751 d_a_short(j,i)=d_a(j,i)
1760 call etotal_long(energia_long)
1761 if (large) call enerprint(potEcomp)
1764 t_elong=t_elong+MPI_Wtime()-tt0
1766 t_elong=t_elong+tcpu()-tt0
1771 if(.not.out1file .and. large) then
1772 write (iout,*) "energia_long",energia_long(0)
1773 write (iout,*) "Initial slow-force accelerations"
1775 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1776 & (d_a(j,i+nres),j=1,3)
1780 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1782 t_enegrad=t_enegrad+tcpu()-tt0
1787 c-----------------------------------------------------------
1788 subroutine random_vel
1789 implicit real*8 (a-h,o-z)
1790 include 'DIMENSIONS'
1791 include 'COMMON.CONTROL'
1792 include 'COMMON.VAR'
1795 include 'COMMON.LANGEVIN'
1797 include 'COMMON.LANGEVIN.lang0'
1799 include 'COMMON.CHAIN'
1800 include 'COMMON.DERIV'
1801 include 'COMMON.GEO'
1802 include 'COMMON.LOCAL'
1803 include 'COMMON.INTERACT'
1804 include 'COMMON.IOUNITS'
1805 include 'COMMON.NAMES'
1806 include 'COMMON.TIME1'
1807 double precision xv,sigv,lowb,highb,vec_afm(3)
1808 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1809 c First generate velocities in the eigenspace of the G matrix
1810 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1817 sigv=dsqrt((Rb*t_bath)/geigen(i))
1820 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1822 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1823 c & " d_t_work_new",d_t_work_new(ii)
1826 C if (SELFGUIDE.gt.0) then
1829 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
1830 C distance=distance+vec_afm(j)**2
1832 C distance=dsqrt(distance)
1834 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
1835 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
1836 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
1837 C & d_t_work_new(j+(afmend-1)*3)
1848 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1851 c write (iout,*) "Ek from eigenvectors",Ek1
1853 c Transform velocities to UNRES coordinate space
1859 d_t_work(ind)=d_t_work(ind)
1860 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1862 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1866 c Transfer to the d_t vector
1868 d_t(j,0)=d_t_work(j)
1874 d_t(j,i)=d_t_work(ind)
1878 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1881 d_t(j,i+nres)=d_t_work(ind)
1886 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1887 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1892 c-----------------------------------------------------------
1893 subroutine sd_verlet_p_setup
1894 c Sets up the parameters of stochastic Verlet algorithm
1895 implicit real*8 (a-h,o-z)
1896 include 'DIMENSIONS'
1900 include 'COMMON.CONTROL'
1901 include 'COMMON.VAR'
1904 include 'COMMON.LANGEVIN'
1906 include 'COMMON.LANGEVIN.lang0'
1908 include 'COMMON.CHAIN'
1909 include 'COMMON.DERIV'
1910 include 'COMMON.GEO'
1911 include 'COMMON.LOCAL'
1912 include 'COMMON.INTERACT'
1913 include 'COMMON.IOUNITS'
1914 include 'COMMON.NAMES'
1915 include 'COMMON.TIME1'
1916 double precision emgdt(MAXRES6),
1917 & pterm,vterm,rho,rhoc,vsig,
1918 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1919 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1920 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1921 logical lprn /.false./
1922 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1923 double precision ktm
1930 c AL 8/17/04 Code adapted from tinker
1932 c Get the frictional and random terms for stochastic dynamics in the
1933 c eigenspace of mass-scaled UNRES friction matrix
1936 gdt = fricgam(i) * d_time
1938 c Stochastic dynamics reduces to simple MD for zero friction
1940 if (gdt .le. zero) then
1941 pfric_vec(i) = 1.0d0
1942 vfric_vec(i) = d_time
1943 afric_vec(i) = 0.5d0 * d_time * d_time
1944 prand_vec(i) = 0.0d0
1945 vrand_vec1(i) = 0.0d0
1946 vrand_vec2(i) = 0.0d0
1948 c Analytical expressions when friction coefficient is large
1951 if (gdt .ge. gdt_radius) then
1954 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1955 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1956 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1957 vterm = 1.0d0 - egdt**2
1958 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1960 c Use series expansions when friction coefficient is small
1971 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1972 & - gdt5/120.0d0 + gdt6/720.0d0
1973 & - gdt7/5040.0d0 + gdt8/40320.0d0
1974 & - gdt9/362880.0d0) / fricgam(i)**2
1975 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1976 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1977 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1978 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1979 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1980 & + 127.0d0*gdt9/90720.0d0
1981 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1982 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1983 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1984 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1985 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1986 & - 17.0d0*gdt2/1280.0d0
1987 & + 17.0d0*gdt3/6144.0d0
1988 & + 40967.0d0*gdt4/34406400.0d0
1989 & - 57203.0d0*gdt5/275251200.0d0
1990 & - 1429487.0d0*gdt6/13212057600.0d0)
1993 c Compute the scaling factors of random terms for the nonzero friction case
1995 ktm = 0.5d0*d_time/fricgam(i)
1996 psig = dsqrt(ktm*pterm) / fricgam(i)
1997 vsig = dsqrt(ktm*vterm)
1998 rhoc = dsqrt(1.0d0 - rho*rho)
2000 vrand_vec1(i) = vsig * rho
2001 vrand_vec2(i) = vsig * rhoc
2006 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2009 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2010 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2014 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2017 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2018 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2019 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2020 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2021 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2022 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2025 t_sdsetup=t_sdsetup+MPI_Wtime()
2027 t_sdsetup=t_sdsetup+tcpu()-tt0
2031 c-------------------------------------------------------------
2032 subroutine eigtransf1(n,ndim,ab,d,c)
2035 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2041 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2047 c-------------------------------------------------------------
2048 subroutine eigtransf(n,ndim,a,b,d,c)
2051 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2057 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2063 c-------------------------------------------------------------
2064 subroutine sd_verlet1
2065 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2066 implicit real*8 (a-h,o-z)
2067 include 'DIMENSIONS'
2068 include 'COMMON.CONTROL'
2069 include 'COMMON.VAR'
2072 include 'COMMON.LANGEVIN'
2074 include 'COMMON.LANGEVIN.lang0'
2076 include 'COMMON.CHAIN'
2077 include 'COMMON.DERIV'
2078 include 'COMMON.GEO'
2079 include 'COMMON.LOCAL'
2080 include 'COMMON.INTERACT'
2081 include 'COMMON.IOUNITS'
2082 include 'COMMON.NAMES'
2083 double precision stochforcvec(MAXRES6)
2084 common /stochcalc/ stochforcvec
2085 logical lprn /.false./
2087 c write (iout,*) "dc_old"
2089 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2090 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2093 dc_work(j)=dc_old(j,0)
2094 d_t_work(j)=d_t_old(j,0)
2095 d_a_work(j)=d_a_old(j,0)
2100 dc_work(ind+j)=dc_old(j,i)
2101 d_t_work(ind+j)=d_t_old(j,i)
2102 d_a_work(ind+j)=d_a_old(j,i)
2107 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2109 dc_work(ind+j)=dc_old(j,i+nres)
2110 d_t_work(ind+j)=d_t_old(j,i+nres)
2111 d_a_work(ind+j)=d_a_old(j,i+nres)
2119 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2123 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2124 & vfric_mat(i,j),afric_mat(i,j),
2125 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2133 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2134 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2135 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2136 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2138 d_t_work_new(i)=ddt1+0.5d0*ddt2
2139 d_t_work(i)=ddt1+ddt2
2144 d_t(j,0)=d_t_work(j)
2149 dc(j,i)=dc_work(ind+j)
2150 d_t(j,i)=d_t_work(ind+j)
2155 if (itype(i).ne.10) then
2158 dc(j,inres)=dc_work(ind+j)
2159 d_t(j,inres)=d_t_work(ind+j)
2166 c--------------------------------------------------------------------------
2167 subroutine sd_verlet2
2168 c Calculating the adjusted velocities for accelerations
2169 implicit real*8 (a-h,o-z)
2170 include 'DIMENSIONS'
2171 include 'COMMON.CONTROL'
2172 include 'COMMON.VAR'
2175 include 'COMMON.LANGEVIN'
2177 include 'COMMON.LANGEVIN.lang0'
2179 include 'COMMON.CHAIN'
2180 include 'COMMON.DERIV'
2181 include 'COMMON.GEO'
2182 include 'COMMON.LOCAL'
2183 include 'COMMON.INTERACT'
2184 include 'COMMON.IOUNITS'
2185 include 'COMMON.NAMES'
2186 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2187 common /stochcalc/ stochforcvec
2189 c Compute the stochastic forces which contribute to velocity change
2191 call stochastic_force(stochforcvecV)
2198 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2199 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2200 & vrand_mat2(i,j)*stochforcvecV(j)
2202 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2206 d_t(j,0)=d_t_work(j)
2211 d_t(j,i)=d_t_work(ind+j)
2216 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2219 d_t(j,inres)=d_t_work(ind+j)
2226 c-----------------------------------------------------------
2227 subroutine sd_verlet_ciccotti_setup
2228 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2230 implicit real*8 (a-h,o-z)
2231 include 'DIMENSIONS'
2235 include 'COMMON.CONTROL'
2236 include 'COMMON.VAR'
2239 include 'COMMON.LANGEVIN'
2241 include 'COMMON.LANGEVIN.lang0'
2243 include 'COMMON.CHAIN'
2244 include 'COMMON.DERIV'
2245 include 'COMMON.GEO'
2246 include 'COMMON.LOCAL'
2247 include 'COMMON.INTERACT'
2248 include 'COMMON.IOUNITS'
2249 include 'COMMON.NAMES'
2250 include 'COMMON.TIME1'
2251 double precision emgdt(MAXRES6),
2252 & pterm,vterm,rho,rhoc,vsig,
2253 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2254 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2255 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2256 logical lprn /.false./
2257 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2258 double precision ktm
2265 c AL 8/17/04 Code adapted from tinker
2267 c Get the frictional and random terms for stochastic dynamics in the
2268 c eigenspace of mass-scaled UNRES friction matrix
2271 write (iout,*) "i",i," fricgam",fricgam(i)
2272 gdt = fricgam(i) * d_time
2274 c Stochastic dynamics reduces to simple MD for zero friction
2276 if (gdt .le. zero) then
2277 pfric_vec(i) = 1.0d0
2278 vfric_vec(i) = d_time
2279 afric_vec(i) = 0.5d0*d_time*d_time
2280 prand_vec(i) = afric_vec(i)
2281 vrand_vec2(i) = vfric_vec(i)
2283 c Analytical expressions when friction coefficient is large
2288 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2289 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2290 prand_vec(i) = afric_vec(i)
2291 vrand_vec2(i) = vfric_vec(i)
2293 c Compute the scaling factors of random terms for the nonzero friction case
2295 c ktm = 0.5d0*d_time/fricgam(i)
2296 c psig = dsqrt(ktm*pterm) / fricgam(i)
2297 c vsig = dsqrt(ktm*vterm)
2298 c prand_vec(i) = psig*afric_vec(i)
2299 c vrand_vec2(i) = vsig*vfric_vec(i)
2304 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2307 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2308 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2312 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2314 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2315 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2316 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2317 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2318 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2320 t_sdsetup=t_sdsetup+MPI_Wtime()
2322 t_sdsetup=t_sdsetup+tcpu()-tt0
2326 c-------------------------------------------------------------
2327 subroutine sd_verlet1_ciccotti
2328 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2329 implicit real*8 (a-h,o-z)
2330 include 'DIMENSIONS'
2334 include 'COMMON.CONTROL'
2335 include 'COMMON.VAR'
2338 include 'COMMON.LANGEVIN'
2340 include 'COMMON.LANGEVIN.lang0'
2342 include 'COMMON.CHAIN'
2343 include 'COMMON.DERIV'
2344 include 'COMMON.GEO'
2345 include 'COMMON.LOCAL'
2346 include 'COMMON.INTERACT'
2347 include 'COMMON.IOUNITS'
2348 include 'COMMON.NAMES'
2349 double precision stochforcvec(MAXRES6)
2350 common /stochcalc/ stochforcvec
2351 logical lprn /.false./
2353 c write (iout,*) "dc_old"
2355 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2356 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2359 dc_work(j)=dc_old(j,0)
2360 d_t_work(j)=d_t_old(j,0)
2361 d_a_work(j)=d_a_old(j,0)
2366 dc_work(ind+j)=dc_old(j,i)
2367 d_t_work(ind+j)=d_t_old(j,i)
2368 d_a_work(ind+j)=d_a_old(j,i)
2373 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2375 dc_work(ind+j)=dc_old(j,i+nres)
2376 d_t_work(ind+j)=d_t_old(j,i+nres)
2377 d_a_work(ind+j)=d_a_old(j,i+nres)
2386 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2390 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2391 & vfric_mat(i,j),afric_mat(i,j),
2392 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2400 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2401 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2402 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2403 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2405 d_t_work_new(i)=ddt1+0.5d0*ddt2
2406 d_t_work(i)=ddt1+ddt2
2411 d_t(j,0)=d_t_work(j)
2416 dc(j,i)=dc_work(ind+j)
2417 d_t(j,i)=d_t_work(ind+j)
2422 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2425 dc(j,inres)=dc_work(ind+j)
2426 d_t(j,inres)=d_t_work(ind+j)
2433 c--------------------------------------------------------------------------
2434 subroutine sd_verlet2_ciccotti
2435 c Calculating the adjusted velocities for accelerations
2436 implicit real*8 (a-h,o-z)
2437 include 'DIMENSIONS'
2438 include 'COMMON.CONTROL'
2439 include 'COMMON.VAR'
2442 include 'COMMON.LANGEVIN'
2444 include 'COMMON.LANGEVIN.lang0'
2446 include 'COMMON.CHAIN'
2447 include 'COMMON.DERIV'
2448 include 'COMMON.GEO'
2449 include 'COMMON.LOCAL'
2450 include 'COMMON.INTERACT'
2451 include 'COMMON.IOUNITS'
2452 include 'COMMON.NAMES'
2453 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2454 common /stochcalc/ stochforcvec
2456 c Compute the stochastic forces which contribute to velocity change
2458 call stochastic_force(stochforcvecV)
2465 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2466 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2467 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2469 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2473 d_t(j,0)=d_t_work(j)
2478 d_t(j,i)=d_t_work(ind+j)
2483 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2486 d_t(j,inres)=d_t_work(ind+j)