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)
1612 if (start_from_model) then
1613 i_model=iran_num(1,constr_homology)
1614 write (iout,*) 'starting from model ',i_model
1617 c(j,i)=chomo(j,i,i_model)
1621 call int_from_cart(.true.,.false.)
1622 call sc_loc_geom(.false.)
1625 dc(j,i)=c(j,i+1)-c(j,i)
1626 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1631 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1632 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1636 c call chainbuild_cart
1637 write (iout,*) "PREMINIM ",preminim
1638 if(iranconf.ne.0 .or. preminim) then
1640 print *, 'Calling OVERLAP_SC'
1641 call overlap_sc(fail)
1645 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1646 print *,'SC_move',nft_sc,etot
1647 if(me.eq.king.or..not.out1file)
1648 & write(iout,*) 'SC_move',nft_sc,etot
1652 print *, 'Calling MINIM_DC'
1653 call minim_dc(etot,iretcode,nfun)
1655 call geom_to_var(nvar,varia)
1656 print *,'Calling MINIMIZE.'
1657 call minimize(etot,varia,iretcode,nfun)
1658 call var_to_geom(nvar,varia)
1660 if(me.eq.king.or..not.out1file) then
1661 write(iout,*) "Minimized energy is",etot
1662 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1663 call etotal(potEcomp)
1664 call enerprint(potEcomp)
1668 call chainbuild_cart
1673 kinetic_T=2.0d0/(dimen3*Rb)*EK
1674 if(me.eq.king.or..not.out1file)then
1684 call etotal(potEcomp)
1685 if (large) call enerprint(potEcomp)
1688 t_etotal=t_etotal+MPI_Wtime()-tt0
1690 t_etotal=t_etotal+tcpu()-tt0
1697 if (amax*d_time .gt. dvmax) then
1698 d_time=d_time*dvmax/amax
1699 if(me.eq.king.or..not.out1file) write (iout,*)
1700 & "Time step reduced to",d_time,
1701 & " because of too large initial acceleration."
1703 C if(me.eq.king.or..not.out1file)then
1704 C write(iout,*) "Potential energy and its components"
1705 C call enerprint(potEcomp)
1706 c write(iout,*) (potEcomp(i),i=0,n_ene)
1708 potE=potEcomp(0)-potEcomp(20)
1711 if (ntwe.ne.0) call statout(itime)
1712 if(me.eq.king.or..not.out1file)
1713 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1714 & " Kinetic energy",EK," Potential energy",potE,
1715 & " Total energy",totE," Maximum acceleration ",
1718 write (iout,*) "Initial coordinates"
1720 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1721 & (c(j,i+nres),j=1,3)
1723 write (iout,*) "Initial dC"
1725 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1726 & (dc(j,i+nres),j=1,3)
1728 write (iout,*) "Initial velocities"
1729 write (iout,"(13x,' backbone ',23x,' side chain')")
1731 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1732 & (d_t(j,i+nres),j=1,3)
1734 write (iout,*) "Initial accelerations"
1736 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1737 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1738 & (d_a(j,i+nres),j=1,3)
1744 d_t_old(j,i)=d_t(j,i)
1745 d_a_old(j,i)=d_a(j,i)
1747 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1756 call etotal_short(energia_short)
1757 if (large) call enerprint(potEcomp)
1760 t_eshort=t_eshort+MPI_Wtime()-tt0
1762 t_eshort=t_eshort+tcpu()-tt0
1767 if(.not.out1file .and. large) then
1768 write (iout,*) "energia_long",energia_long(0),
1769 & " energia_short",energia_short(0),
1770 & " total",energia_long(0)+energia_short(0)
1771 write (iout,*) "Initial fast-force accelerations"
1773 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1774 & (d_a(j,i+nres),j=1,3)
1777 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1780 d_a_short(j,i)=d_a(j,i)
1789 call etotal_long(energia_long)
1790 if (large) call enerprint(potEcomp)
1793 t_elong=t_elong+MPI_Wtime()-tt0
1795 t_elong=t_elong+tcpu()-tt0
1800 if(.not.out1file .and. large) then
1801 write (iout,*) "energia_long",energia_long(0)
1802 write (iout,*) "Initial slow-force accelerations"
1804 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1805 & (d_a(j,i+nres),j=1,3)
1809 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1811 t_enegrad=t_enegrad+tcpu()-tt0
1816 c-----------------------------------------------------------
1817 subroutine random_vel
1818 implicit real*8 (a-h,o-z)
1819 include 'DIMENSIONS'
1820 include 'COMMON.CONTROL'
1821 include 'COMMON.VAR'
1824 include 'COMMON.LANGEVIN'
1826 include 'COMMON.LANGEVIN.lang0'
1828 include 'COMMON.CHAIN'
1829 include 'COMMON.DERIV'
1830 include 'COMMON.GEO'
1831 include 'COMMON.LOCAL'
1832 include 'COMMON.INTERACT'
1833 include 'COMMON.IOUNITS'
1834 include 'COMMON.NAMES'
1835 include 'COMMON.TIME1'
1836 double precision xv,sigv,lowb,highb,vec_afm(3)
1837 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1838 c First generate velocities in the eigenspace of the G matrix
1839 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1846 sigv=dsqrt((Rb*t_bath)/geigen(i))
1849 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1851 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1852 c & " d_t_work_new",d_t_work_new(ii)
1855 C if (SELFGUIDE.gt.0) then
1858 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
1859 C distance=distance+vec_afm(j)**2
1861 C distance=dsqrt(distance)
1863 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
1864 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
1865 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
1866 C & d_t_work_new(j+(afmend-1)*3)
1877 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1880 c write (iout,*) "Ek from eigenvectors",Ek1
1882 c Transform velocities to UNRES coordinate space
1888 d_t_work(ind)=d_t_work(ind)
1889 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1891 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1895 c Transfer to the d_t vector
1897 d_t(j,0)=d_t_work(j)
1903 d_t(j,i)=d_t_work(ind)
1907 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1910 d_t(j,i+nres)=d_t_work(ind)
1915 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1916 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1921 c-----------------------------------------------------------
1922 subroutine sd_verlet_p_setup
1923 c Sets up the parameters of stochastic Verlet algorithm
1924 implicit real*8 (a-h,o-z)
1925 include 'DIMENSIONS'
1929 include 'COMMON.CONTROL'
1930 include 'COMMON.VAR'
1933 include 'COMMON.LANGEVIN'
1935 include 'COMMON.LANGEVIN.lang0'
1937 include 'COMMON.CHAIN'
1938 include 'COMMON.DERIV'
1939 include 'COMMON.GEO'
1940 include 'COMMON.LOCAL'
1941 include 'COMMON.INTERACT'
1942 include 'COMMON.IOUNITS'
1943 include 'COMMON.NAMES'
1944 include 'COMMON.TIME1'
1945 double precision emgdt(MAXRES6),
1946 & pterm,vterm,rho,rhoc,vsig,
1947 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1948 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1949 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1950 logical lprn /.false./
1951 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1952 double precision ktm
1959 c AL 8/17/04 Code adapted from tinker
1961 c Get the frictional and random terms for stochastic dynamics in the
1962 c eigenspace of mass-scaled UNRES friction matrix
1965 gdt = fricgam(i) * d_time
1967 c Stochastic dynamics reduces to simple MD for zero friction
1969 if (gdt .le. zero) then
1970 pfric_vec(i) = 1.0d0
1971 vfric_vec(i) = d_time
1972 afric_vec(i) = 0.5d0 * d_time * d_time
1973 prand_vec(i) = 0.0d0
1974 vrand_vec1(i) = 0.0d0
1975 vrand_vec2(i) = 0.0d0
1977 c Analytical expressions when friction coefficient is large
1980 if (gdt .ge. gdt_radius) then
1983 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1984 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1985 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1986 vterm = 1.0d0 - egdt**2
1987 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1989 c Use series expansions when friction coefficient is small
2000 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2001 & - gdt5/120.0d0 + gdt6/720.0d0
2002 & - gdt7/5040.0d0 + gdt8/40320.0d0
2003 & - gdt9/362880.0d0) / fricgam(i)**2
2004 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2005 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2006 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2007 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2008 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2009 & + 127.0d0*gdt9/90720.0d0
2010 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2011 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2012 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2013 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2014 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2015 & - 17.0d0*gdt2/1280.0d0
2016 & + 17.0d0*gdt3/6144.0d0
2017 & + 40967.0d0*gdt4/34406400.0d0
2018 & - 57203.0d0*gdt5/275251200.0d0
2019 & - 1429487.0d0*gdt6/13212057600.0d0)
2022 c Compute the scaling factors of random terms for the nonzero friction case
2024 ktm = 0.5d0*d_time/fricgam(i)
2025 psig = dsqrt(ktm*pterm) / fricgam(i)
2026 vsig = dsqrt(ktm*vterm)
2027 rhoc = dsqrt(1.0d0 - rho*rho)
2029 vrand_vec1(i) = vsig * rho
2030 vrand_vec2(i) = vsig * rhoc
2035 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2038 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2039 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2043 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2046 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2047 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2048 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2049 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2050 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2051 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2054 t_sdsetup=t_sdsetup+MPI_Wtime()
2056 t_sdsetup=t_sdsetup+tcpu()-tt0
2060 c-------------------------------------------------------------
2061 subroutine eigtransf1(n,ndim,ab,d,c)
2064 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2070 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2076 c-------------------------------------------------------------
2077 subroutine eigtransf(n,ndim,a,b,d,c)
2080 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2086 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2092 c-------------------------------------------------------------
2093 subroutine sd_verlet1
2094 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2095 implicit real*8 (a-h,o-z)
2096 include 'DIMENSIONS'
2097 include 'COMMON.CONTROL'
2098 include 'COMMON.VAR'
2101 include 'COMMON.LANGEVIN'
2103 include 'COMMON.LANGEVIN.lang0'
2105 include 'COMMON.CHAIN'
2106 include 'COMMON.DERIV'
2107 include 'COMMON.GEO'
2108 include 'COMMON.LOCAL'
2109 include 'COMMON.INTERACT'
2110 include 'COMMON.IOUNITS'
2111 include 'COMMON.NAMES'
2112 double precision stochforcvec(MAXRES6)
2113 common /stochcalc/ stochforcvec
2114 logical lprn /.false./
2116 c write (iout,*) "dc_old"
2118 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2119 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2122 dc_work(j)=dc_old(j,0)
2123 d_t_work(j)=d_t_old(j,0)
2124 d_a_work(j)=d_a_old(j,0)
2129 dc_work(ind+j)=dc_old(j,i)
2130 d_t_work(ind+j)=d_t_old(j,i)
2131 d_a_work(ind+j)=d_a_old(j,i)
2136 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2138 dc_work(ind+j)=dc_old(j,i+nres)
2139 d_t_work(ind+j)=d_t_old(j,i+nres)
2140 d_a_work(ind+j)=d_a_old(j,i+nres)
2148 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2152 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2153 & vfric_mat(i,j),afric_mat(i,j),
2154 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2162 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2163 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2164 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2165 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2167 d_t_work_new(i)=ddt1+0.5d0*ddt2
2168 d_t_work(i)=ddt1+ddt2
2173 d_t(j,0)=d_t_work(j)
2178 dc(j,i)=dc_work(ind+j)
2179 d_t(j,i)=d_t_work(ind+j)
2184 if (itype(i).ne.10) then
2187 dc(j,inres)=dc_work(ind+j)
2188 d_t(j,inres)=d_t_work(ind+j)
2195 c--------------------------------------------------------------------------
2196 subroutine sd_verlet2
2197 c Calculating the adjusted velocities for accelerations
2198 implicit real*8 (a-h,o-z)
2199 include 'DIMENSIONS'
2200 include 'COMMON.CONTROL'
2201 include 'COMMON.VAR'
2204 include 'COMMON.LANGEVIN'
2206 include 'COMMON.LANGEVIN.lang0'
2208 include 'COMMON.CHAIN'
2209 include 'COMMON.DERIV'
2210 include 'COMMON.GEO'
2211 include 'COMMON.LOCAL'
2212 include 'COMMON.INTERACT'
2213 include 'COMMON.IOUNITS'
2214 include 'COMMON.NAMES'
2215 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2216 common /stochcalc/ stochforcvec
2218 c Compute the stochastic forces which contribute to velocity change
2220 call stochastic_force(stochforcvecV)
2227 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2228 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2229 & vrand_mat2(i,j)*stochforcvecV(j)
2231 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2235 d_t(j,0)=d_t_work(j)
2240 d_t(j,i)=d_t_work(ind+j)
2245 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2248 d_t(j,inres)=d_t_work(ind+j)
2255 c-----------------------------------------------------------
2256 subroutine sd_verlet_ciccotti_setup
2257 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2259 implicit real*8 (a-h,o-z)
2260 include 'DIMENSIONS'
2264 include 'COMMON.CONTROL'
2265 include 'COMMON.VAR'
2268 include 'COMMON.LANGEVIN'
2270 include 'COMMON.LANGEVIN.lang0'
2272 include 'COMMON.CHAIN'
2273 include 'COMMON.DERIV'
2274 include 'COMMON.GEO'
2275 include 'COMMON.LOCAL'
2276 include 'COMMON.INTERACT'
2277 include 'COMMON.IOUNITS'
2278 include 'COMMON.NAMES'
2279 include 'COMMON.TIME1'
2280 double precision emgdt(MAXRES6),
2281 & pterm,vterm,rho,rhoc,vsig,
2282 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2283 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2284 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2285 logical lprn /.false./
2286 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2287 double precision ktm
2294 c AL 8/17/04 Code adapted from tinker
2296 c Get the frictional and random terms for stochastic dynamics in the
2297 c eigenspace of mass-scaled UNRES friction matrix
2300 write (iout,*) "i",i," fricgam",fricgam(i)
2301 gdt = fricgam(i) * d_time
2303 c Stochastic dynamics reduces to simple MD for zero friction
2305 if (gdt .le. zero) then
2306 pfric_vec(i) = 1.0d0
2307 vfric_vec(i) = d_time
2308 afric_vec(i) = 0.5d0*d_time*d_time
2309 prand_vec(i) = afric_vec(i)
2310 vrand_vec2(i) = vfric_vec(i)
2312 c Analytical expressions when friction coefficient is large
2317 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2318 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2319 prand_vec(i) = afric_vec(i)
2320 vrand_vec2(i) = vfric_vec(i)
2322 c Compute the scaling factors of random terms for the nonzero friction case
2324 c ktm = 0.5d0*d_time/fricgam(i)
2325 c psig = dsqrt(ktm*pterm) / fricgam(i)
2326 c vsig = dsqrt(ktm*vterm)
2327 c prand_vec(i) = psig*afric_vec(i)
2328 c vrand_vec2(i) = vsig*vfric_vec(i)
2333 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2336 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2337 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2341 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2343 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2344 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2345 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2346 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2347 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2349 t_sdsetup=t_sdsetup+MPI_Wtime()
2351 t_sdsetup=t_sdsetup+tcpu()-tt0
2355 c-------------------------------------------------------------
2356 subroutine sd_verlet1_ciccotti
2357 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2358 implicit real*8 (a-h,o-z)
2359 include 'DIMENSIONS'
2363 include 'COMMON.CONTROL'
2364 include 'COMMON.VAR'
2367 include 'COMMON.LANGEVIN'
2369 include 'COMMON.LANGEVIN.lang0'
2371 include 'COMMON.CHAIN'
2372 include 'COMMON.DERIV'
2373 include 'COMMON.GEO'
2374 include 'COMMON.LOCAL'
2375 include 'COMMON.INTERACT'
2376 include 'COMMON.IOUNITS'
2377 include 'COMMON.NAMES'
2378 double precision stochforcvec(MAXRES6)
2379 common /stochcalc/ stochforcvec
2380 logical lprn /.false./
2382 c write (iout,*) "dc_old"
2384 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2385 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2388 dc_work(j)=dc_old(j,0)
2389 d_t_work(j)=d_t_old(j,0)
2390 d_a_work(j)=d_a_old(j,0)
2395 dc_work(ind+j)=dc_old(j,i)
2396 d_t_work(ind+j)=d_t_old(j,i)
2397 d_a_work(ind+j)=d_a_old(j,i)
2402 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2404 dc_work(ind+j)=dc_old(j,i+nres)
2405 d_t_work(ind+j)=d_t_old(j,i+nres)
2406 d_a_work(ind+j)=d_a_old(j,i+nres)
2415 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2419 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2420 & vfric_mat(i,j),afric_mat(i,j),
2421 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2429 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2430 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2431 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2432 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2434 d_t_work_new(i)=ddt1+0.5d0*ddt2
2435 d_t_work(i)=ddt1+ddt2
2440 d_t(j,0)=d_t_work(j)
2445 dc(j,i)=dc_work(ind+j)
2446 d_t(j,i)=d_t_work(ind+j)
2451 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2454 dc(j,inres)=dc_work(ind+j)
2455 d_t(j,inres)=d_t_work(ind+j)
2462 c--------------------------------------------------------------------------
2463 subroutine sd_verlet2_ciccotti
2464 c Calculating the adjusted velocities for accelerations
2465 implicit real*8 (a-h,o-z)
2466 include 'DIMENSIONS'
2467 include 'COMMON.CONTROL'
2468 include 'COMMON.VAR'
2471 include 'COMMON.LANGEVIN'
2473 include 'COMMON.LANGEVIN.lang0'
2475 include 'COMMON.CHAIN'
2476 include 'COMMON.DERIV'
2477 include 'COMMON.GEO'
2478 include 'COMMON.LOCAL'
2479 include 'COMMON.INTERACT'
2480 include 'COMMON.IOUNITS'
2481 include 'COMMON.NAMES'
2482 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2483 common /stochcalc/ stochforcvec
2485 c Compute the stochastic forces which contribute to velocity change
2487 call stochastic_force(stochforcvecV)
2494 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2495 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2496 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2498 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2502 d_t(j,0)=d_t_work(j)
2507 d_t(j,i)=d_t_work(ind+j)
2512 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2515 d_t(j,inres)=d_t_work(ind+j)