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 write (iout,*) "PREMINIM ",preminim
1614 if(iranconf.ne.0 .or. preminim) then
1616 print *, 'Calling OVERLAP_SC'
1617 call overlap_sc(fail)
1621 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1622 print *,'SC_move',nft_sc,etot
1623 if(me.eq.king.or..not.out1file)
1624 & write(iout,*) 'SC_move',nft_sc,etot
1628 print *, 'Calling MINIM_DC'
1629 call minim_dc(etot,iretcode,nfun)
1631 call geom_to_var(nvar,varia)
1632 print *,'Calling MINIMIZE.'
1633 call minimize(etot,varia,iretcode,nfun)
1634 call var_to_geom(nvar,varia)
1636 if(me.eq.king.or..not.out1file) then
1637 write(iout,*) "Minimized energy is",etot
1638 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1639 call etotal(potEcomp)
1640 call enerprint(potEcomp)
1644 call chainbuild_cart
1649 kinetic_T=2.0d0/(dimen3*Rb)*EK
1650 if(me.eq.king.or..not.out1file)then
1660 call etotal(potEcomp)
1661 if (large) call enerprint(potEcomp)
1664 t_etotal=t_etotal+MPI_Wtime()-tt0
1666 t_etotal=t_etotal+tcpu()-tt0
1673 if (amax*d_time .gt. dvmax) then
1674 d_time=d_time*dvmax/amax
1675 if(me.eq.king.or..not.out1file) write (iout,*)
1676 & "Time step reduced to",d_time,
1677 & " because of too large initial acceleration."
1679 C if(me.eq.king.or..not.out1file)then
1680 C write(iout,*) "Potential energy and its components"
1681 C call enerprint(potEcomp)
1682 c write(iout,*) (potEcomp(i),i=0,n_ene)
1684 potE=potEcomp(0)-potEcomp(20)
1687 if (ntwe.ne.0) call statout(itime)
1688 if(me.eq.king.or..not.out1file)
1689 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1690 & " Kinetic energy",EK," Potential energy",potE,
1691 & " Total energy",totE," Maximum acceleration ",
1694 write (iout,*) "Initial coordinates"
1696 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1697 & (c(j,i+nres),j=1,3)
1699 write (iout,*) "Initial dC"
1701 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1702 & (dc(j,i+nres),j=1,3)
1704 write (iout,*) "Initial velocities"
1705 write (iout,"(13x,' backbone ',23x,' side chain')")
1707 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1708 & (d_t(j,i+nres),j=1,3)
1710 write (iout,*) "Initial accelerations"
1712 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1713 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1714 & (d_a(j,i+nres),j=1,3)
1720 d_t_old(j,i)=d_t(j,i)
1721 d_a_old(j,i)=d_a(j,i)
1723 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1732 call etotal_short(energia_short)
1733 if (large) call enerprint(potEcomp)
1736 t_eshort=t_eshort+MPI_Wtime()-tt0
1738 t_eshort=t_eshort+tcpu()-tt0
1743 if(.not.out1file .and. large) then
1744 write (iout,*) "energia_long",energia_long(0),
1745 & " energia_short",energia_short(0),
1746 & " total",energia_long(0)+energia_short(0)
1747 write (iout,*) "Initial fast-force accelerations"
1749 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1750 & (d_a(j,i+nres),j=1,3)
1753 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1756 d_a_short(j,i)=d_a(j,i)
1765 call etotal_long(energia_long)
1766 if (large) call enerprint(potEcomp)
1769 t_elong=t_elong+MPI_Wtime()-tt0
1771 t_elong=t_elong+tcpu()-tt0
1776 if(.not.out1file .and. large) then
1777 write (iout,*) "energia_long",energia_long(0)
1778 write (iout,*) "Initial slow-force accelerations"
1780 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1781 & (d_a(j,i+nres),j=1,3)
1785 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1787 t_enegrad=t_enegrad+tcpu()-tt0
1792 c-----------------------------------------------------------
1793 subroutine random_vel
1794 implicit real*8 (a-h,o-z)
1795 include 'DIMENSIONS'
1796 include 'COMMON.CONTROL'
1797 include 'COMMON.VAR'
1800 include 'COMMON.LANGEVIN'
1802 include 'COMMON.LANGEVIN.lang0'
1804 include 'COMMON.CHAIN'
1805 include 'COMMON.DERIV'
1806 include 'COMMON.GEO'
1807 include 'COMMON.LOCAL'
1808 include 'COMMON.INTERACT'
1809 include 'COMMON.IOUNITS'
1810 include 'COMMON.NAMES'
1811 include 'COMMON.TIME1'
1812 double precision xv,sigv,lowb,highb,vec_afm(3)
1813 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1814 c First generate velocities in the eigenspace of the G matrix
1815 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1822 sigv=dsqrt((Rb*t_bath)/geigen(i))
1825 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1827 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1828 c & " d_t_work_new",d_t_work_new(ii)
1831 C if (SELFGUIDE.gt.0) then
1834 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
1835 C distance=distance+vec_afm(j)**2
1837 C distance=dsqrt(distance)
1839 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
1840 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
1841 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
1842 C & d_t_work_new(j+(afmend-1)*3)
1853 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1856 c write (iout,*) "Ek from eigenvectors",Ek1
1858 c Transform velocities to UNRES coordinate space
1864 d_t_work(ind)=d_t_work(ind)
1865 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1867 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1871 c Transfer to the d_t vector
1873 d_t(j,0)=d_t_work(j)
1879 d_t(j,i)=d_t_work(ind)
1883 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1886 d_t(j,i+nres)=d_t_work(ind)
1891 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1892 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1897 c-----------------------------------------------------------
1898 subroutine sd_verlet_p_setup
1899 c Sets up the parameters of stochastic Verlet algorithm
1900 implicit real*8 (a-h,o-z)
1901 include 'DIMENSIONS'
1905 include 'COMMON.CONTROL'
1906 include 'COMMON.VAR'
1909 include 'COMMON.LANGEVIN'
1911 include 'COMMON.LANGEVIN.lang0'
1913 include 'COMMON.CHAIN'
1914 include 'COMMON.DERIV'
1915 include 'COMMON.GEO'
1916 include 'COMMON.LOCAL'
1917 include 'COMMON.INTERACT'
1918 include 'COMMON.IOUNITS'
1919 include 'COMMON.NAMES'
1920 include 'COMMON.TIME1'
1921 double precision emgdt(MAXRES6),
1922 & pterm,vterm,rho,rhoc,vsig,
1923 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1924 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1925 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1926 logical lprn /.false./
1927 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1928 double precision ktm
1935 c AL 8/17/04 Code adapted from tinker
1937 c Get the frictional and random terms for stochastic dynamics in the
1938 c eigenspace of mass-scaled UNRES friction matrix
1941 gdt = fricgam(i) * d_time
1943 c Stochastic dynamics reduces to simple MD for zero friction
1945 if (gdt .le. zero) then
1946 pfric_vec(i) = 1.0d0
1947 vfric_vec(i) = d_time
1948 afric_vec(i) = 0.5d0 * d_time * d_time
1949 prand_vec(i) = 0.0d0
1950 vrand_vec1(i) = 0.0d0
1951 vrand_vec2(i) = 0.0d0
1953 c Analytical expressions when friction coefficient is large
1956 if (gdt .ge. gdt_radius) then
1959 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1960 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1961 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1962 vterm = 1.0d0 - egdt**2
1963 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1965 c Use series expansions when friction coefficient is small
1976 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1977 & - gdt5/120.0d0 + gdt6/720.0d0
1978 & - gdt7/5040.0d0 + gdt8/40320.0d0
1979 & - gdt9/362880.0d0) / fricgam(i)**2
1980 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1981 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1982 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1983 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1984 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1985 & + 127.0d0*gdt9/90720.0d0
1986 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1987 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1988 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1989 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1990 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1991 & - 17.0d0*gdt2/1280.0d0
1992 & + 17.0d0*gdt3/6144.0d0
1993 & + 40967.0d0*gdt4/34406400.0d0
1994 & - 57203.0d0*gdt5/275251200.0d0
1995 & - 1429487.0d0*gdt6/13212057600.0d0)
1998 c Compute the scaling factors of random terms for the nonzero friction case
2000 ktm = 0.5d0*d_time/fricgam(i)
2001 psig = dsqrt(ktm*pterm) / fricgam(i)
2002 vsig = dsqrt(ktm*vterm)
2003 rhoc = dsqrt(1.0d0 - rho*rho)
2005 vrand_vec1(i) = vsig * rho
2006 vrand_vec2(i) = vsig * rhoc
2011 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2014 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2015 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2019 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2022 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2023 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2024 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2025 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2026 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2027 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2030 t_sdsetup=t_sdsetup+MPI_Wtime()
2032 t_sdsetup=t_sdsetup+tcpu()-tt0
2036 c-------------------------------------------------------------
2037 subroutine eigtransf1(n,ndim,ab,d,c)
2040 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2046 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2052 c-------------------------------------------------------------
2053 subroutine eigtransf(n,ndim,a,b,d,c)
2056 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2062 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2068 c-------------------------------------------------------------
2069 subroutine sd_verlet1
2070 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2071 implicit real*8 (a-h,o-z)
2072 include 'DIMENSIONS'
2073 include 'COMMON.CONTROL'
2074 include 'COMMON.VAR'
2077 include 'COMMON.LANGEVIN'
2079 include 'COMMON.LANGEVIN.lang0'
2081 include 'COMMON.CHAIN'
2082 include 'COMMON.DERIV'
2083 include 'COMMON.GEO'
2084 include 'COMMON.LOCAL'
2085 include 'COMMON.INTERACT'
2086 include 'COMMON.IOUNITS'
2087 include 'COMMON.NAMES'
2088 double precision stochforcvec(MAXRES6)
2089 common /stochcalc/ stochforcvec
2090 logical lprn /.false./
2092 c write (iout,*) "dc_old"
2094 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2095 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2098 dc_work(j)=dc_old(j,0)
2099 d_t_work(j)=d_t_old(j,0)
2100 d_a_work(j)=d_a_old(j,0)
2105 dc_work(ind+j)=dc_old(j,i)
2106 d_t_work(ind+j)=d_t_old(j,i)
2107 d_a_work(ind+j)=d_a_old(j,i)
2112 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2114 dc_work(ind+j)=dc_old(j,i+nres)
2115 d_t_work(ind+j)=d_t_old(j,i+nres)
2116 d_a_work(ind+j)=d_a_old(j,i+nres)
2124 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2128 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2129 & vfric_mat(i,j),afric_mat(i,j),
2130 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2138 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2139 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2140 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2141 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2143 d_t_work_new(i)=ddt1+0.5d0*ddt2
2144 d_t_work(i)=ddt1+ddt2
2149 d_t(j,0)=d_t_work(j)
2154 dc(j,i)=dc_work(ind+j)
2155 d_t(j,i)=d_t_work(ind+j)
2160 if (itype(i).ne.10) then
2163 dc(j,inres)=dc_work(ind+j)
2164 d_t(j,inres)=d_t_work(ind+j)
2171 c--------------------------------------------------------------------------
2172 subroutine sd_verlet2
2173 c Calculating the adjusted velocities for accelerations
2174 implicit real*8 (a-h,o-z)
2175 include 'DIMENSIONS'
2176 include 'COMMON.CONTROL'
2177 include 'COMMON.VAR'
2180 include 'COMMON.LANGEVIN'
2182 include 'COMMON.LANGEVIN.lang0'
2184 include 'COMMON.CHAIN'
2185 include 'COMMON.DERIV'
2186 include 'COMMON.GEO'
2187 include 'COMMON.LOCAL'
2188 include 'COMMON.INTERACT'
2189 include 'COMMON.IOUNITS'
2190 include 'COMMON.NAMES'
2191 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2192 common /stochcalc/ stochforcvec
2194 c Compute the stochastic forces which contribute to velocity change
2196 call stochastic_force(stochforcvecV)
2203 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2204 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2205 & vrand_mat2(i,j)*stochforcvecV(j)
2207 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2211 d_t(j,0)=d_t_work(j)
2216 d_t(j,i)=d_t_work(ind+j)
2221 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2224 d_t(j,inres)=d_t_work(ind+j)
2231 c-----------------------------------------------------------
2232 subroutine sd_verlet_ciccotti_setup
2233 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2235 implicit real*8 (a-h,o-z)
2236 include 'DIMENSIONS'
2240 include 'COMMON.CONTROL'
2241 include 'COMMON.VAR'
2244 include 'COMMON.LANGEVIN'
2246 include 'COMMON.LANGEVIN.lang0'
2248 include 'COMMON.CHAIN'
2249 include 'COMMON.DERIV'
2250 include 'COMMON.GEO'
2251 include 'COMMON.LOCAL'
2252 include 'COMMON.INTERACT'
2253 include 'COMMON.IOUNITS'
2254 include 'COMMON.NAMES'
2255 include 'COMMON.TIME1'
2256 double precision emgdt(MAXRES6),
2257 & pterm,vterm,rho,rhoc,vsig,
2258 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2259 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2260 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2261 logical lprn /.false./
2262 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2263 double precision ktm
2270 c AL 8/17/04 Code adapted from tinker
2272 c Get the frictional and random terms for stochastic dynamics in the
2273 c eigenspace of mass-scaled UNRES friction matrix
2276 write (iout,*) "i",i," fricgam",fricgam(i)
2277 gdt = fricgam(i) * d_time
2279 c Stochastic dynamics reduces to simple MD for zero friction
2281 if (gdt .le. zero) then
2282 pfric_vec(i) = 1.0d0
2283 vfric_vec(i) = d_time
2284 afric_vec(i) = 0.5d0*d_time*d_time
2285 prand_vec(i) = afric_vec(i)
2286 vrand_vec2(i) = vfric_vec(i)
2288 c Analytical expressions when friction coefficient is large
2293 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2294 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2295 prand_vec(i) = afric_vec(i)
2296 vrand_vec2(i) = vfric_vec(i)
2298 c Compute the scaling factors of random terms for the nonzero friction case
2300 c ktm = 0.5d0*d_time/fricgam(i)
2301 c psig = dsqrt(ktm*pterm) / fricgam(i)
2302 c vsig = dsqrt(ktm*vterm)
2303 c prand_vec(i) = psig*afric_vec(i)
2304 c vrand_vec2(i) = vsig*vfric_vec(i)
2309 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2312 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2313 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2317 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2319 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2320 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2321 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2322 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2323 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2325 t_sdsetup=t_sdsetup+MPI_Wtime()
2327 t_sdsetup=t_sdsetup+tcpu()-tt0
2331 c-------------------------------------------------------------
2332 subroutine sd_verlet1_ciccotti
2333 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2334 implicit real*8 (a-h,o-z)
2335 include 'DIMENSIONS'
2339 include 'COMMON.CONTROL'
2340 include 'COMMON.VAR'
2343 include 'COMMON.LANGEVIN'
2345 include 'COMMON.LANGEVIN.lang0'
2347 include 'COMMON.CHAIN'
2348 include 'COMMON.DERIV'
2349 include 'COMMON.GEO'
2350 include 'COMMON.LOCAL'
2351 include 'COMMON.INTERACT'
2352 include 'COMMON.IOUNITS'
2353 include 'COMMON.NAMES'
2354 double precision stochforcvec(MAXRES6)
2355 common /stochcalc/ stochforcvec
2356 logical lprn /.false./
2358 c write (iout,*) "dc_old"
2360 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2361 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2364 dc_work(j)=dc_old(j,0)
2365 d_t_work(j)=d_t_old(j,0)
2366 d_a_work(j)=d_a_old(j,0)
2371 dc_work(ind+j)=dc_old(j,i)
2372 d_t_work(ind+j)=d_t_old(j,i)
2373 d_a_work(ind+j)=d_a_old(j,i)
2378 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2380 dc_work(ind+j)=dc_old(j,i+nres)
2381 d_t_work(ind+j)=d_t_old(j,i+nres)
2382 d_a_work(ind+j)=d_a_old(j,i+nres)
2391 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2395 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2396 & vfric_mat(i,j),afric_mat(i,j),
2397 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2405 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2406 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2407 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2408 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2410 d_t_work_new(i)=ddt1+0.5d0*ddt2
2411 d_t_work(i)=ddt1+ddt2
2416 d_t(j,0)=d_t_work(j)
2421 dc(j,i)=dc_work(ind+j)
2422 d_t(j,i)=d_t_work(ind+j)
2427 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2430 dc(j,inres)=dc_work(ind+j)
2431 d_t(j,inres)=d_t_work(ind+j)
2438 c--------------------------------------------------------------------------
2439 subroutine sd_verlet2_ciccotti
2440 c Calculating the adjusted velocities for accelerations
2441 implicit real*8 (a-h,o-z)
2442 include 'DIMENSIONS'
2443 include 'COMMON.CONTROL'
2444 include 'COMMON.VAR'
2447 include 'COMMON.LANGEVIN'
2449 include 'COMMON.LANGEVIN.lang0'
2451 include 'COMMON.CHAIN'
2452 include 'COMMON.DERIV'
2453 include 'COMMON.GEO'
2454 include 'COMMON.LOCAL'
2455 include 'COMMON.INTERACT'
2456 include 'COMMON.IOUNITS'
2457 include 'COMMON.NAMES'
2458 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2459 common /stochcalc/ stochforcvec
2461 c Compute the stochastic forces which contribute to velocity change
2463 call stochastic_force(stochforcvecV)
2470 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2471 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2472 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2474 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2478 d_t(j,0)=d_t_work(j)
2483 d_t(j,i)=d_t_work(ind+j)
2488 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2491 d_t(j,inres)=d_t_work(ind+j)