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 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
523 if (d_time.ne.d_time0) then
526 if (lang.eq.2 .or. lang.eq.3) then
527 if (large) write (iout,'(a)')
528 & "Restore original matrices for stochastic step"
529 c write (iout,*) "Calling sd_verlet_setup: 2"
530 c Restore the matrices of tinker integrator if the time step has been restored
533 pfric_mat(i,j)=pfric0_mat(i,j,0)
534 afric_mat(i,j)=afric0_mat(i,j,0)
535 vfric_mat(i,j)=vfric0_mat(i,j,0)
536 prand_mat(i,j)=prand0_mat(i,j,0)
537 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
538 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
547 c Calculate the kinetic and the total energy and the kinetic temperature
552 c write (iout,*) "step",itime," EK",EK," EK1",EK1
554 c Couple the system to Berendsen bath if needed
555 if (tbf .and. lang.eq.0) then
558 kinetic_T=2.0d0/(dimen3*Rb)*EK
559 c Backup the coordinates, velocities, and accelerations
563 d_t_old(j,i)=d_t(j,i)
564 d_a_old(j,i)=d_a(j,i)
568 if (mod(itime,ntwe).eq.0 .and. large) then
569 write (iout,*) "Velocities, step 2"
571 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
572 & (d_t(j,i+nres),j=1,3)
578 c-------------------------------------------------------------------------------
579 subroutine RESPA_step(itime)
580 c-------------------------------------------------------------------------------
581 c Perform a single RESPA step.
582 c-------------------------------------------------------------------------------
583 implicit real*8 (a-h,o-z)
587 integer IERROR,ERRCODE
589 include 'COMMON.SETUP'
590 include 'COMMON.CONTROL'
594 include 'COMMON.LANGEVIN'
596 include 'COMMON.LANGEVIN.lang0'
598 include 'COMMON.CHAIN'
599 include 'COMMON.DERIV'
601 include 'COMMON.LOCAL'
602 include 'COMMON.INTERACT'
603 include 'COMMON.IOUNITS'
604 include 'COMMON.NAMES'
605 include 'COMMON.TIME1'
606 double precision energia_short(0:n_ene),
607 & energia_long(0:n_ene)
608 double precision cm(3),L(3),vcm(3),incr(3)
609 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
610 & d_a_old0(3,0:maxres2)
611 logical PRINT_AMTS_MSG /.false./
612 integer ilen,count,rstcount
615 integer maxcount_scale /10/
617 double precision stochforcvec(MAXRES6)
618 common /stochcalc/ stochforcvec
621 common /cipiszcze/ itt
624 if (large.and. mod(itime,ntwe).eq.0) then
625 write (iout,*) "***************** RESPA itime",itime
626 write (iout,*) "Cartesian and internal coordinates: step 0"
628 call pdbout(0.0d0,"cipiszcze",iout)
630 write (iout,*) "Accelerations from long-range forces"
632 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
633 & (d_a(j,i+nres),j=1,3)
635 write (iout,*) "Velocities, step 0"
637 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
638 & (d_t(j,i+nres),j=1,3)
643 c Perform the initial RESPA step (increment velocities)
644 c write (iout,*) "*********************** RESPA ini"
647 if (mod(itime,ntwe).eq.0 .and. large) then
648 write (iout,*) "Velocities, end"
650 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
651 & (d_t(j,i+nres),j=1,3)
655 c Compute the short-range forces
661 C 7/2/2009 commented out
663 c call etotal_short(energia_short)
666 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
669 d_a(j,i)=d_a_short(j,i)
673 if (large.and. mod(itime,ntwe).eq.0) then
674 write (iout,*) "energia_short",energia_short(0)
675 write (iout,*) "Accelerations from short-range forces"
677 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
678 & (d_a(j,i+nres),j=1,3)
683 t_enegrad=t_enegrad+MPI_Wtime()-tt0
685 t_enegrad=t_enegrad+tcpu()-tt0
690 d_t_old(j,i)=d_t(j,i)
691 d_a_old(j,i)=d_a(j,i)
694 c 6/30/08 A-MTS: attempt at increasing the split number
697 dc_old0(j,i)=dc_old(j,i)
698 d_t_old0(j,i)=d_t_old(j,i)
699 d_a_old0(j,i)=d_a_old(j,i)
702 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
703 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
710 c write (iout,*) "itime",itime," ntime_split",ntime_split
711 c Split the time step
712 d_time=d_time0/ntime_split
713 c Perform the short-range RESPA steps (velocity Verlet increments of
714 c positions and velocities using short-range forces)
715 c write (iout,*) "*********************** RESPA split"
716 do itsplit=1,ntime_split
719 else if (lang.eq.2 .or. lang.eq.3) then
721 call stochastic_force(stochforcvec)
724 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
726 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
731 c First step of the velocity Verlet algorithm
736 else if (lang.eq.3) then
738 call sd_verlet1_ciccotti
740 else if (lang.eq.1) then
745 c Build the chain from the newly calculated coordinates
747 if (rattle) call rattle1
749 if (large.and. mod(itime,ntwe).eq.0) then
750 write (iout,*) "***** ITSPLIT",itsplit
751 write (iout,*) "Cartesian and internal coordinates: step 1"
752 call pdbout(0.0d0,"cipiszcze",iout)
755 write (iout,*) "Velocities, step 1"
757 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
758 & (d_t(j,i+nres),j=1,3)
767 c Calculate energy and forces
769 call etotal_short(energia_short)
770 if (large.and. mod(itime,ntwe).eq.0)
771 & call enerprint(energia_short)
774 t_eshort=t_eshort+MPI_Wtime()-tt0
776 t_eshort=t_eshort+tcpu()-tt0
780 c Get the new accelerations
782 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
785 d_a_short(j,i)=d_a(j,i)
789 if (large.and. mod(itime,ntwe).eq.0) then
790 write (iout,*)"energia_short",energia_short(0)
791 write (iout,*) "Accelerations from short-range forces"
793 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
794 & (d_a(j,i+nres),j=1,3)
799 c Determine maximum acceleration and scale down the timestep if needed
801 amax=amax/ntime_split**2
802 call predict_edrift(epdrift)
803 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
804 & write (iout,*) "amax",amax," damax",damax,
805 & " epdrift",epdrift," epdriftmax",epdriftmax
806 c Exit loop and try with increased split number if the change of
807 c acceleration is too big
808 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
809 if (ntime_split.lt.maxtime_split) then
811 ntime_split=ntime_split*2
814 dc_old(j,i)=dc_old0(j,i)
815 d_t_old(j,i)=d_t_old0(j,i)
816 d_a_old(j,i)=d_a_old0(j,i)
819 if (PRINT_AMTS_MSG) then
820 write (iout,*) "acceleration/energy drift too large",amax,
821 & epdrift," split increased to ",ntime_split," itime",itime,
827 & "Uh-hu. Bumpy landscape. Maximum splitting number",
829 & " already reached!!! Trying to carry on!"
833 t_enegrad=t_enegrad+MPI_Wtime()-tt0
835 t_enegrad=t_enegrad+tcpu()-tt0
837 c Second step of the velocity Verlet algorithm
842 else if (lang.eq.3) then
844 call sd_verlet2_ciccotti
846 else if (lang.eq.1) then
851 if (rattle) call rattle2
852 c Backup the coordinates, velocities, and accelerations
856 d_t_old(j,i)=d_t(j,i)
857 d_a_old(j,i)=d_a(j,i)
864 c Restore the time step
866 c Compute long-range forces
873 call etotal_long(energia_long)
874 if (large.and. mod(itime,ntwe).eq.0)
875 & call enerprint(energia_long)
878 t_elong=t_elong+MPI_Wtime()-tt0
880 t_elong=t_elong+tcpu()-tt0
886 t_enegrad=t_enegrad+MPI_Wtime()-tt0
888 t_enegrad=t_enegrad+tcpu()-tt0
890 c Compute accelerations from long-range forces
892 if (large.and. mod(itime,ntwe).eq.0) then
893 write (iout,*) "energia_long",energia_long(0)
894 write (iout,*) "Cartesian and internal coordinates: step 2"
896 call pdbout(0.0d0,"cipiszcze",iout)
898 write (iout,*) "Accelerations from long-range forces"
900 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
901 & (d_a(j,i+nres),j=1,3)
903 write (iout,*) "Velocities, step 2"
905 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
906 & (d_t(j,i+nres),j=1,3)
910 c Compute the final RESPA step (increment velocities)
911 c write (iout,*) "*********************** RESPA fin"
913 c Compute the complete potential energy
915 potEcomp(i)=energia_short(i)+energia_long(i)
917 potE=potEcomp(0)-potEcomp(20)
918 c potE=energia_short(0)+energia_long(0)
920 c Calculate the kinetic and the total energy and the kinetic temperature
923 c Couple the system to Berendsen bath if needed
924 if (tbf .and. lang.eq.0) then
927 kinetic_T=2.0d0/(dimen3*Rb)*EK
928 c Backup the coordinates, velocities, and accelerations
930 if (mod(itime,ntwe).eq.0 .and. large) then
931 write (iout,*) "Velocities, end"
933 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
934 & (d_t(j,i+nres),j=1,3)
940 c---------------------------------------------------------------------
942 c First and last RESPA step (incrementing velocities using long-range
944 implicit real*8 (a-h,o-z)
946 include 'COMMON.CONTROL'
949 include 'COMMON.CHAIN'
950 include 'COMMON.DERIV'
952 include 'COMMON.LOCAL'
953 include 'COMMON.INTERACT'
954 include 'COMMON.IOUNITS'
955 include 'COMMON.NAMES'
957 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
961 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
965 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
968 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
974 c-----------------------------------------------------------------
976 c Applying velocity Verlet algorithm - step 1 to coordinates
977 implicit real*8 (a-h,o-z)
979 include 'COMMON.CONTROL'
982 include 'COMMON.CHAIN'
983 include 'COMMON.DERIV'
985 include 'COMMON.LOCAL'
986 include 'COMMON.INTERACT'
987 include 'COMMON.IOUNITS'
988 include 'COMMON.NAMES'
989 double precision adt,adt2
992 write (iout,*) "VELVERLET1 START: DC"
994 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
995 & (dc(j,i+nres),j=1,3)
999 adt=d_a_old(j,0)*d_time
1001 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1002 d_t_new(j,0)=d_t_old(j,0)+adt2
1003 d_t(j,0)=d_t_old(j,0)+adt
1009 adt=d_a_old(j,i)*d_time
1011 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1012 d_t_new(j,i)=d_t_old(j,i)+adt2
1013 d_t(j,i)=d_t_old(j,i)+adt
1018 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1021 adt=d_a_old(j,inres)*d_time
1023 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1024 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1025 d_t(j,inres)=d_t_old(j,inres)+adt
1030 write (iout,*) "VELVERLET1 END: DC"
1032 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1033 & (dc(j,i+nres),j=1,3)
1038 c---------------------------------------------------------------------
1040 c Step 2 of the velocity Verlet algorithm: update velocities
1041 implicit real*8 (a-h,o-z)
1042 include 'DIMENSIONS'
1043 include 'COMMON.CONTROL'
1044 include 'COMMON.VAR'
1046 include 'COMMON.CHAIN'
1047 include 'COMMON.DERIV'
1048 include 'COMMON.GEO'
1049 include 'COMMON.LOCAL'
1050 include 'COMMON.INTERACT'
1051 include 'COMMON.IOUNITS'
1052 include 'COMMON.NAMES'
1054 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1058 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1062 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1065 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1071 c-----------------------------------------------------------------
1072 subroutine sddir_precalc
1073 c Applying velocity Verlet algorithm - step 1 to coordinates
1074 implicit real*8 (a-h,o-z)
1075 include 'DIMENSIONS'
1079 include 'COMMON.CONTROL'
1080 include 'COMMON.VAR'
1083 include 'COMMON.LANGEVIN'
1085 include 'COMMON.LANGEVIN.lang0'
1087 include 'COMMON.CHAIN'
1088 include 'COMMON.DERIV'
1089 include 'COMMON.GEO'
1090 include 'COMMON.LOCAL'
1091 include 'COMMON.INTERACT'
1092 include 'COMMON.IOUNITS'
1093 include 'COMMON.NAMES'
1094 include 'COMMON.TIME1'
1095 double precision stochforcvec(MAXRES6)
1096 common /stochcalc/ stochforcvec
1098 c Compute friction and stochastic forces
1107 time_fric=time_fric+MPI_Wtime()-time00
1110 time_fric=time_fric+tcpu()-time00
1113 call stochastic_force(stochforcvec)
1115 time_stoch=time_stoch+MPI_Wtime()-time00
1117 time_stoch=time_stoch+tcpu()-time00
1120 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1121 c forces (d_as_work)
1123 call ginv_mult(fric_work, d_af_work)
1124 call ginv_mult(stochforcvec, d_as_work)
1127 c---------------------------------------------------------------------
1128 subroutine sddir_verlet1
1129 c Applying velocity Verlet algorithm - step 1 to velocities
1130 implicit real*8 (a-h,o-z)
1131 include 'DIMENSIONS'
1132 include 'COMMON.CONTROL'
1133 include 'COMMON.VAR'
1136 include 'COMMON.LANGEVIN'
1138 include 'COMMON.LANGEVIN.lang0'
1140 include 'COMMON.CHAIN'
1141 include 'COMMON.DERIV'
1142 include 'COMMON.GEO'
1143 include 'COMMON.LOCAL'
1144 include 'COMMON.INTERACT'
1145 include 'COMMON.IOUNITS'
1146 include 'COMMON.NAMES'
1147 c Revised 3/31/05 AL: correlation between random contributions to
1148 c position and velocity increments included.
1149 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1150 double precision adt,adt2
1152 c Add the contribution from BOTH friction and stochastic force to the
1153 c coordinates, but ONLY the contribution from the friction forces to velocities
1156 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1157 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1158 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1159 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1160 d_t(j,0)=d_t_old(j,0)+adt
1165 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1166 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1167 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1168 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1169 d_t(j,i)=d_t_old(j,i)+adt
1174 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1177 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1178 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1179 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1180 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1181 d_t(j,inres)=d_t_old(j,inres)+adt
1188 c---------------------------------------------------------------------
1189 subroutine sddir_verlet2
1190 c Calculating the adjusted velocities for accelerations
1191 implicit real*8 (a-h,o-z)
1192 include 'DIMENSIONS'
1193 include 'COMMON.CONTROL'
1194 include 'COMMON.VAR'
1197 include 'COMMON.LANGEVIN'
1199 include 'COMMON.LANGEVIN.lang0'
1201 include 'COMMON.CHAIN'
1202 include 'COMMON.DERIV'
1203 include 'COMMON.GEO'
1204 include 'COMMON.LOCAL'
1205 include 'COMMON.INTERACT'
1206 include 'COMMON.IOUNITS'
1207 include 'COMMON.NAMES'
1208 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1209 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1210 c Revised 3/31/05 AL: correlation between random contributions to
1211 c position and velocity increments included.
1212 c The correlation coefficients are calculated at low-friction limit.
1213 c Also, friction forces are now not calculated with new velocities.
1215 c call friction_force
1216 call stochastic_force(stochforcvec)
1218 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1219 c forces (d_as_work)
1221 call ginv_mult(stochforcvec, d_as_work1)
1227 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1228 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1233 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1234 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1239 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1242 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1243 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1244 & +cos60*d_as_work1(ind+j))*d_time
1251 c---------------------------------------------------------------------
1252 subroutine max_accel
1254 c Find the maximum difference in the accelerations of the the sites
1255 c at the beginning and the end of the time step.
1257 implicit real*8 (a-h,o-z)
1258 include 'DIMENSIONS'
1259 include 'COMMON.CONTROL'
1260 include 'COMMON.VAR'
1262 include 'COMMON.CHAIN'
1263 include 'COMMON.DERIV'
1264 include 'COMMON.GEO'
1265 include 'COMMON.LOCAL'
1266 include 'COMMON.INTERACT'
1267 include 'COMMON.IOUNITS'
1268 double precision aux(3),accel(3),accel_old(3),dacc
1270 c aux(j)=d_a(j,0)-d_a_old(j,0)
1271 accel_old(j)=d_a_old(j,0)
1278 c 7/3/08 changed to asymmetric difference
1280 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1281 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1282 accel(j)=accel(j)+0.5d0*d_a(j,i)
1283 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1284 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1285 dacc=dabs(accel(j)-accel_old(j))
1286 c write (iout,*) i,dacc
1287 if (dacc.gt.amax) amax=dacc
1295 accel_old(j)=d_a_old(j,0)
1300 accel_old(j)=accel_old(j)+d_a_old(j,1)
1301 accel(j)=accel(j)+d_a(j,1)
1305 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1307 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1308 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1309 accel(j)=accel(j)+d_a(j,i+nres)
1313 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1314 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1315 dacc=dabs(accel(j)-accel_old(j))
1316 c write (iout,*) "side-chain",i,dacc
1317 if (dacc.gt.amax) amax=dacc
1321 accel_old(j)=accel_old(j)+d_a_old(j,i)
1322 accel(j)=accel(j)+d_a(j,i)
1323 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1328 c---------------------------------------------------------------------
1329 subroutine predict_edrift(epdrift)
1331 c Predict the drift of the potential energy
1333 implicit real*8 (a-h,o-z)
1334 include 'DIMENSIONS'
1335 include 'COMMON.CONTROL'
1336 include 'COMMON.VAR'
1338 include 'COMMON.CHAIN'
1339 include 'COMMON.DERIV'
1340 include 'COMMON.GEO'
1341 include 'COMMON.LOCAL'
1342 include 'COMMON.INTERACT'
1343 include 'COMMON.IOUNITS'
1344 include 'COMMON.MUCA'
1345 double precision epdrift,epdriftij
1346 c Drift of the potential energy
1352 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1353 if (lmuca) epdriftij=epdriftij*factor
1354 c write (iout,*) "back",i,j,epdriftij
1355 if (epdriftij.gt.epdrift) epdrift=epdriftij
1359 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1362 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1363 if (lmuca) epdriftij=epdriftij*factor
1364 c write (iout,*) "side",i,j,epdriftij
1365 if (epdriftij.gt.epdrift) epdrift=epdriftij
1369 epdrift=0.5d0*epdrift*d_time*d_time
1370 c write (iout,*) "epdrift",epdrift
1373 c-----------------------------------------------------------------------
1374 subroutine verlet_bath
1376 c Coupling to the thermostat by using the Berendsen algorithm
1378 implicit real*8 (a-h,o-z)
1379 include 'DIMENSIONS'
1380 include 'COMMON.CONTROL'
1381 include 'COMMON.VAR'
1383 include 'COMMON.CHAIN'
1384 include 'COMMON.DERIV'
1385 include 'COMMON.GEO'
1386 include 'COMMON.LOCAL'
1387 include 'COMMON.INTERACT'
1388 include 'COMMON.IOUNITS'
1389 include 'COMMON.NAMES'
1390 double precision T_half,fact
1392 T_half=2.0d0/(dimen3*Rb)*EK
1393 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1394 c write(iout,*) "T_half", T_half
1395 c write(iout,*) "EK", EK
1396 c write(iout,*) "fact", fact
1398 d_t(j,0)=fact*d_t(j,0)
1402 d_t(j,i)=fact*d_t(j,i)
1406 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1409 d_t(j,inres)=fact*d_t(j,inres)
1415 c---------------------------------------------------------
1417 c Set up the initial conditions of a MD simulation
1418 implicit real*8 (a-h,o-z)
1419 include 'DIMENSIONS'
1423 integer IERROR,ERRCODE
1425 include 'COMMON.SETUP'
1426 include 'COMMON.CONTROL'
1427 include 'COMMON.VAR'
1430 include 'COMMON.LANGEVIN'
1432 include 'COMMON.LANGEVIN.lang0'
1434 include 'COMMON.CHAIN'
1435 include 'COMMON.DERIV'
1436 include 'COMMON.GEO'
1437 include 'COMMON.LOCAL'
1438 include 'COMMON.INTERACT'
1439 include 'COMMON.IOUNITS'
1440 include 'COMMON.NAMES'
1441 include 'COMMON.REMD'
1442 real*8 energia_long(0:n_ene),
1443 & energia_short(0:n_ene),vcm(3),incr(3)
1444 double precision cm(3),L(3),xv,sigv,lowb,highb
1445 double precision varia(maxvar)
1453 c write(iout,*) "d_time", d_time
1454 c Compute the standard deviations of stochastic forces for Langevin dynamics
1455 c if the friction coefficients do not depend on surface area
1456 if (lang.gt.0 .and. .not.surfarea) then
1458 stdforcp(i)=stdfp*dsqrt(gamp)
1461 stdforcsc(i)=stdfsc(iabs(itype(i)))
1462 & *dsqrt(gamsc(iabs(itype(i))))
1465 c Open the pdb file for snapshotshots
1468 if (ilen(tmpdir).gt.0)
1469 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1470 & liczba(:ilen(liczba))//".pdb")
1472 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1476 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1477 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1478 & liczba(:ilen(liczba))//".x")
1479 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1482 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1483 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1484 & liczba(:ilen(liczba))//".cx")
1485 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1491 if (ilen(tmpdir).gt.0)
1492 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1493 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1495 if (ilen(tmpdir).gt.0)
1496 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1497 cartname=prefix(:ilen(prefix))//"_MD.cx"
1501 write (qstr,'(256(1h ))')
1504 iq = qinfrag(i,iset)*10
1505 iw = wfrag(i,iset)/100
1507 if(me.eq.king.or..not.out1file)
1508 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1509 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1514 iq = qinpair(i,iset)*10
1515 iw = wpair(i,iset)/100
1517 if(me.eq.king.or..not.out1file)
1518 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1519 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1523 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1525 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1527 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1529 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1533 if (restart1file) then
1535 & inquire(file=mremd_rst_name,exist=file_exist)
1537 write (*,*) me," Before broadcast: file_exist",file_exist
1538 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1540 write (*,*) me," After broadcast: file_exist",file_exist
1541 c inquire(file=mremd_rst_name,exist=file_exist)
1543 if(me.eq.king.or..not.out1file)
1544 & write(iout,*) "Initial state read by master and distributed"
1546 if (ilen(tmpdir).gt.0)
1547 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1548 & //liczba(:ilen(liczba))//'.rst')
1549 inquire(file=rest2name,exist=file_exist)
1552 if(.not.restart1file) then
1553 if(me.eq.king.or..not.out1file)
1554 & write(iout,*) "Initial state will be read from file ",
1555 & rest2name(:ilen(rest2name))
1558 call rescale_weights(t_bath)
1560 if(me.eq.king.or..not.out1file)then
1561 if (restart1file) then
1562 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1565 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1568 write(iout,*) "Initial velocities randomly generated"
1574 c Generate initial velocities
1575 if(me.eq.king.or..not.out1file)
1576 & write(iout,*) "Initial velocities randomly generated"
1580 c rest2name = prefix(:ilen(prefix))//'.rst'
1581 if(me.eq.king.or..not.out1file)then
1582 write (iout,*) "Initial velocities"
1584 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1585 & (d_t(j,i+nres),j=1,3)
1587 c Zeroing the total angular momentum of the system
1588 write(iout,*) "Calling the zero-angular
1589 & momentum subroutine"
1592 c Getting the potential energy and forces and velocities and accelerations
1594 c write (iout,*) "velocity of the center of the mass:"
1595 c write (iout,*) (vcm(j),j=1,3)
1597 d_t(j,0)=d_t(j,0)-vcm(j)
1599 c Removing the velocity of the center of mass
1601 if(me.eq.king.or..not.out1file)then
1602 write (iout,*) "vcm right after adjustment:"
1603 write (iout,*) (vcm(j),j=1,3)
1607 if(iranconf.ne.0) then
1609 print *, 'Calling OVERLAP_SC'
1610 call overlap_sc(fail)
1614 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1615 print *,'SC_move',nft_sc,etot
1616 if(me.eq.king.or..not.out1file)
1617 & write(iout,*) 'SC_move',nft_sc,etot
1621 print *, 'Calling MINIM_DC'
1622 call minim_dc(etot,iretcode,nfun)
1624 call geom_to_var(nvar,varia)
1625 print *,'Calling MINIMIZE.'
1626 call minimize(etot,varia,iretcode,nfun)
1627 call var_to_geom(nvar,varia)
1629 if(me.eq.king.or..not.out1file)
1630 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1633 call chainbuild_cart
1638 kinetic_T=2.0d0/(dimen3*Rb)*EK
1639 if(me.eq.king.or..not.out1file)then
1649 call etotal(potEcomp)
1650 if (large) call enerprint(potEcomp)
1653 t_etotal=t_etotal+MPI_Wtime()-tt0
1655 t_etotal=t_etotal+tcpu()-tt0
1662 if (amax*d_time .gt. dvmax) then
1663 d_time=d_time*dvmax/amax
1664 if(me.eq.king.or..not.out1file) write (iout,*)
1665 & "Time step reduced to",d_time,
1666 & " because of too large initial acceleration."
1668 if(me.eq.king.or..not.out1file)then
1669 write(iout,*) "Potential energy and its components"
1670 call enerprint(potEcomp)
1671 c write(iout,*) (potEcomp(i),i=0,n_ene)
1673 potE=potEcomp(0)-potEcomp(20)
1676 if (ntwe.ne.0) call statout(itime)
1677 if(me.eq.king.or..not.out1file)
1678 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1679 & " Kinetic energy",EK," Potential energy",potE,
1680 & " Total energy",totE," Maximum acceleration ",
1683 write (iout,*) "Initial coordinates"
1685 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1686 & (c(j,i+nres),j=1,3)
1688 write (iout,*) "Initial dC"
1690 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1691 & (dc(j,i+nres),j=1,3)
1693 write (iout,*) "Initial velocities"
1694 write (iout,"(13x,' backbone ',23x,' side chain')")
1696 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1697 & (d_t(j,i+nres),j=1,3)
1699 write (iout,*) "Initial accelerations"
1701 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1702 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1703 & (d_a(j,i+nres),j=1,3)
1709 d_t_old(j,i)=d_t(j,i)
1710 d_a_old(j,i)=d_a(j,i)
1712 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1721 call etotal_short(energia_short)
1722 if (large) call enerprint(potEcomp)
1725 t_eshort=t_eshort+MPI_Wtime()-tt0
1727 t_eshort=t_eshort+tcpu()-tt0
1732 if(.not.out1file .and. large) then
1733 write (iout,*) "energia_long",energia_long(0),
1734 & " energia_short",energia_short(0),
1735 & " total",energia_long(0)+energia_short(0)
1736 write (iout,*) "Initial fast-force accelerations"
1738 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1739 & (d_a(j,i+nres),j=1,3)
1742 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1745 d_a_short(j,i)=d_a(j,i)
1754 call etotal_long(energia_long)
1755 if (large) call enerprint(potEcomp)
1758 t_elong=t_elong+MPI_Wtime()-tt0
1760 t_elong=t_elong+tcpu()-tt0
1765 if(.not.out1file .and. large) then
1766 write (iout,*) "energia_long",energia_long(0)
1767 write (iout,*) "Initial slow-force accelerations"
1769 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1770 & (d_a(j,i+nres),j=1,3)
1774 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1776 t_enegrad=t_enegrad+tcpu()-tt0
1781 c-----------------------------------------------------------
1782 subroutine random_vel
1783 implicit real*8 (a-h,o-z)
1784 include 'DIMENSIONS'
1785 include 'COMMON.CONTROL'
1786 include 'COMMON.VAR'
1789 include 'COMMON.LANGEVIN'
1791 include 'COMMON.LANGEVIN.lang0'
1793 include 'COMMON.CHAIN'
1794 include 'COMMON.DERIV'
1795 include 'COMMON.GEO'
1796 include 'COMMON.LOCAL'
1797 include 'COMMON.INTERACT'
1798 include 'COMMON.IOUNITS'
1799 include 'COMMON.NAMES'
1800 include 'COMMON.TIME1'
1801 double precision xv,sigv,lowb,highb
1802 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1803 c First generate velocities in the eigenspace of the G matrix
1804 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1811 sigv=dsqrt((Rb*t_bath)/geigen(i))
1814 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1815 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1816 c & " d_t_work_new",d_t_work_new(ii)
1825 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1828 c write (iout,*) "Ek from eigenvectors",Ek1
1830 c Transform velocities to UNRES coordinate space
1836 d_t_work(ind)=d_t_work(ind)
1837 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1839 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1843 c Transfer to the d_t vector
1845 d_t(j,0)=d_t_work(j)
1851 d_t(j,i)=d_t_work(ind)
1855 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1858 d_t(j,i+nres)=d_t_work(ind)
1863 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1864 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1869 c-----------------------------------------------------------
1870 subroutine sd_verlet_p_setup
1871 c Sets up the parameters of stochastic Verlet algorithm
1872 implicit real*8 (a-h,o-z)
1873 include 'DIMENSIONS'
1877 include 'COMMON.CONTROL'
1878 include 'COMMON.VAR'
1881 include 'COMMON.LANGEVIN'
1883 include 'COMMON.LANGEVIN.lang0'
1885 include 'COMMON.CHAIN'
1886 include 'COMMON.DERIV'
1887 include 'COMMON.GEO'
1888 include 'COMMON.LOCAL'
1889 include 'COMMON.INTERACT'
1890 include 'COMMON.IOUNITS'
1891 include 'COMMON.NAMES'
1892 include 'COMMON.TIME1'
1893 double precision emgdt(MAXRES6),
1894 & pterm,vterm,rho,rhoc,vsig,
1895 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1896 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1897 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1898 logical lprn /.false./
1899 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1900 double precision ktm
1907 c AL 8/17/04 Code adapted from tinker
1909 c Get the frictional and random terms for stochastic dynamics in the
1910 c eigenspace of mass-scaled UNRES friction matrix
1913 gdt = fricgam(i) * d_time
1915 c Stochastic dynamics reduces to simple MD for zero friction
1917 if (gdt .le. zero) then
1918 pfric_vec(i) = 1.0d0
1919 vfric_vec(i) = d_time
1920 afric_vec(i) = 0.5d0 * d_time * d_time
1921 prand_vec(i) = 0.0d0
1922 vrand_vec1(i) = 0.0d0
1923 vrand_vec2(i) = 0.0d0
1925 c Analytical expressions when friction coefficient is large
1928 if (gdt .ge. gdt_radius) then
1931 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1932 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1933 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1934 vterm = 1.0d0 - egdt**2
1935 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1937 c Use series expansions when friction coefficient is small
1948 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1949 & - gdt5/120.0d0 + gdt6/720.0d0
1950 & - gdt7/5040.0d0 + gdt8/40320.0d0
1951 & - gdt9/362880.0d0) / fricgam(i)**2
1952 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1953 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1954 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1955 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1956 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1957 & + 127.0d0*gdt9/90720.0d0
1958 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1959 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1960 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1961 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1962 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1963 & - 17.0d0*gdt2/1280.0d0
1964 & + 17.0d0*gdt3/6144.0d0
1965 & + 40967.0d0*gdt4/34406400.0d0
1966 & - 57203.0d0*gdt5/275251200.0d0
1967 & - 1429487.0d0*gdt6/13212057600.0d0)
1970 c Compute the scaling factors of random terms for the nonzero friction case
1972 ktm = 0.5d0*d_time/fricgam(i)
1973 psig = dsqrt(ktm*pterm) / fricgam(i)
1974 vsig = dsqrt(ktm*vterm)
1975 rhoc = dsqrt(1.0d0 - rho*rho)
1977 vrand_vec1(i) = vsig * rho
1978 vrand_vec2(i) = vsig * rhoc
1983 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1986 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1987 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1991 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1994 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1995 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1996 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1997 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1998 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1999 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2002 t_sdsetup=t_sdsetup+MPI_Wtime()
2004 t_sdsetup=t_sdsetup+tcpu()-tt0
2008 c-------------------------------------------------------------
2009 subroutine eigtransf1(n,ndim,ab,d,c)
2012 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2018 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2024 c-------------------------------------------------------------
2025 subroutine eigtransf(n,ndim,a,b,d,c)
2028 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2034 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2040 c-------------------------------------------------------------
2041 subroutine sd_verlet1
2042 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2043 implicit real*8 (a-h,o-z)
2044 include 'DIMENSIONS'
2045 include 'COMMON.CONTROL'
2046 include 'COMMON.VAR'
2049 include 'COMMON.LANGEVIN'
2051 include 'COMMON.LANGEVIN.lang0'
2053 include 'COMMON.CHAIN'
2054 include 'COMMON.DERIV'
2055 include 'COMMON.GEO'
2056 include 'COMMON.LOCAL'
2057 include 'COMMON.INTERACT'
2058 include 'COMMON.IOUNITS'
2059 include 'COMMON.NAMES'
2060 double precision stochforcvec(MAXRES6)
2061 common /stochcalc/ stochforcvec
2062 logical lprn /.false./
2064 c write (iout,*) "dc_old"
2066 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2067 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2070 dc_work(j)=dc_old(j,0)
2071 d_t_work(j)=d_t_old(j,0)
2072 d_a_work(j)=d_a_old(j,0)
2077 dc_work(ind+j)=dc_old(j,i)
2078 d_t_work(ind+j)=d_t_old(j,i)
2079 d_a_work(ind+j)=d_a_old(j,i)
2084 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2086 dc_work(ind+j)=dc_old(j,i+nres)
2087 d_t_work(ind+j)=d_t_old(j,i+nres)
2088 d_a_work(ind+j)=d_a_old(j,i+nres)
2096 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2100 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2101 & vfric_mat(i,j),afric_mat(i,j),
2102 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2110 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2111 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2112 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2113 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2115 d_t_work_new(i)=ddt1+0.5d0*ddt2
2116 d_t_work(i)=ddt1+ddt2
2121 d_t(j,0)=d_t_work(j)
2126 dc(j,i)=dc_work(ind+j)
2127 d_t(j,i)=d_t_work(ind+j)
2132 if (itype(i).ne.10) then
2135 dc(j,inres)=dc_work(ind+j)
2136 d_t(j,inres)=d_t_work(ind+j)
2143 c--------------------------------------------------------------------------
2144 subroutine sd_verlet2
2145 c Calculating the adjusted velocities for accelerations
2146 implicit real*8 (a-h,o-z)
2147 include 'DIMENSIONS'
2148 include 'COMMON.CONTROL'
2149 include 'COMMON.VAR'
2152 include 'COMMON.LANGEVIN'
2154 include 'COMMON.LANGEVIN.lang0'
2156 include 'COMMON.CHAIN'
2157 include 'COMMON.DERIV'
2158 include 'COMMON.GEO'
2159 include 'COMMON.LOCAL'
2160 include 'COMMON.INTERACT'
2161 include 'COMMON.IOUNITS'
2162 include 'COMMON.NAMES'
2163 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2164 common /stochcalc/ stochforcvec
2166 c Compute the stochastic forces which contribute to velocity change
2168 call stochastic_force(stochforcvecV)
2175 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2176 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2177 & vrand_mat2(i,j)*stochforcvecV(j)
2179 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2183 d_t(j,0)=d_t_work(j)
2188 d_t(j,i)=d_t_work(ind+j)
2193 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2196 d_t(j,inres)=d_t_work(ind+j)
2203 c-----------------------------------------------------------
2204 subroutine sd_verlet_ciccotti_setup
2205 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2207 implicit real*8 (a-h,o-z)
2208 include 'DIMENSIONS'
2212 include 'COMMON.CONTROL'
2213 include 'COMMON.VAR'
2216 include 'COMMON.LANGEVIN'
2218 include 'COMMON.LANGEVIN.lang0'
2220 include 'COMMON.CHAIN'
2221 include 'COMMON.DERIV'
2222 include 'COMMON.GEO'
2223 include 'COMMON.LOCAL'
2224 include 'COMMON.INTERACT'
2225 include 'COMMON.IOUNITS'
2226 include 'COMMON.NAMES'
2227 include 'COMMON.TIME1'
2228 double precision emgdt(MAXRES6),
2229 & pterm,vterm,rho,rhoc,vsig,
2230 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2231 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2232 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2233 logical lprn /.false./
2234 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2235 double precision ktm
2242 c AL 8/17/04 Code adapted from tinker
2244 c Get the frictional and random terms for stochastic dynamics in the
2245 c eigenspace of mass-scaled UNRES friction matrix
2248 write (iout,*) "i",i," fricgam",fricgam(i)
2249 gdt = fricgam(i) * d_time
2251 c Stochastic dynamics reduces to simple MD for zero friction
2253 if (gdt .le. zero) then
2254 pfric_vec(i) = 1.0d0
2255 vfric_vec(i) = d_time
2256 afric_vec(i) = 0.5d0*d_time*d_time
2257 prand_vec(i) = afric_vec(i)
2258 vrand_vec2(i) = vfric_vec(i)
2260 c Analytical expressions when friction coefficient is large
2265 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2266 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2267 prand_vec(i) = afric_vec(i)
2268 vrand_vec2(i) = vfric_vec(i)
2270 c Compute the scaling factors of random terms for the nonzero friction case
2272 c ktm = 0.5d0*d_time/fricgam(i)
2273 c psig = dsqrt(ktm*pterm) / fricgam(i)
2274 c vsig = dsqrt(ktm*vterm)
2275 c prand_vec(i) = psig*afric_vec(i)
2276 c vrand_vec2(i) = vsig*vfric_vec(i)
2281 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2284 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2285 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2289 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2291 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2292 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2293 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2294 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2295 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2297 t_sdsetup=t_sdsetup+MPI_Wtime()
2299 t_sdsetup=t_sdsetup+tcpu()-tt0
2303 c-------------------------------------------------------------
2304 subroutine sd_verlet1_ciccotti
2305 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2306 implicit real*8 (a-h,o-z)
2307 include 'DIMENSIONS'
2311 include 'COMMON.CONTROL'
2312 include 'COMMON.VAR'
2315 include 'COMMON.LANGEVIN'
2317 include 'COMMON.LANGEVIN.lang0'
2319 include 'COMMON.CHAIN'
2320 include 'COMMON.DERIV'
2321 include 'COMMON.GEO'
2322 include 'COMMON.LOCAL'
2323 include 'COMMON.INTERACT'
2324 include 'COMMON.IOUNITS'
2325 include 'COMMON.NAMES'
2326 double precision stochforcvec(MAXRES6)
2327 common /stochcalc/ stochforcvec
2328 logical lprn /.false./
2330 c write (iout,*) "dc_old"
2332 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2333 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2336 dc_work(j)=dc_old(j,0)
2337 d_t_work(j)=d_t_old(j,0)
2338 d_a_work(j)=d_a_old(j,0)
2343 dc_work(ind+j)=dc_old(j,i)
2344 d_t_work(ind+j)=d_t_old(j,i)
2345 d_a_work(ind+j)=d_a_old(j,i)
2350 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2352 dc_work(ind+j)=dc_old(j,i+nres)
2353 d_t_work(ind+j)=d_t_old(j,i+nres)
2354 d_a_work(ind+j)=d_a_old(j,i+nres)
2363 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2367 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2368 & vfric_mat(i,j),afric_mat(i,j),
2369 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2377 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2378 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2379 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2380 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2382 d_t_work_new(i)=ddt1+0.5d0*ddt2
2383 d_t_work(i)=ddt1+ddt2
2388 d_t(j,0)=d_t_work(j)
2393 dc(j,i)=dc_work(ind+j)
2394 d_t(j,i)=d_t_work(ind+j)
2399 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2402 dc(j,inres)=dc_work(ind+j)
2403 d_t(j,inres)=d_t_work(ind+j)
2410 c--------------------------------------------------------------------------
2411 subroutine sd_verlet2_ciccotti
2412 c Calculating the adjusted velocities for accelerations
2413 implicit real*8 (a-h,o-z)
2414 include 'DIMENSIONS'
2415 include 'COMMON.CONTROL'
2416 include 'COMMON.VAR'
2419 include 'COMMON.LANGEVIN'
2421 include 'COMMON.LANGEVIN.lang0'
2423 include 'COMMON.CHAIN'
2424 include 'COMMON.DERIV'
2425 include 'COMMON.GEO'
2426 include 'COMMON.LOCAL'
2427 include 'COMMON.INTERACT'
2428 include 'COMMON.IOUNITS'
2429 include 'COMMON.NAMES'
2430 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2431 common /stochcalc/ stochforcvec
2433 c Compute the stochastic forces which contribute to velocity change
2435 call stochastic_force(stochforcvecV)
2442 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2443 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2444 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2446 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2450 d_t(j,0)=d_t_work(j)
2455 d_t(j,i)=d_t_work(ind+j)
2460 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2463 d_t(j,inres)=d_t_work(ind+j)