2 c------------------------------------------------
3 c The driver for molecular dynamics subroutines
4 c------------------------------------------------
5 implicit real*8 (a-h,o-z)
11 include 'COMMON.SETUP'
12 include 'COMMON.CONTROL'
16 include 'COMMON.LANGEVIN'
18 include 'COMMON.LANGEVIN.lang0'
20 include 'COMMON.CHAIN'
21 include 'COMMON.DERIV'
23 include 'COMMON.LOCAL'
24 include 'COMMON.INTERACT'
25 include 'COMMON.IOUNITS'
26 include 'COMMON.NAMES'
27 include 'COMMON.TIME1'
28 include 'COMMON.HAIRPIN'
29 double precision cm(3),L(3),vcm(3)
31 double precision v_work(maxres6),v_transf(maxres6)
41 if (ilen(tmpdir).gt.0)
42 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
43 & //liczba(:ilen(liczba))//'.rst')
45 if (ilen(tmpdir).gt.0)
46 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
53 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
59 c Determine the inverse of the inertia matrix.
60 call setup_MD_matrices
64 t_MDsetup = MPI_Wtime()-tt0
66 t_MDsetup = tcpu()-tt0
69 c Entering the MD loop
75 if (lang.eq.2 .or. lang.eq.3) then
79 call sd_verlet_p_setup
81 call sd_verlet_ciccotti_setup
85 pfric0_mat(i,j,0)=pfric_mat(i,j)
86 afric0_mat(i,j,0)=afric_mat(i,j)
87 vfric0_mat(i,j,0)=vfric_mat(i,j)
88 prand0_mat(i,j,0)=prand_mat(i,j)
89 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
90 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
99 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
101 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
105 else if (lang.eq.1 .or. lang.eq.4) then
109 t_langsetup=MPI_Wtime()-tt0
112 t_langsetup=tcpu()-tt0
115 do itime=1,n_timestep
117 if (large.and. mod(itime,ntwe).eq.0)
118 & write (iout,*) "itime",itime
120 if (lang.gt.0 .and. surfarea .and.
121 & mod(itime,reset_fricmat).eq.0) then
122 if (lang.eq.2 .or. lang.eq.3) then
126 call sd_verlet_p_setup
128 call sd_verlet_ciccotti_setup
132 pfric0_mat(i,j,0)=pfric_mat(i,j)
133 afric0_mat(i,j,0)=afric_mat(i,j)
134 vfric0_mat(i,j,0)=vfric_mat(i,j)
135 prand0_mat(i,j,0)=prand_mat(i,j)
136 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
137 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
142 flag_stoch(i)=.false.
145 else if (lang.eq.1 .or. lang.eq.4) then
148 write (iout,'(a,i10)')
149 & "Friction matrix reset based on surface area, itime",itime
151 if (reset_vel .and. tbf .and. lang.eq.0
152 & .and. mod(itime,count_reset_vel).eq.0) then
154 write(iout,'(a,f20.2)')
155 & "Velocities reset to random values, time",totT
158 d_t_old(j,i)=d_t(j,i)
162 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
166 d_t(j,0)=d_t(j,0)-vcm(j)
169 kinetic_T=2.0d0/(dimen3*Rb)*EK
170 scalfac=dsqrt(T_bath/kinetic_T)
171 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
174 d_t_old(j,i)=scalfac*d_t(j,i)
180 c Time-reversible RESPA algorithm
181 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
182 call RESPA_step(itime)
184 c Variable time step algorithm.
185 call velverlet_step(itime)
189 call brown_step(itime)
191 print *,"Brown dynamics not here!"
193 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
199 if (mod(itime,ntwe).eq.0) then
201 C call enerprint(potEcomp)
202 C print *,itime,'AFM',Eafmforc,etot
216 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
219 v_work(ind)=d_t(j,i+nres)
224 write (66,'(80f10.5)')
225 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
229 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
231 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
233 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
236 if (mod(itime,ntwx).eq.0) then
237 write(iout,*) 'time=',itime
238 C call check_ecartint
240 write (tytul,'("time",f8.2)') totT
242 call hairpin(.true.,nharp,iharp)
243 call secondary2(.true.)
244 call pdbout(potE,tytul,ipdb)
249 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
250 open(irest2,file=rest2name,status='unknown')
251 write(irest2,*) totT,EK,potE,totE,t_bath
253 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
256 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
267 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
269 & 'MD calculations setup:',t_MDsetup,
270 & 'Energy & gradient evaluation:',t_enegrad,
271 & 'Stochastic MD setup:',t_langsetup,
272 & 'Stochastic MD step setup:',t_sdsetup,
274 write (iout,'(/28(1h=),a25,27(1h=))')
275 & ' End of MD calculation '
277 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
279 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
280 & " time_fricmatmult",time_fricmatmult," time_fsample ",
285 c-------------------------------------------------------------------------------
286 subroutine velverlet_step(itime)
287 c-------------------------------------------------------------------------------
288 c Perform a single velocity Verlet step; the time step can be rescaled if
289 c increments in accelerations exceed the threshold
290 c-------------------------------------------------------------------------------
291 implicit real*8 (a-h,o-z)
295 integer ierror,ierrcode
297 include 'COMMON.SETUP'
298 include 'COMMON.CONTROL'
302 include 'COMMON.LANGEVIN'
304 include 'COMMON.LANGEVIN.lang0'
306 include 'COMMON.CHAIN'
307 include 'COMMON.DERIV'
309 include 'COMMON.LOCAL'
310 include 'COMMON.INTERACT'
311 include 'COMMON.IOUNITS'
312 include 'COMMON.NAMES'
313 include 'COMMON.TIME1'
314 include 'COMMON.MUCA'
315 double precision vcm(3),incr(3)
316 double precision cm(3),L(3)
317 integer ilen,count,rstcount
320 integer maxcount_scale /20/
322 double precision stochforcvec(MAXRES6)
323 common /stochcalc/ stochforcvec
331 else if (lang.eq.2 .or. lang.eq.3) then
333 call stochastic_force(stochforcvec)
336 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
338 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
345 icount_scale=icount_scale+1
346 if (icount_scale.gt.maxcount_scale) then
348 & "ERROR: too many attempts at scaling down the time step. ",
349 & "amax=",amax,"epdrift=",epdrift,
350 & "damax=",damax,"edriftmax=",edriftmax,
354 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
358 c First step of the velocity Verlet algorithm
363 else if (lang.eq.3) then
365 call sd_verlet1_ciccotti
367 else if (lang.eq.1) then
372 c Build the chain from the newly calculated coordinates
374 if (rattle) call rattle1
376 if (large.and. mod(itime,ntwe).eq.0) then
377 write (iout,*) "Cartesian and internal coordinates: step 1"
382 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
383 & (dc(j,i+nres),j=1,3)
385 write (iout,*) "Accelerations"
387 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
388 & (d_a(j,i+nres),j=1,3)
390 write (iout,*) "Velocities, step 1"
392 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
393 & (d_t(j,i+nres),j=1,3)
402 c Calculate energy and forces
404 call etotal(potEcomp)
405 if (large.and. mod(itime,ntwe).eq.0)
406 & call enerprint(potEcomp)
409 t_etotal=t_etotal+MPI_Wtime()-tt0
411 t_etotal=t_etotal+tcpu()-tt0
414 potE=potEcomp(0)-potEcomp(20)
416 c Get the new accelerations
419 t_enegrad=t_enegrad+MPI_Wtime()-tt0
421 t_enegrad=t_enegrad+tcpu()-tt0
423 c Determine maximum acceleration and scale down the timestep if needed
425 amax=amax/(itime_scal+1)**2
426 call predict_edrift(epdrift)
427 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
428 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
430 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
432 itime_scal=itime_scal+ifac_time
433 c fac_time=dmin1(damax/amax,0.5d0)
434 fac_time=0.5d0**ifac_time
435 d_time=d_time*fac_time
436 if (lang.eq.2 .or. lang.eq.3) then
438 c write (iout,*) "Calling sd_verlet_setup: 1"
439 c Rescale the stochastic forces and recalculate or restore
440 c the matrices of tinker integrator
441 if (itime_scal.gt.maxflag_stoch) then
442 if (large) write (iout,'(a,i5,a)')
443 & "Calculate matrices for stochastic step;",
444 & " itime_scal ",itime_scal
446 call sd_verlet_p_setup
448 call sd_verlet_ciccotti_setup
450 write (iout,'(2a,i3,a,i3,1h.)')
451 & "Warning: cannot store matrices for stochastic",
452 & " integration because the index",itime_scal,
453 & " is greater than",maxflag_stoch
454 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
455 & " integration Langevin algorithm for better efficiency."
456 else if (flag_stoch(itime_scal)) then
457 if (large) write (iout,'(a,i5,a,l1)')
458 & "Restore matrices for stochastic step; itime_scal ",
459 & itime_scal," flag ",flag_stoch(itime_scal)
462 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
463 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
464 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
465 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
466 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
467 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
471 if (large) write (iout,'(2a,i5,a,l1)')
472 & "Calculate & store matrices for stochastic step;",
473 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
475 call sd_verlet_p_setup
477 call sd_verlet_ciccotti_setup
479 flag_stoch(ifac_time)=.true.
482 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
483 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
484 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
485 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
486 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
487 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
491 fac_time=1.0d0/dsqrt(fac_time)
493 stochforcvec(i)=fac_time*stochforcvec(i)
496 else if (lang.eq.1) then
497 c Rescale the accelerations due to stochastic forces
498 fac_time=1.0d0/dsqrt(fac_time)
500 d_as_work(i)=d_as_work(i)*fac_time
503 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
504 & "itime",itime," Timestep scaled down to ",
505 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
507 c Second step of the velocity Verlet algorithm
512 else if (lang.eq.3) then
514 call sd_verlet2_ciccotti
516 else if (lang.eq.1) then
521 if (rattle) call rattle2
524 C print *,totTafm,"TU?"
525 if (d_time.ne.d_time0) then
528 if (lang.eq.2 .or. lang.eq.3) then
529 if (large) write (iout,'(a)')
530 & "Restore original matrices for stochastic step"
531 c write (iout,*) "Calling sd_verlet_setup: 2"
532 c Restore the matrices of tinker integrator if the time step has been restored
535 pfric_mat(i,j)=pfric0_mat(i,j,0)
536 afric_mat(i,j)=afric0_mat(i,j,0)
537 vfric_mat(i,j)=vfric0_mat(i,j,0)
538 prand_mat(i,j)=prand0_mat(i,j,0)
539 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
540 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
549 c Calculate the kinetic and the total energy and the kinetic temperature
554 c write (iout,*) "step",itime," EK",EK," EK1",EK1
556 c Couple the system to Berendsen bath if needed
557 if (tbf .and. lang.eq.0) then
560 kinetic_T=2.0d0/(dimen3*Rb)*EK
561 c Backup the coordinates, velocities, and accelerations
565 d_t_old(j,i)=d_t(j,i)
566 d_a_old(j,i)=d_a(j,i)
570 if (mod(itime,ntwe).eq.0 .and. large) then
571 write (iout,*) "Velocities, step 2"
573 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
574 & (d_t(j,i+nres),j=1,3)
580 c-------------------------------------------------------------------------------
581 subroutine RESPA_step(itime)
582 c-------------------------------------------------------------------------------
583 c Perform a single RESPA step.
584 c-------------------------------------------------------------------------------
585 implicit real*8 (a-h,o-z)
589 integer IERROR,ERRCODE
591 include 'COMMON.SETUP'
592 include 'COMMON.CONTROL'
596 include 'COMMON.LANGEVIN'
598 include 'COMMON.LANGEVIN.lang0'
600 include 'COMMON.CHAIN'
601 include 'COMMON.DERIV'
603 include 'COMMON.LOCAL'
604 include 'COMMON.INTERACT'
605 include 'COMMON.IOUNITS'
606 include 'COMMON.NAMES'
607 include 'COMMON.TIME1'
608 double precision energia_short(0:n_ene),
609 & energia_long(0:n_ene)
610 double precision cm(3),L(3),vcm(3),incr(3)
611 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
612 & d_a_old0(3,0:maxres2)
613 logical PRINT_AMTS_MSG /.false./
614 integer ilen,count,rstcount
617 integer maxcount_scale /10/
619 double precision stochforcvec(MAXRES6)
620 common /stochcalc/ stochforcvec
623 common /cipiszcze/ itt
626 if (large.and. mod(itime,ntwe).eq.0) then
627 write (iout,*) "***************** RESPA itime",itime
628 write (iout,*) "Cartesian and internal coordinates: step 0"
630 call pdbout(0.0d0,"cipiszcze",iout)
632 write (iout,*) "Accelerations from long-range forces"
634 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
635 & (d_a(j,i+nres),j=1,3)
637 write (iout,*) "Velocities, step 0"
639 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
640 & (d_t(j,i+nres),j=1,3)
645 c Perform the initial RESPA step (increment velocities)
646 c write (iout,*) "*********************** RESPA ini"
649 if (mod(itime,ntwe).eq.0 .and. large) then
650 write (iout,*) "Velocities, end"
652 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
653 & (d_t(j,i+nres),j=1,3)
657 c Compute the short-range forces
663 C 7/2/2009 commented out
665 c call etotal_short(energia_short)
668 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
671 d_a(j,i)=d_a_short(j,i)
675 if (large.and. mod(itime,ntwe).eq.0) then
676 write (iout,*) "energia_short",energia_short(0)
677 write (iout,*) "Accelerations from short-range forces"
679 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
680 & (d_a(j,i+nres),j=1,3)
685 t_enegrad=t_enegrad+MPI_Wtime()-tt0
687 t_enegrad=t_enegrad+tcpu()-tt0
692 d_t_old(j,i)=d_t(j,i)
693 d_a_old(j,i)=d_a(j,i)
696 c 6/30/08 A-MTS: attempt at increasing the split number
699 dc_old0(j,i)=dc_old(j,i)
700 d_t_old0(j,i)=d_t_old(j,i)
701 d_a_old0(j,i)=d_a_old(j,i)
704 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
705 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
712 c write (iout,*) "itime",itime," ntime_split",ntime_split
713 c Split the time step
714 d_time=d_time0/ntime_split
715 c Perform the short-range RESPA steps (velocity Verlet increments of
716 c positions and velocities using short-range forces)
717 c write (iout,*) "*********************** RESPA split"
718 do itsplit=1,ntime_split
721 else if (lang.eq.2 .or. lang.eq.3) then
723 call stochastic_force(stochforcvec)
726 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
728 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
733 c First step of the velocity Verlet algorithm
738 else if (lang.eq.3) then
740 call sd_verlet1_ciccotti
742 else if (lang.eq.1) then
747 c Build the chain from the newly calculated coordinates
749 if (rattle) call rattle1
751 if (large.and. mod(itime,ntwe).eq.0) then
752 write (iout,*) "***** ITSPLIT",itsplit
753 write (iout,*) "Cartesian and internal coordinates: step 1"
754 call pdbout(0.0d0,"cipiszcze",iout)
757 write (iout,*) "Velocities, step 1"
759 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
760 & (d_t(j,i+nres),j=1,3)
769 c Calculate energy and forces
771 call etotal_short(energia_short)
772 if (large.and. mod(itime,ntwe).eq.0)
773 & call enerprint(energia_short)
776 t_eshort=t_eshort+MPI_Wtime()-tt0
778 t_eshort=t_eshort+tcpu()-tt0
782 c Get the new accelerations
784 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
787 d_a_short(j,i)=d_a(j,i)
791 if (large.and. mod(itime,ntwe).eq.0) then
792 write (iout,*)"energia_short",energia_short(0)
793 write (iout,*) "Accelerations from short-range forces"
795 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
796 & (d_a(j,i+nres),j=1,3)
801 c Determine maximum acceleration and scale down the timestep if needed
803 amax=amax/ntime_split**2
804 call predict_edrift(epdrift)
805 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
806 & write (iout,*) "amax",amax," damax",damax,
807 & " epdrift",epdrift," epdriftmax",epdriftmax
808 c Exit loop and try with increased split number if the change of
809 c acceleration is too big
810 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
811 if (ntime_split.lt.maxtime_split) then
813 ntime_split=ntime_split*2
816 dc_old(j,i)=dc_old0(j,i)
817 d_t_old(j,i)=d_t_old0(j,i)
818 d_a_old(j,i)=d_a_old0(j,i)
821 if (PRINT_AMTS_MSG) then
822 write (iout,*) "acceleration/energy drift too large",amax,
823 & epdrift," split increased to ",ntime_split," itime",itime,
829 & "Uh-hu. Bumpy landscape. Maximum splitting number",
831 & " already reached!!! Trying to carry on!"
835 t_enegrad=t_enegrad+MPI_Wtime()-tt0
837 t_enegrad=t_enegrad+tcpu()-tt0
839 c Second step of the velocity Verlet algorithm
844 else if (lang.eq.3) then
846 call sd_verlet2_ciccotti
848 else if (lang.eq.1) then
853 if (rattle) call rattle2
854 c Backup the coordinates, velocities, and accelerations
858 d_t_old(j,i)=d_t(j,i)
859 d_a_old(j,i)=d_a(j,i)
866 c Restore the time step
868 c Compute long-range forces
875 call etotal_long(energia_long)
876 if (large.and. mod(itime,ntwe).eq.0)
877 & call enerprint(energia_long)
880 t_elong=t_elong+MPI_Wtime()-tt0
882 t_elong=t_elong+tcpu()-tt0
888 t_enegrad=t_enegrad+MPI_Wtime()-tt0
890 t_enegrad=t_enegrad+tcpu()-tt0
892 c Compute accelerations from long-range forces
894 if (large.and. mod(itime,ntwe).eq.0) then
895 write (iout,*) "energia_long",energia_long(0)
896 write (iout,*) "Cartesian and internal coordinates: step 2"
898 call pdbout(0.0d0,"cipiszcze",iout)
900 write (iout,*) "Accelerations from long-range forces"
902 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
903 & (d_a(j,i+nres),j=1,3)
905 write (iout,*) "Velocities, step 2"
907 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
908 & (d_t(j,i+nres),j=1,3)
912 c Compute the final RESPA step (increment velocities)
913 c write (iout,*) "*********************** RESPA fin"
915 c Compute the complete potential energy
917 potEcomp(i)=energia_short(i)+energia_long(i)
919 potE=potEcomp(0)-potEcomp(20)
920 c potE=energia_short(0)+energia_long(0)
923 c Calculate the kinetic and the total energy and the kinetic temperature
926 c Couple the system to Berendsen bath if needed
927 if (tbf .and. lang.eq.0) then
930 kinetic_T=2.0d0/(dimen3*Rb)*EK
931 c Backup the coordinates, velocities, and accelerations
933 if (mod(itime,ntwe).eq.0 .and. large) then
934 write (iout,*) "Velocities, end"
936 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
937 & (d_t(j,i+nres),j=1,3)
943 c---------------------------------------------------------------------
945 c First and last RESPA step (incrementing velocities using long-range
947 implicit real*8 (a-h,o-z)
949 include 'COMMON.CONTROL'
952 include 'COMMON.CHAIN'
953 include 'COMMON.DERIV'
955 include 'COMMON.LOCAL'
956 include 'COMMON.INTERACT'
957 include 'COMMON.IOUNITS'
958 include 'COMMON.NAMES'
960 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
964 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
968 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
971 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
977 c-----------------------------------------------------------------
979 c Applying velocity Verlet algorithm - step 1 to coordinates
980 implicit real*8 (a-h,o-z)
982 include 'COMMON.CONTROL'
985 include 'COMMON.CHAIN'
986 include 'COMMON.DERIV'
988 include 'COMMON.LOCAL'
989 include 'COMMON.INTERACT'
990 include 'COMMON.IOUNITS'
991 include 'COMMON.NAMES'
992 double precision adt,adt2
995 write (iout,*) "VELVERLET1 START: DC"
997 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
998 & (dc(j,i+nres),j=1,3)
1002 adt=d_a_old(j,0)*d_time
1004 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1005 d_t_new(j,0)=d_t_old(j,0)+adt2
1006 d_t(j,0)=d_t_old(j,0)+adt
1012 adt=d_a_old(j,i)*d_time
1014 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1015 d_t_new(j,i)=d_t_old(j,i)+adt2
1016 d_t(j,i)=d_t_old(j,i)+adt
1021 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1024 adt=d_a_old(j,inres)*d_time
1026 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1027 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1028 d_t(j,inres)=d_t_old(j,inres)+adt
1033 write (iout,*) "VELVERLET1 END: DC"
1035 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1036 & (dc(j,i+nres),j=1,3)
1041 c---------------------------------------------------------------------
1043 c Step 2 of the velocity Verlet algorithm: update velocities
1044 implicit real*8 (a-h,o-z)
1045 include 'DIMENSIONS'
1046 include 'COMMON.CONTROL'
1047 include 'COMMON.VAR'
1049 include 'COMMON.CHAIN'
1050 include 'COMMON.DERIV'
1051 include 'COMMON.GEO'
1052 include 'COMMON.LOCAL'
1053 include 'COMMON.INTERACT'
1054 include 'COMMON.IOUNITS'
1055 include 'COMMON.NAMES'
1057 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1061 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1065 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1068 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1074 c-----------------------------------------------------------------
1075 subroutine sddir_precalc
1076 c Applying velocity Verlet algorithm - step 1 to coordinates
1077 implicit real*8 (a-h,o-z)
1078 include 'DIMENSIONS'
1082 include 'COMMON.CONTROL'
1083 include 'COMMON.VAR'
1086 include 'COMMON.LANGEVIN'
1088 include 'COMMON.LANGEVIN.lang0'
1090 include 'COMMON.CHAIN'
1091 include 'COMMON.DERIV'
1092 include 'COMMON.GEO'
1093 include 'COMMON.LOCAL'
1094 include 'COMMON.INTERACT'
1095 include 'COMMON.IOUNITS'
1096 include 'COMMON.NAMES'
1097 include 'COMMON.TIME1'
1098 double precision stochforcvec(MAXRES6)
1099 common /stochcalc/ stochforcvec
1101 c Compute friction and stochastic forces
1110 time_fric=time_fric+MPI_Wtime()-time00
1113 time_fric=time_fric+tcpu()-time00
1116 call stochastic_force(stochforcvec)
1118 time_stoch=time_stoch+MPI_Wtime()-time00
1120 time_stoch=time_stoch+tcpu()-time00
1123 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1124 c forces (d_as_work)
1126 call ginv_mult(fric_work, d_af_work)
1127 call ginv_mult(stochforcvec, d_as_work)
1130 c---------------------------------------------------------------------
1131 subroutine sddir_verlet1
1132 c Applying velocity Verlet algorithm - step 1 to velocities
1133 implicit real*8 (a-h,o-z)
1134 include 'DIMENSIONS'
1135 include 'COMMON.CONTROL'
1136 include 'COMMON.VAR'
1139 include 'COMMON.LANGEVIN'
1141 include 'COMMON.LANGEVIN.lang0'
1143 include 'COMMON.CHAIN'
1144 include 'COMMON.DERIV'
1145 include 'COMMON.GEO'
1146 include 'COMMON.LOCAL'
1147 include 'COMMON.INTERACT'
1148 include 'COMMON.IOUNITS'
1149 include 'COMMON.NAMES'
1150 c Revised 3/31/05 AL: correlation between random contributions to
1151 c position and velocity increments included.
1152 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1153 double precision adt,adt2
1155 c Add the contribution from BOTH friction and stochastic force to the
1156 c coordinates, but ONLY the contribution from the friction forces to velocities
1159 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1160 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1161 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1162 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1163 d_t(j,0)=d_t_old(j,0)+adt
1168 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1169 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1170 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1171 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1172 d_t(j,i)=d_t_old(j,i)+adt
1177 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1180 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1181 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1182 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1183 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1184 d_t(j,inres)=d_t_old(j,inres)+adt
1191 c---------------------------------------------------------------------
1192 subroutine sddir_verlet2
1193 c Calculating the adjusted velocities for accelerations
1194 implicit real*8 (a-h,o-z)
1195 include 'DIMENSIONS'
1196 include 'COMMON.CONTROL'
1197 include 'COMMON.VAR'
1200 include 'COMMON.LANGEVIN'
1202 include 'COMMON.LANGEVIN.lang0'
1204 include 'COMMON.CHAIN'
1205 include 'COMMON.DERIV'
1206 include 'COMMON.GEO'
1207 include 'COMMON.LOCAL'
1208 include 'COMMON.INTERACT'
1209 include 'COMMON.IOUNITS'
1210 include 'COMMON.NAMES'
1211 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1212 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1213 c Revised 3/31/05 AL: correlation between random contributions to
1214 c position and velocity increments included.
1215 c The correlation coefficients are calculated at low-friction limit.
1216 c Also, friction forces are now not calculated with new velocities.
1218 c call friction_force
1219 call stochastic_force(stochforcvec)
1221 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1222 c forces (d_as_work)
1224 call ginv_mult(stochforcvec, d_as_work1)
1230 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1231 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1236 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1237 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1242 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1245 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1246 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1247 & +cos60*d_as_work1(ind+j))*d_time
1254 c---------------------------------------------------------------------
1255 subroutine max_accel
1257 c Find the maximum difference in the accelerations of the the sites
1258 c at the beginning and the end of the time step.
1260 implicit real*8 (a-h,o-z)
1261 include 'DIMENSIONS'
1262 include 'COMMON.CONTROL'
1263 include 'COMMON.VAR'
1265 include 'COMMON.CHAIN'
1266 include 'COMMON.DERIV'
1267 include 'COMMON.GEO'
1268 include 'COMMON.LOCAL'
1269 include 'COMMON.INTERACT'
1270 include 'COMMON.IOUNITS'
1271 double precision aux(3),accel(3),accel_old(3),dacc
1273 c aux(j)=d_a(j,0)-d_a_old(j,0)
1274 accel_old(j)=d_a_old(j,0)
1281 c 7/3/08 changed to asymmetric difference
1283 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1284 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1285 accel(j)=accel(j)+0.5d0*d_a(j,i)
1286 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1287 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1288 dacc=dabs(accel(j)-accel_old(j))
1289 c write (iout,*) i,dacc
1290 if (dacc.gt.amax) amax=dacc
1298 accel_old(j)=d_a_old(j,0)
1303 accel_old(j)=accel_old(j)+d_a_old(j,1)
1304 accel(j)=accel(j)+d_a(j,1)
1308 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1310 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1311 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1312 accel(j)=accel(j)+d_a(j,i+nres)
1316 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1317 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1318 dacc=dabs(accel(j)-accel_old(j))
1319 c write (iout,*) "side-chain",i,dacc
1320 if (dacc.gt.amax) amax=dacc
1324 accel_old(j)=accel_old(j)+d_a_old(j,i)
1325 accel(j)=accel(j)+d_a(j,i)
1326 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1331 c---------------------------------------------------------------------
1332 subroutine predict_edrift(epdrift)
1334 c Predict the drift of the potential energy
1336 implicit real*8 (a-h,o-z)
1337 include 'DIMENSIONS'
1338 include 'COMMON.CONTROL'
1339 include 'COMMON.VAR'
1341 include 'COMMON.CHAIN'
1342 include 'COMMON.DERIV'
1343 include 'COMMON.GEO'
1344 include 'COMMON.LOCAL'
1345 include 'COMMON.INTERACT'
1346 include 'COMMON.IOUNITS'
1347 include 'COMMON.MUCA'
1348 double precision epdrift,epdriftij
1349 c Drift of the potential energy
1355 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1356 if (lmuca) epdriftij=epdriftij*factor
1357 c write (iout,*) "back",i,j,epdriftij
1358 if (epdriftij.gt.epdrift) epdrift=epdriftij
1362 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1365 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1366 if (lmuca) epdriftij=epdriftij*factor
1367 c write (iout,*) "side",i,j,epdriftij
1368 if (epdriftij.gt.epdrift) epdrift=epdriftij
1372 epdrift=0.5d0*epdrift*d_time*d_time
1373 c write (iout,*) "epdrift",epdrift
1376 c-----------------------------------------------------------------------
1377 subroutine verlet_bath
1379 c Coupling to the thermostat by using the Berendsen algorithm
1381 implicit real*8 (a-h,o-z)
1382 include 'DIMENSIONS'
1383 include 'COMMON.CONTROL'
1384 include 'COMMON.VAR'
1386 include 'COMMON.CHAIN'
1387 include 'COMMON.DERIV'
1388 include 'COMMON.GEO'
1389 include 'COMMON.LOCAL'
1390 include 'COMMON.INTERACT'
1391 include 'COMMON.IOUNITS'
1392 include 'COMMON.NAMES'
1393 double precision T_half,fact
1395 T_half=2.0d0/(dimen3*Rb)*EK
1396 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1397 c write(iout,*) "T_half", T_half
1398 c write(iout,*) "EK", EK
1399 c write(iout,*) "fact", fact
1401 d_t(j,0)=fact*d_t(j,0)
1405 d_t(j,i)=fact*d_t(j,i)
1409 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1412 d_t(j,inres)=fact*d_t(j,inres)
1418 c---------------------------------------------------------
1420 c Set up the initial conditions of a MD simulation
1421 implicit real*8 (a-h,o-z)
1422 include 'DIMENSIONS'
1426 integer IERROR,ERRCODE
1428 include 'COMMON.SETUP'
1429 include 'COMMON.CONTROL'
1430 include 'COMMON.VAR'
1433 include 'COMMON.LANGEVIN'
1435 include 'COMMON.LANGEVIN.lang0'
1437 include 'COMMON.CHAIN'
1438 include 'COMMON.DERIV'
1439 include 'COMMON.GEO'
1440 include 'COMMON.LOCAL'
1441 include 'COMMON.INTERACT'
1442 include 'COMMON.IOUNITS'
1443 include 'COMMON.NAMES'
1444 include 'COMMON.REMD'
1445 real*8 energia_long(0:n_ene),
1446 & energia_short(0:n_ene),vcm(3),incr(3)
1447 double precision cm(3),L(3),xv,sigv,lowb,highb
1448 double precision varia(maxvar)
1456 c write(iout,*) "d_time", d_time
1457 c Compute the standard deviations of stochastic forces for Langevin dynamics
1458 c if the friction coefficients do not depend on surface area
1459 if (lang.gt.0 .and. .not.surfarea) then
1461 stdforcp(i)=stdfp*dsqrt(gamp)
1464 stdforcsc(i)=stdfsc(iabs(itype(i)))
1465 & *dsqrt(gamsc(iabs(itype(i))))
1468 c Open the pdb file for snapshotshots
1471 if (ilen(tmpdir).gt.0)
1472 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1473 & liczba(:ilen(liczba))//".pdb")
1475 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1479 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1480 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1481 & liczba(:ilen(liczba))//".x")
1482 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1485 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1486 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1487 & liczba(:ilen(liczba))//".cx")
1488 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1494 if (ilen(tmpdir).gt.0)
1495 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1496 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1498 if (ilen(tmpdir).gt.0)
1499 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1500 cartname=prefix(:ilen(prefix))//"_MD.cx"
1504 write (qstr,'(256(1h ))')
1507 iq = qinfrag(i,iset)*10
1508 iw = wfrag(i,iset)/100
1510 if(me.eq.king.or..not.out1file)
1511 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1512 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1517 iq = qinpair(i,iset)*10
1518 iw = wpair(i,iset)/100
1520 if(me.eq.king.or..not.out1file)
1521 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1522 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1526 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1528 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1530 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1532 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1536 if (restart1file) then
1538 & inquire(file=mremd_rst_name,exist=file_exist)
1540 write (*,*) me," Before broadcast: file_exist",file_exist
1541 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1543 write (*,*) me," After broadcast: file_exist",file_exist
1544 c inquire(file=mremd_rst_name,exist=file_exist)
1546 if(me.eq.king.or..not.out1file)
1547 & write(iout,*) "Initial state read by master and distributed"
1549 if (ilen(tmpdir).gt.0)
1550 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1551 & //liczba(:ilen(liczba))//'.rst')
1552 inquire(file=rest2name,exist=file_exist)
1555 if(.not.restart1file) then
1556 if(me.eq.king.or..not.out1file)
1557 & write(iout,*) "Initial state will be read from file ",
1558 & rest2name(:ilen(rest2name))
1561 call rescale_weights(t_bath)
1563 if(me.eq.king.or..not.out1file)then
1564 if (restart1file) then
1565 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1568 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1571 write(iout,*) "Initial velocities randomly generated"
1578 c Generate initial velocities
1579 if(me.eq.king.or..not.out1file)
1580 & write(iout,*) "Initial velocities randomly generated"
1583 CtotTafm is the variable for AFM time which eclipsed during
1586 c rest2name = prefix(:ilen(prefix))//'.rst'
1587 if(me.eq.king.or..not.out1file)then
1588 write (iout,*) "Initial velocities"
1590 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1591 & (d_t(j,i+nres),j=1,3)
1593 c Zeroing the total angular momentum of the system
1594 write(iout,*) "Calling the zero-angular
1595 & momentum subroutine"
1598 c Getting the potential energy and forces and velocities and accelerations
1600 c write (iout,*) "velocity of the center of the mass:"
1601 c write (iout,*) (vcm(j),j=1,3)
1603 d_t(j,0)=d_t(j,0)-vcm(j)
1605 c Removing the velocity of the center of mass
1607 if(me.eq.king.or..not.out1file)then
1608 write (iout,*) "vcm right after adjustment:"
1609 write (iout,*) (vcm(j),j=1,3)
1613 if(iranconf.ne.0) then
1615 print *, 'Calling OVERLAP_SC'
1616 call overlap_sc(fail)
1620 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1621 print *,'SC_move',nft_sc,etot
1622 if(me.eq.king.or..not.out1file)
1623 & write(iout,*) 'SC_move',nft_sc,etot
1627 print *, 'Calling MINIM_DC'
1628 call minim_dc(etot,iretcode,nfun)
1630 call geom_to_var(nvar,varia)
1631 print *,'Calling MINIMIZE.'
1632 call minimize(etot,varia,iretcode,nfun)
1633 call var_to_geom(nvar,varia)
1635 if(me.eq.king.or..not.out1file)
1636 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1639 call chainbuild_cart
1644 kinetic_T=2.0d0/(dimen3*Rb)*EK
1645 if(me.eq.king.or..not.out1file)then
1655 call etotal(potEcomp)
1656 call enerprint(potEcomp)
1657 if (large) call enerprint(potEcomp)
1660 t_etotal=t_etotal+MPI_Wtime()-tt0
1662 t_etotal=t_etotal+tcpu()-tt0
1669 if (amax*d_time .gt. dvmax) then
1670 d_time=d_time*dvmax/amax
1671 if(me.eq.king.or..not.out1file) write (iout,*)
1672 & "Time step reduced to",d_time,
1673 & " because of too large initial acceleration."
1675 C if(me.eq.king.or..not.out1file)then
1676 C write(iout,*) "Potential energy and its components"
1677 C call enerprint(potEcomp)
1678 c write(iout,*) (potEcomp(i),i=0,n_ene)
1680 potE=potEcomp(0)-potEcomp(20)
1683 if (ntwe.ne.0) call statout(itime)
1684 if(me.eq.king.or..not.out1file)
1685 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1686 & " Kinetic energy",EK," Potential energy",potE,
1687 & " Total energy",totE," Maximum acceleration ",
1690 write (iout,*) "Initial coordinates"
1692 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1693 & (c(j,i+nres),j=1,3)
1695 write (iout,*) "Initial dC"
1697 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1698 & (dc(j,i+nres),j=1,3)
1700 write (iout,*) "Initial velocities"
1701 write (iout,"(13x,' backbone ',23x,' side chain')")
1703 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1704 & (d_t(j,i+nres),j=1,3)
1706 write (iout,*) "Initial accelerations"
1708 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1709 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1710 & (d_a(j,i+nres),j=1,3)
1716 d_t_old(j,i)=d_t(j,i)
1717 d_a_old(j,i)=d_a(j,i)
1719 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1728 call etotal_short(energia_short)
1729 if (large) call enerprint(potEcomp)
1732 t_eshort=t_eshort+MPI_Wtime()-tt0
1734 t_eshort=t_eshort+tcpu()-tt0
1739 if(.not.out1file .and. large) then
1740 write (iout,*) "energia_long",energia_long(0),
1741 & " energia_short",energia_short(0),
1742 & " total",energia_long(0)+energia_short(0)
1743 write (iout,*) "Initial fast-force accelerations"
1745 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1746 & (d_a(j,i+nres),j=1,3)
1749 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1752 d_a_short(j,i)=d_a(j,i)
1761 call etotal_long(energia_long)
1762 if (large) call enerprint(potEcomp)
1765 t_elong=t_elong+MPI_Wtime()-tt0
1767 t_elong=t_elong+tcpu()-tt0
1772 if(.not.out1file .and. large) then
1773 write (iout,*) "energia_long",energia_long(0)
1774 write (iout,*) "Initial slow-force accelerations"
1776 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1777 & (d_a(j,i+nres),j=1,3)
1781 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1783 t_enegrad=t_enegrad+tcpu()-tt0
1788 c-----------------------------------------------------------
1789 subroutine random_vel
1790 implicit real*8 (a-h,o-z)
1791 include 'DIMENSIONS'
1792 include 'COMMON.CONTROL'
1793 include 'COMMON.VAR'
1796 include 'COMMON.LANGEVIN'
1798 include 'COMMON.LANGEVIN.lang0'
1800 include 'COMMON.CHAIN'
1801 include 'COMMON.DERIV'
1802 include 'COMMON.GEO'
1803 include 'COMMON.LOCAL'
1804 include 'COMMON.INTERACT'
1805 include 'COMMON.IOUNITS'
1806 include 'COMMON.NAMES'
1807 include 'COMMON.TIME1'
1808 double precision xv,sigv,lowb,highb,vec_afm(3)
1809 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1810 c First generate velocities in the eigenspace of the G matrix
1811 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1818 sigv=dsqrt((Rb*t_bath)/geigen(i))
1821 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1823 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1824 c & " d_t_work_new",d_t_work_new(ii)
1827 C if (SELFGUIDE.gt.0) then
1830 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
1831 C distance=distance+vec_afm(j)**2
1833 C distance=dsqrt(distance)
1835 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
1836 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
1837 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
1838 C & d_t_work_new(j+(afmend-1)*3)
1849 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1852 c write (iout,*) "Ek from eigenvectors",Ek1
1854 c Transform velocities to UNRES coordinate space
1860 d_t_work(ind)=d_t_work(ind)
1861 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1863 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1867 c Transfer to the d_t vector
1869 d_t(j,0)=d_t_work(j)
1875 d_t(j,i)=d_t_work(ind)
1879 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1882 d_t(j,i+nres)=d_t_work(ind)
1887 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1888 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1893 c-----------------------------------------------------------
1894 subroutine sd_verlet_p_setup
1895 c Sets up the parameters of stochastic Verlet algorithm
1896 implicit real*8 (a-h,o-z)
1897 include 'DIMENSIONS'
1901 include 'COMMON.CONTROL'
1902 include 'COMMON.VAR'
1905 include 'COMMON.LANGEVIN'
1907 include 'COMMON.LANGEVIN.lang0'
1909 include 'COMMON.CHAIN'
1910 include 'COMMON.DERIV'
1911 include 'COMMON.GEO'
1912 include 'COMMON.LOCAL'
1913 include 'COMMON.INTERACT'
1914 include 'COMMON.IOUNITS'
1915 include 'COMMON.NAMES'
1916 include 'COMMON.TIME1'
1917 double precision emgdt(MAXRES6),
1918 & pterm,vterm,rho,rhoc,vsig,
1919 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1920 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1921 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1922 logical lprn /.false./
1923 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1924 double precision ktm
1931 c AL 8/17/04 Code adapted from tinker
1933 c Get the frictional and random terms for stochastic dynamics in the
1934 c eigenspace of mass-scaled UNRES friction matrix
1937 gdt = fricgam(i) * d_time
1939 c Stochastic dynamics reduces to simple MD for zero friction
1941 if (gdt .le. zero) then
1942 pfric_vec(i) = 1.0d0
1943 vfric_vec(i) = d_time
1944 afric_vec(i) = 0.5d0 * d_time * d_time
1945 prand_vec(i) = 0.0d0
1946 vrand_vec1(i) = 0.0d0
1947 vrand_vec2(i) = 0.0d0
1949 c Analytical expressions when friction coefficient is large
1952 if (gdt .ge. gdt_radius) then
1955 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1956 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1957 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1958 vterm = 1.0d0 - egdt**2
1959 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1961 c Use series expansions when friction coefficient is small
1972 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1973 & - gdt5/120.0d0 + gdt6/720.0d0
1974 & - gdt7/5040.0d0 + gdt8/40320.0d0
1975 & - gdt9/362880.0d0) / fricgam(i)**2
1976 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1977 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1978 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1979 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1980 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1981 & + 127.0d0*gdt9/90720.0d0
1982 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1983 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1984 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1985 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1986 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1987 & - 17.0d0*gdt2/1280.0d0
1988 & + 17.0d0*gdt3/6144.0d0
1989 & + 40967.0d0*gdt4/34406400.0d0
1990 & - 57203.0d0*gdt5/275251200.0d0
1991 & - 1429487.0d0*gdt6/13212057600.0d0)
1994 c Compute the scaling factors of random terms for the nonzero friction case
1996 ktm = 0.5d0*d_time/fricgam(i)
1997 psig = dsqrt(ktm*pterm) / fricgam(i)
1998 vsig = dsqrt(ktm*vterm)
1999 rhoc = dsqrt(1.0d0 - rho*rho)
2001 vrand_vec1(i) = vsig * rho
2002 vrand_vec2(i) = vsig * rhoc
2007 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2010 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2011 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2015 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2018 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2019 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2020 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2021 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2022 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2023 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2026 t_sdsetup=t_sdsetup+MPI_Wtime()
2028 t_sdsetup=t_sdsetup+tcpu()-tt0
2032 c-------------------------------------------------------------
2033 subroutine eigtransf1(n,ndim,ab,d,c)
2036 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2042 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2048 c-------------------------------------------------------------
2049 subroutine eigtransf(n,ndim,a,b,d,c)
2052 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2058 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2064 c-------------------------------------------------------------
2065 subroutine sd_verlet1
2066 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2067 implicit real*8 (a-h,o-z)
2068 include 'DIMENSIONS'
2069 include 'COMMON.CONTROL'
2070 include 'COMMON.VAR'
2073 include 'COMMON.LANGEVIN'
2075 include 'COMMON.LANGEVIN.lang0'
2077 include 'COMMON.CHAIN'
2078 include 'COMMON.DERIV'
2079 include 'COMMON.GEO'
2080 include 'COMMON.LOCAL'
2081 include 'COMMON.INTERACT'
2082 include 'COMMON.IOUNITS'
2083 include 'COMMON.NAMES'
2084 double precision stochforcvec(MAXRES6)
2085 common /stochcalc/ stochforcvec
2086 logical lprn /.false./
2088 c write (iout,*) "dc_old"
2090 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2091 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2094 dc_work(j)=dc_old(j,0)
2095 d_t_work(j)=d_t_old(j,0)
2096 d_a_work(j)=d_a_old(j,0)
2101 dc_work(ind+j)=dc_old(j,i)
2102 d_t_work(ind+j)=d_t_old(j,i)
2103 d_a_work(ind+j)=d_a_old(j,i)
2108 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2110 dc_work(ind+j)=dc_old(j,i+nres)
2111 d_t_work(ind+j)=d_t_old(j,i+nres)
2112 d_a_work(ind+j)=d_a_old(j,i+nres)
2120 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2124 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2125 & vfric_mat(i,j),afric_mat(i,j),
2126 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2134 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2135 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2136 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2137 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2139 d_t_work_new(i)=ddt1+0.5d0*ddt2
2140 d_t_work(i)=ddt1+ddt2
2145 d_t(j,0)=d_t_work(j)
2150 dc(j,i)=dc_work(ind+j)
2151 d_t(j,i)=d_t_work(ind+j)
2156 if (itype(i).ne.10) then
2159 dc(j,inres)=dc_work(ind+j)
2160 d_t(j,inres)=d_t_work(ind+j)
2167 c--------------------------------------------------------------------------
2168 subroutine sd_verlet2
2169 c Calculating the adjusted velocities for accelerations
2170 implicit real*8 (a-h,o-z)
2171 include 'DIMENSIONS'
2172 include 'COMMON.CONTROL'
2173 include 'COMMON.VAR'
2176 include 'COMMON.LANGEVIN'
2178 include 'COMMON.LANGEVIN.lang0'
2180 include 'COMMON.CHAIN'
2181 include 'COMMON.DERIV'
2182 include 'COMMON.GEO'
2183 include 'COMMON.LOCAL'
2184 include 'COMMON.INTERACT'
2185 include 'COMMON.IOUNITS'
2186 include 'COMMON.NAMES'
2187 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2188 common /stochcalc/ stochforcvec
2190 c Compute the stochastic forces which contribute to velocity change
2192 call stochastic_force(stochforcvecV)
2199 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2200 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2201 & vrand_mat2(i,j)*stochforcvecV(j)
2203 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2207 d_t(j,0)=d_t_work(j)
2212 d_t(j,i)=d_t_work(ind+j)
2217 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2220 d_t(j,inres)=d_t_work(ind+j)
2227 c-----------------------------------------------------------
2228 subroutine sd_verlet_ciccotti_setup
2229 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2231 implicit real*8 (a-h,o-z)
2232 include 'DIMENSIONS'
2236 include 'COMMON.CONTROL'
2237 include 'COMMON.VAR'
2240 include 'COMMON.LANGEVIN'
2242 include 'COMMON.LANGEVIN.lang0'
2244 include 'COMMON.CHAIN'
2245 include 'COMMON.DERIV'
2246 include 'COMMON.GEO'
2247 include 'COMMON.LOCAL'
2248 include 'COMMON.INTERACT'
2249 include 'COMMON.IOUNITS'
2250 include 'COMMON.NAMES'
2251 include 'COMMON.TIME1'
2252 double precision emgdt(MAXRES6),
2253 & pterm,vterm,rho,rhoc,vsig,
2254 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2255 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2256 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2257 logical lprn /.false./
2258 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2259 double precision ktm
2266 c AL 8/17/04 Code adapted from tinker
2268 c Get the frictional and random terms for stochastic dynamics in the
2269 c eigenspace of mass-scaled UNRES friction matrix
2272 write (iout,*) "i",i," fricgam",fricgam(i)
2273 gdt = fricgam(i) * d_time
2275 c Stochastic dynamics reduces to simple MD for zero friction
2277 if (gdt .le. zero) then
2278 pfric_vec(i) = 1.0d0
2279 vfric_vec(i) = d_time
2280 afric_vec(i) = 0.5d0*d_time*d_time
2281 prand_vec(i) = afric_vec(i)
2282 vrand_vec2(i) = vfric_vec(i)
2284 c Analytical expressions when friction coefficient is large
2289 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2290 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2291 prand_vec(i) = afric_vec(i)
2292 vrand_vec2(i) = vfric_vec(i)
2294 c Compute the scaling factors of random terms for the nonzero friction case
2296 c ktm = 0.5d0*d_time/fricgam(i)
2297 c psig = dsqrt(ktm*pterm) / fricgam(i)
2298 c vsig = dsqrt(ktm*vterm)
2299 c prand_vec(i) = psig*afric_vec(i)
2300 c vrand_vec2(i) = vsig*vfric_vec(i)
2305 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2308 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2309 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2313 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2315 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2316 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2317 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2318 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2319 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2321 t_sdsetup=t_sdsetup+MPI_Wtime()
2323 t_sdsetup=t_sdsetup+tcpu()-tt0
2327 c-------------------------------------------------------------
2328 subroutine sd_verlet1_ciccotti
2329 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2330 implicit real*8 (a-h,o-z)
2331 include 'DIMENSIONS'
2335 include 'COMMON.CONTROL'
2336 include 'COMMON.VAR'
2339 include 'COMMON.LANGEVIN'
2341 include 'COMMON.LANGEVIN.lang0'
2343 include 'COMMON.CHAIN'
2344 include 'COMMON.DERIV'
2345 include 'COMMON.GEO'
2346 include 'COMMON.LOCAL'
2347 include 'COMMON.INTERACT'
2348 include 'COMMON.IOUNITS'
2349 include 'COMMON.NAMES'
2350 double precision stochforcvec(MAXRES6)
2351 common /stochcalc/ stochforcvec
2352 logical lprn /.false./
2354 c write (iout,*) "dc_old"
2356 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2357 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2360 dc_work(j)=dc_old(j,0)
2361 d_t_work(j)=d_t_old(j,0)
2362 d_a_work(j)=d_a_old(j,0)
2367 dc_work(ind+j)=dc_old(j,i)
2368 d_t_work(ind+j)=d_t_old(j,i)
2369 d_a_work(ind+j)=d_a_old(j,i)
2374 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2376 dc_work(ind+j)=dc_old(j,i+nres)
2377 d_t_work(ind+j)=d_t_old(j,i+nres)
2378 d_a_work(ind+j)=d_a_old(j,i+nres)
2387 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2391 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2392 & vfric_mat(i,j),afric_mat(i,j),
2393 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2401 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2402 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2403 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2404 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2406 d_t_work_new(i)=ddt1+0.5d0*ddt2
2407 d_t_work(i)=ddt1+ddt2
2412 d_t(j,0)=d_t_work(j)
2417 dc(j,i)=dc_work(ind+j)
2418 d_t(j,i)=d_t_work(ind+j)
2423 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2426 dc(j,inres)=dc_work(ind+j)
2427 d_t(j,inres)=d_t_work(ind+j)
2434 c--------------------------------------------------------------------------
2435 subroutine sd_verlet2_ciccotti
2436 c Calculating the adjusted velocities for accelerations
2437 implicit real*8 (a-h,o-z)
2438 include 'DIMENSIONS'
2439 include 'COMMON.CONTROL'
2440 include 'COMMON.VAR'
2443 include 'COMMON.LANGEVIN'
2445 include 'COMMON.LANGEVIN.lang0'
2447 include 'COMMON.CHAIN'
2448 include 'COMMON.DERIV'
2449 include 'COMMON.GEO'
2450 include 'COMMON.LOCAL'
2451 include 'COMMON.INTERACT'
2452 include 'COMMON.IOUNITS'
2453 include 'COMMON.NAMES'
2454 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2455 common /stochcalc/ stochforcvec
2457 c Compute the stochastic forces which contribute to velocity change
2459 call stochastic_force(stochforcvecV)
2466 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2467 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2468 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2470 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2474 d_t(j,0)=d_t_work(j)
2479 d_t(j,i)=d_t_work(ind+j)
2484 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2487 d_t(j,inres)=d_t_work(ind+j)