2 c------------------------------------------------
3 c The driver for molecular dynamics subroutines
4 c------------------------------------------------
11 include 'COMMON.SETUP'
12 include 'COMMON.CONTROL'
16 include 'COMMON.LAGRANGE.5diag'
18 include 'COMMON.LAGRANGE'
21 include 'COMMON.LANGEVIN'
24 include 'COMMON.LANGEVIN.lang0.5diag'
26 include 'COMMON.LANGEVIN.lang0'
29 include 'COMMON.CHAIN'
30 include 'COMMON.DERIV'
32 include 'COMMON.LOCAL'
33 include 'COMMON.INTERACT'
34 include 'COMMON.IOUNITS'
35 include 'COMMON.NAMES'
36 include 'COMMON.TIME1'
37 include 'COMMON.HAIRPIN'
38 double precision cm(3),L(3),vcm(3)
40 double precision v_work(maxres6),v_transf(maxres6)
48 integer i,j,icount_scale,itime_scal
49 integer nharp,iharp(4,maxres/3)
50 double precision scalfac
54 if (ilen(tmpdir).gt.0)
55 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
56 & //liczba(:ilen(liczba))//'.rst')
58 if (ilen(tmpdir).gt.0)
59 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
61 write (iout,*) "MD lang",lang
67 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
68 write (iout,*) "Restart frequency:",irest_freq
74 c Determine the inverse of the inertia matrix.
75 call setup_MD_matrices
79 t_MDsetup = MPI_Wtime()-tt0
81 t_MDsetup = tcpu()-tt0
84 c Entering the MD loop
90 if (lang.eq.2 .or. lang.eq.3) then
94 call sd_verlet_p_setup
96 call sd_verlet_ciccotti_setup
100 pfric0_mat(i,j,0)=pfric_mat(i,j)
101 afric0_mat(i,j,0)=afric_mat(i,j)
102 vfric0_mat(i,j,0)=vfric_mat(i,j)
103 prand0_mat(i,j,0)=prand_mat(i,j)
104 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
105 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
110 flag_stoch(i)=.false.
114 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
116 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
120 else if (lang.eq.1 .or. lang.eq.4) then
124 t_langsetup=MPI_Wtime()-tt0
127 t_langsetup=tcpu()-tt0
130 do itime=1,n_timestep
132 if (large.and. mod(itime,ntwe).eq.0)
133 & write (iout,*) "itime",itime
135 if (lang.gt.0 .and. surfarea .and.
136 & mod(itime,reset_fricmat).eq.0) then
137 if (lang.eq.2 .or. lang.eq.3) then
141 call sd_verlet_p_setup
143 call sd_verlet_ciccotti_setup
147 pfric0_mat(i,j,0)=pfric_mat(i,j)
148 afric0_mat(i,j,0)=afric_mat(i,j)
149 vfric0_mat(i,j,0)=vfric_mat(i,j)
150 prand0_mat(i,j,0)=prand_mat(i,j)
151 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
152 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
157 flag_stoch(i)=.false.
160 else if (lang.eq.1 .or. lang.eq.4) then
163 write (iout,'(a,i10)')
164 & "Friction matrix reset based on surface area, itime",itime
166 if (reset_vel .and. tbf .and. lang.eq.0
167 & .and. mod(itime,count_reset_vel).eq.0) then
169 write(iout,'(a,f20.2)')
170 & "Velocities reset to random values, time",totT
173 d_t_old(j,i)=d_t(j,i)
177 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
181 d_t(j,0)=d_t(j,0)-vcm(j)
184 kinetic_T=2.0d0/(dimen3*Rb)*EK
185 scalfac=dsqrt(T_bath/kinetic_T)
186 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
189 d_t_old(j,i)=scalfac*d_t(j,i)
195 c Time-reversible RESPA algorithm
196 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
197 call RESPA_step(itime)
199 c Variable time step algorithm.
200 call velverlet_step(itime)
204 call brown_step(itime)
206 print *,"Brown dynamics not here!"
208 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
215 if (mod(itime,ntwe).eq.0) then
217 C call enerprint(potEcomp)
218 C print *,itime,'AFM',Eafmforc,etot
232 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
235 v_work(ind)=d_t(j,i+nres)
240 write (66,'(80f10.5)')
241 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
245 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
247 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
249 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
252 if (mod(itime,ntwx).eq.0) then
253 c write(iout,*) 'time=',itime
254 C call check_ecartint
256 write (tytul,'("time",f8.2)') totT
258 call hairpin(.true.,nharp,iharp)
259 call secondary2(.true.)
260 call pdbout(potE,tytul,ipdb)
265 if (rstcount.eq.irest_freq.or.itime.eq.n_timestep) then
266 open(irest2,file=rest2name,status='unknown')
267 write(irest2,*) totT,EK,potE,totE,t_bath
269 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
272 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
283 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
285 & 'MD calculations setup:',t_MDsetup,
286 & 'Energy & gradient evaluation:',t_enegrad,
287 & 'Stochastic MD setup:',t_langsetup,
288 & 'Stochastic MD step setup:',t_sdsetup,
290 write (iout,'(/28(1h=),a25,27(1h=))')
291 & ' End of MD calculation '
293 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
294 & " eshort",t_eshort," list",time_list
295 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
296 & " time_fricmatmult",time_fricmatmult," time_fsample ",
301 c-------------------------------------------------------------------------------
302 subroutine velverlet_step(itime)
303 c-------------------------------------------------------------------------------
304 c Perform a single velocity Verlet step; the time step can be rescaled if
305 c increments in accelerations exceed the threshold
306 c-------------------------------------------------------------------------------
311 integer ierror,ierrcode,errcode
313 include 'COMMON.SETUP'
314 include 'COMMON.CONTROL'
318 include 'COMMON.LAGRANGE.5diag'
320 include 'COMMON.LAGRANGE'
323 include 'COMMON.LANGEVIN'
326 include 'COMMON.LANGEVIN.lang0.5diag'
328 include 'COMMON.LANGEVIN.lang0'
331 include 'COMMON.CHAIN'
332 include 'COMMON.DERIV'
334 include 'COMMON.LOCAL'
335 include 'COMMON.INTERACT'
336 include 'COMMON.IOUNITS'
337 include 'COMMON.NAMES'
338 include 'COMMON.TIME1'
339 include 'COMMON.MUCA'
340 double precision vcm(3),incr(3)
341 double precision cm(3),L(3)
342 integer ilen,count,rstcount
345 integer maxcount_scale /20/
347 double precision stochforcvec(MAXRES6)
348 common /stochcalc/ stochforcvec
349 integer itime,icount_scale,itime_scal,ifac_time,i,j,itt
351 double precision epdrift,fac_time
358 else if (lang.eq.2 .or. lang.eq.3) then
360 call stochastic_force(stochforcvec)
363 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
365 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
372 icount_scale=icount_scale+1
373 if (icount_scale.gt.maxcount_scale) then
375 & "ERROR: too many attempts at scaling down the time step. ",
376 & "amax=",amax,"epdrift=",epdrift,
377 & "damax=",damax,"edriftmax=",edriftmax,
381 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
385 c First step of the velocity Verlet algorithm
390 else if (lang.eq.3) then
392 call sd_verlet1_ciccotti
394 else if (lang.eq.1) then
399 c Build the chain from the newly calculated coordinates
401 if (rattle) call rattle1
403 if (large.and. mod(itime,ntwe).eq.0) then
404 write (iout,*) "Cartesian and internal coordinates: step 1"
409 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
410 & (dc(j,i+nres),j=1,3)
412 write (iout,*) "Accelerations"
414 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
415 & (d_a(j,i+nres),j=1,3)
417 write (iout,*) "Velocities, step 1"
419 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
420 & (d_t(j,i+nres),j=1,3)
429 c Calculate energy and forces
431 call etotal(potEcomp)
432 ! AL 4/17/17: Reduce the steps if NaNs occurred.
433 if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0))) then
438 if (large.and. mod(itime,ntwe).eq.0)
439 & call enerprint(potEcomp)
442 t_etotal=t_etotal+MPI_Wtime()-tt0
444 t_etotal=t_etotal+tcpu()-tt0
447 potE=potEcomp(0)-potEcomp(27)
449 c Get the new accelerations
452 t_enegrad=t_enegrad+MPI_Wtime()-tt0
454 t_enegrad=t_enegrad+tcpu()-tt0
456 c Determine maximum acceleration and scale down the timestep if needed
458 amax=amax/(itime_scal+1)**2
459 call predict_edrift(epdrift)
460 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
461 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
463 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
465 itime_scal=itime_scal+ifac_time
466 c fac_time=dmin1(damax/amax,0.5d0)
467 fac_time=0.5d0**ifac_time
468 d_time=d_time*fac_time
469 if (lang.eq.2 .or. lang.eq.3) then
471 c write (iout,*) "Calling sd_verlet_setup: 1"
472 c Rescale the stochastic forces and recalculate or restore
473 c the matrices of tinker integrator
474 if (itime_scal.gt.maxflag_stoch) then
475 if (large) write (iout,'(a,i5,a)')
476 & "Calculate matrices for stochastic step;",
477 & " itime_scal ",itime_scal
479 call sd_verlet_p_setup
481 call sd_verlet_ciccotti_setup
483 write (iout,'(2a,i3,a,i3,1h.)')
484 & "Warning: cannot store matrices for stochastic",
485 & " integration because the index",itime_scal,
486 & " is greater than",maxflag_stoch
487 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
488 & " integration Langevin algorithm for better efficiency."
489 else if (flag_stoch(itime_scal)) then
490 if (large) write (iout,'(a,i5,a,l1)')
491 & "Restore matrices for stochastic step; itime_scal ",
492 & itime_scal," flag ",flag_stoch(itime_scal)
495 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
496 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
497 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
498 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
499 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
500 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
504 if (large) write (iout,'(2a,i5,a,l1)')
505 & "Calculate & store matrices for stochastic step;",
506 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
508 call sd_verlet_p_setup
510 call sd_verlet_ciccotti_setup
512 flag_stoch(ifac_time)=.true.
515 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
516 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
517 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
518 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
519 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
520 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
524 fac_time=1.0d0/dsqrt(fac_time)
526 stochforcvec(i)=fac_time*stochforcvec(i)
529 else if (lang.eq.1) then
530 c Rescale the accelerations due to stochastic forces
531 fac_time=1.0d0/dsqrt(fac_time)
533 d_as_work(i)=d_as_work(i)*fac_time
536 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
537 & "itime",itime," Timestep scaled down to ",
538 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
540 c Second step of the velocity Verlet algorithm
545 else if (lang.eq.3) then
547 call sd_verlet2_ciccotti
549 else if (lang.eq.1) then
554 if (rattle) call rattle2
557 C print *,totTafm,"TU?"
558 if (d_time.ne.d_time0) then
561 if (lang.eq.2 .or. lang.eq.3) then
562 if (large) write (iout,'(a)')
563 & "Restore original matrices for stochastic step"
564 c write (iout,*) "Calling sd_verlet_setup: 2"
565 c Restore the matrices of tinker integrator if the time step has been restored
568 pfric_mat(i,j)=pfric0_mat(i,j,0)
569 afric_mat(i,j)=afric0_mat(i,j,0)
570 vfric_mat(i,j)=vfric0_mat(i,j,0)
571 prand_mat(i,j)=prand0_mat(i,j,0)
572 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
573 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
582 c Calculate the kinetic and the total energy and the kinetic temperature
587 c write (iout,*) "step",itime," EK",EK," EK1",EK1
589 c Couple the system to Berendsen bath if needed
590 if (tbf .and. lang.eq.0) then
593 kinetic_T=2.0d0/(dimen3*Rb)*EK
594 c Backup the coordinates, velocities, and accelerations
598 d_t_old(j,i)=d_t(j,i)
599 d_a_old(j,i)=d_a(j,i)
603 if (mod(itime,ntwe).eq.0 .and. large) then
604 write (iout,*) "Velocities, step 2"
606 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
607 & (d_t(j,i+nres),j=1,3)
613 c-------------------------------------------------------------------------------
614 subroutine RESPA_step(itime)
615 c-------------------------------------------------------------------------------
616 c Perform a single RESPA step.
617 c-------------------------------------------------------------------------------
622 integer IERROR,ERRCODE
624 include 'COMMON.SETUP'
625 include 'COMMON.CONTROL'
629 include 'COMMON.LAGRANGE.5diag'
631 include 'COMMON.LAGRANGE'
634 include 'COMMON.LANGEVIN'
637 include 'COMMON.LANGEVIN.lang0.5diag'
639 include 'COMMON.LANGEVIN.lang0'
642 include 'COMMON.CHAIN'
643 include 'COMMON.DERIV'
645 include 'COMMON.LOCAL'
646 include 'COMMON.INTERACT'
647 include 'COMMON.IOUNITS'
648 include 'COMMON.NAMES'
649 include 'COMMON.TIME1'
651 common /shortcheck/ lprint_short
652 double precision energia_short(0:n_ene),
653 & energia_long(0:n_ene)
654 double precision cm(3),L(3),vcm(3),incr(3)
655 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
656 & d_a_old0(3,0:maxres2)
658 double precision fac_time
659 logical PRINT_AMTS_MSG /.false./
660 integer ilen,count,rstcount
663 integer maxcount_scale /10/
665 double precision stochforcvec(MAXRES6)
666 common /stochcalc/ stochforcvec
670 common /cipiszcze/ itt
672 double precision epdrift,epdriftmax
676 if (large.and. mod(itime,ntwe).eq.0) then
677 write (iout,*) "***************** RESPA itime",itime
678 write (iout,*) "Cartesian and internal coordinates: step 0"
683 write (iout,*) "Accelerations from long-range forces"
685 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
686 & (d_a(j,i+nres),j=1,3)
688 write (iout,*) "Velocities, step 0"
690 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
691 & (d_t(j,i+nres),j=1,3)
696 c Perform the initial RESPA step (increment velocities)
697 c write (iout,*) "*********************** RESPA ini"
700 if (mod(itime,ntwe).eq.0 .and. large) then
701 write (iout,*) "Velocities, end"
703 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
704 & (d_t(j,i+nres),j=1,3)
708 c Compute the short-range forces
714 C 7/2/2009 commented out
716 c call etotal_short(energia_short)
719 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
722 d_a(j,i)=d_a_short(j,i)
726 if (large.and. mod(itime,ntwe).eq.0) then
727 write (iout,*) "energia_short",energia_short(0)
728 write (iout,*) "Accelerations from short-range forces"
730 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
731 & (d_a(j,i+nres),j=1,3)
736 t_enegrad=t_enegrad+MPI_Wtime()-tt0
738 t_enegrad=t_enegrad+tcpu()-tt0
743 d_t_old(j,i)=d_t(j,i)
744 d_a_old(j,i)=d_a(j,i)
747 c 6/30/08 A-MTS: attempt at increasing the split number
750 dc_old0(j,i)=dc_old(j,i)
751 d_t_old0(j,i)=d_t_old(j,i)
752 d_a_old0(j,i)=d_a_old(j,i)
755 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
756 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
763 c write (iout,*) "itime",itime," ntime_split",ntime_split
764 c Split the time step
765 d_time=d_time0/ntime_split
766 c Perform the short-range RESPA steps (velocity Verlet increments of
767 c positions and velocities using short-range forces)
768 c write (iout,*) "*********************** RESPA split"
769 do itsplit=1,ntime_split
772 else if (lang.eq.2 .or. lang.eq.3) then
774 call stochastic_force(stochforcvec)
777 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
779 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
784 c First step of the velocity Verlet algorithm
789 else if (lang.eq.3) then
791 call sd_verlet1_ciccotti
793 else if (lang.eq.1) then
798 c Build the chain from the newly calculated coordinates
800 if (rattle) call rattle1
802 if (large.and. mod(itime,ntwe).eq.0) then
803 write (iout,*) "***** ITSPLIT",itsplit
804 write (iout,*) "Cartesian and internal coordinates: step 1"
809 write (iout,*) "Velocities, step 1"
811 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
812 & (d_t(j,i+nres),j=1,3)
821 c Calculate energy and forces
822 c if (large.and. mod(itime,ntwe).eq.0) lprint_short=.true.
824 call etotal_short(energia_short)
825 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
826 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) )
828 if (PRINT_AMTS_MSG) write (iout,*)
829 & "Infinities/NaNs in energia_short",
830 & energia_short(0),"; increasing ntime_split to",ntime_split
831 ntime_split=ntime_split*2
832 if (ntime_split.gt.maxtime_split) then
834 write (iout,*) "Cannot rescue the run; aborting job.",
835 & " Retry with a smaller time step"
837 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
839 write (iout,*) "Cannot rescue the run; terminating.",
840 & " Retry with a smaller time step"
846 if (large.and. mod(itime,ntwe).eq.0)
847 & call enerprint(energia_short)
850 t_eshort=t_eshort+MPI_Wtime()-tt0
852 t_eshort=t_eshort+tcpu()-tt0
856 c Get the new accelerations
858 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
861 d_a_short(j,i)=d_a(j,i)
865 if (large.and. mod(itime,ntwe).eq.0) then
866 write (iout,*)"energia_short",energia_short(0)
867 write (iout,*) "Accelerations from short-range forces"
869 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
870 & (d_a(j,i+nres),j=1,3)
875 c Determine maximum acceleration and scale down the timestep if needed
877 amax=amax/ntime_split**2
878 call predict_edrift(epdrift)
879 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
880 & write (iout,*) "amax",amax," damax",damax,
881 & " epdrift",epdrift," epdriftmax",epdriftmax
882 c Exit loop and try with increased split number if the change of
883 c acceleration is too big
884 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
885 if (ntime_split.lt.maxtime_split) then
887 ntime_split=ntime_split*2
890 dc_old(j,i)=dc_old0(j,i)
891 d_t_old(j,i)=d_t_old0(j,i)
892 d_a_old(j,i)=d_a_old0(j,i)
895 if (PRINT_AMTS_MSG) then
896 write (iout,*) "acceleration/energy drift too large",amax,
897 & epdrift," split increased to ",ntime_split," itime",itime,
903 & "Uh-hu. Bumpy landscape. Maximum splitting number",
905 & " already reached!!! Trying to carry on!"
909 t_enegrad=t_enegrad+MPI_Wtime()-tt0
911 t_enegrad=t_enegrad+tcpu()-tt0
913 c Second step of the velocity Verlet algorithm
918 else if (lang.eq.3) then
920 call sd_verlet2_ciccotti
922 else if (lang.eq.1) then
927 if (rattle) call rattle2
928 c Backup the coordinates, velocities, and accelerations
932 d_t_old(j,i)=d_t(j,i)
933 d_a_old(j,i)=d_a(j,i)
940 c Restore the time step
942 c Compute long-range forces
949 call etotal_long(energia_long)
950 ! AL 4/17/2017 Handling NaNs
951 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
954 & "Infinitied/NaNs in energia_long, Aborting MPI job."
955 call enerprint(energia_long)
957 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
959 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
960 call enerprint(energia_long)
965 c lprint_short=.false.
966 if (large.and. mod(itime,ntwe).eq.0)
967 & call enerprint(energia_long)
970 t_elong=t_elong+MPI_Wtime()-tt0
972 t_elong=t_elong+tcpu()-tt0
978 t_enegrad=t_enegrad+MPI_Wtime()-tt0
980 t_enegrad=t_enegrad+tcpu()-tt0
982 c Compute accelerations from long-range forces
984 if (large.and. mod(itime,ntwe).eq.0) then
985 write (iout,*) "energia_long",energia_long(0)
986 write (iout,*) "Cartesian and internal coordinates: step 2"
991 write (iout,*) "Accelerations from long-range forces"
993 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
994 & (d_a(j,i+nres),j=1,3)
996 write (iout,*) "Velocities, step 2"
998 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
999 & (d_t(j,i+nres),j=1,3)
1003 c Compute the final RESPA step (increment velocities)
1004 c write (iout,*) "*********************** RESPA fin"
1006 c Compute the complete potential energy
1008 potEcomp(i)=energia_short(i)+energia_long(i)
1010 potE=potEcomp(0)-potEcomp(27)
1012 if (large.and. mod(itime,ntwe).eq.0) then
1013 call enerprint(potEcomp)
1014 write (iout,*) "potE",potE
1017 c potE=energia_short(0)+energia_long(0)
1020 c Calculate the kinetic and the total energy and the kinetic temperature
1023 c Couple the system to Berendsen bath if needed
1024 if (tbf .and. lang.eq.0) then
1027 kinetic_T=2.0d0/(dimen3*Rb)*EK
1028 c Backup the coordinates, velocities, and accelerations
1030 if (mod(itime,ntwe).eq.0 .and. large) then
1031 write (iout,*) "Velocities, end"
1033 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1034 & (d_t(j,i+nres),j=1,3)
1040 c---------------------------------------------------------------------
1041 subroutine RESPA_vel
1042 c First and last RESPA step (incrementing velocities using long-range
1045 include 'DIMENSIONS'
1046 include 'COMMON.CONTROL'
1047 include 'COMMON.VAR'
1050 include 'COMMON.LAGRANGE.5diag'
1052 include 'COMMON.LAGRANGE'
1054 include 'COMMON.CHAIN'
1055 include 'COMMON.DERIV'
1056 include 'COMMON.GEO'
1057 include 'COMMON.LOCAL'
1058 include 'COMMON.INTERACT'
1059 include 'COMMON.IOUNITS'
1060 include 'COMMON.NAMES'
1063 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1067 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1071 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1074 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1080 c-----------------------------------------------------------------
1082 c Applying velocity Verlet algorithm - step 1 to coordinates
1084 include 'DIMENSIONS'
1085 include 'COMMON.CONTROL'
1086 include 'COMMON.VAR'
1089 include 'COMMON.LAGRANGE.5diag'
1091 include 'COMMON.LAGRANGE'
1093 include 'COMMON.CHAIN'
1094 include 'COMMON.DERIV'
1095 include 'COMMON.GEO'
1096 include 'COMMON.LOCAL'
1097 include 'COMMON.INTERACT'
1098 include 'COMMON.IOUNITS'
1099 include 'COMMON.NAMES'
1100 double precision adt,adt2
1103 write (iout,*) "VELVERLET1 START: DC"
1105 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1106 & (dc(j,i+nres),j=1,3)
1110 adt=d_a_old(j,0)*d_time
1112 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1113 d_t_new(j,0)=d_t_old(j,0)+adt2
1114 d_t(j,0)=d_t_old(j,0)+adt
1120 write (iout,*) "i",i," d_a_old",(d_a_old(j,i),j=1,3)
1123 adt=d_a_old(j,i)*d_time
1125 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1126 d_t_new(j,i)=d_t_old(j,i)+adt2
1127 d_t(j,i)=d_t_old(j,i)+adt
1132 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1135 adt=d_a_old(j,inres)*d_time
1137 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1138 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1139 d_t(j,inres)=d_t_old(j,inres)+adt
1144 write (iout,*) "VELVERLET1 END: DC"
1146 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1147 & (dc(j,i+nres),j=1,3)
1152 c---------------------------------------------------------------------
1154 c Step 2 of the velocity Verlet algorithm: update velocities
1156 include 'DIMENSIONS'
1157 include 'COMMON.CONTROL'
1158 include 'COMMON.VAR'
1161 include 'COMMON.LAGRANGE.5diag'
1163 include 'COMMON.LAGRANGE'
1165 include 'COMMON.CHAIN'
1166 include 'COMMON.DERIV'
1167 include 'COMMON.GEO'
1168 include 'COMMON.LOCAL'
1169 include 'COMMON.INTERACT'
1170 include 'COMMON.IOUNITS'
1171 include 'COMMON.NAMES'
1174 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1178 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1182 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1185 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1191 c-----------------------------------------------------------------
1192 subroutine sddir_precalc
1193 c Applying velocity Verlet algorithm - step 1 to coordinates
1195 include 'DIMENSIONS'
1199 include 'COMMON.CONTROL'
1200 include 'COMMON.VAR'
1203 include 'COMMON.LAGRANGE.5diag'
1205 include 'COMMON.LAGRANGE'
1208 include 'COMMON.LANGEVIN'
1211 include 'COMMON.LANGEVIN.lang0.5diag'
1213 include 'COMMON.LANGEVIN.lang0'
1216 include 'COMMON.CHAIN'
1217 include 'COMMON.DERIV'
1218 include 'COMMON.GEO'
1219 include 'COMMON.LOCAL'
1220 include 'COMMON.INTERACT'
1221 include 'COMMON.IOUNITS'
1222 include 'COMMON.NAMES'
1223 include 'COMMON.TIME1'
1224 double precision time00
1225 double precision stochforcvec(MAXRES6)
1226 common /stochcalc/ stochforcvec
1229 c Compute friction and stochastic forces
1237 c write (iout,*) "After friction_force"
1239 time_fric=time_fric+MPI_Wtime()-time00
1242 time_fric=time_fric+tcpu()-time00
1245 call stochastic_force(stochforcvec)
1246 c write (iout,*) "After stochastic_force"
1248 time_stoch=time_stoch+MPI_Wtime()-time00
1250 time_stoch=time_stoch+tcpu()-time00
1253 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1254 c forces (d_as_work)
1257 c write (iout,*) "friction accelerations"
1258 call fivediaginv_mult(dimen,fric_work, d_af_work)
1259 c write (iout,*) "stochastic acceleratios"
1260 call fivediaginv_mult(dimen,stochforcvec, d_as_work)
1262 call ginv_mult(fric_work, d_af_work)
1263 call ginv_mult(stochforcvec, d_as_work)
1266 write (iout,*) "d_af_work"
1267 write (iout,'(3f10.5)') (d_af_work(i),i=1,dimen3)
1268 write (iout,*) "d_as_work"
1269 write (iout,'(3f10.5)') (d_as_work(i),i=1,dimen3)
1270 write (iout,*) "Leaving sddir_precalc"
1274 c---------------------------------------------------------------------
1275 subroutine sddir_verlet1
1276 c Applying velocity Verlet algorithm - step 1 to velocities
1278 include 'DIMENSIONS'
1279 include 'COMMON.CONTROL'
1280 include 'COMMON.VAR'
1283 include 'COMMON.LAGRANGE.5diag'
1285 include 'COMMON.LAGRANGE'
1288 include 'COMMON.LANGEVIN'
1291 include 'COMMON.LANGEVIN.lang0.5diag'
1293 include 'COMMON.LANGEVIN.lang0'
1296 include 'COMMON.CHAIN'
1297 include 'COMMON.DERIV'
1298 include 'COMMON.GEO'
1299 include 'COMMON.LOCAL'
1300 include 'COMMON.INTERACT'
1301 include 'COMMON.IOUNITS'
1302 include 'COMMON.NAMES'
1303 c Revised 3/31/05 AL: correlation between random contributions to
1304 c position and velocity increments included.
1305 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1306 double precision adt,adt2
1307 integer i,j,ind,inres
1309 c Add the contribution from BOTH friction and stochastic force to the
1310 c coordinates, but ONLY the contribution from the friction forces to velocities
1313 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1314 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1315 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1316 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1317 d_t(j,0)=d_t_old(j,0)+adt
1322 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1323 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1324 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1325 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1326 d_t(j,i)=d_t_old(j,i)+adt
1331 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1334 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1335 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1336 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1337 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1338 d_t(j,inres)=d_t_old(j,inres)+adt
1345 c---------------------------------------------------------------------
1346 subroutine sddir_verlet2
1347 c Calculating the adjusted velocities for accelerations
1349 include 'DIMENSIONS'
1350 include 'COMMON.CONTROL'
1351 include 'COMMON.VAR'
1354 include 'COMMON.LAGRANGE.5diag'
1356 include 'COMMON.LAGRANGE'
1359 include 'COMMON.LANGEVIN'
1362 include 'COMMON.LANGEVIN.lang0.5diag'
1364 include 'COMMON.LANGEVIN.lang0'
1367 include 'COMMON.CHAIN'
1368 include 'COMMON.DERIV'
1369 include 'COMMON.GEO'
1370 include 'COMMON.LOCAL'
1371 include 'COMMON.INTERACT'
1372 include 'COMMON.IOUNITS'
1373 include 'COMMON.NAMES'
1374 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1375 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1376 integer i,j,inres,ind
1377 c Revised 3/31/05 AL: correlation between random contributions to
1378 c position and velocity increments included.
1379 c The correlation coefficients are calculated at low-friction limit.
1380 c Also, friction forces are now not calculated with new velocities.
1382 c call friction_force
1383 call stochastic_force(stochforcvec)
1385 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1386 c forces (d_as_work)
1389 call fivediaginv_mult(maxres6,stochforcvec, d_as_work1)
1391 call ginv_mult(stochforcvec, d_as_work1)
1397 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1398 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1403 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1404 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1409 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1412 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1413 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1414 & +cos60*d_as_work1(ind+j))*d_time
1421 c---------------------------------------------------------------------
1422 subroutine max_accel
1424 c Find the maximum difference in the accelerations of the the sites
1425 c at the beginning and the end of the time step.
1428 include 'DIMENSIONS'
1429 include 'COMMON.CONTROL'
1430 include 'COMMON.VAR'
1433 include 'COMMON.LAGRANGE.5diag'
1435 include 'COMMON.LAGRANGE'
1437 include 'COMMON.CHAIN'
1438 include 'COMMON.DERIV'
1439 include 'COMMON.GEO'
1440 include 'COMMON.LOCAL'
1441 include 'COMMON.INTERACT'
1442 include 'COMMON.IOUNITS'
1444 double precision aux(3),accel(3),accel_old(3),dacc
1446 c aux(j)=d_a(j,0)-d_a_old(j,0)
1447 accel_old(j)=d_a_old(j,0)
1454 c 7/3/08 changed to asymmetric difference
1456 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1457 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1458 accel(j)=accel(j)+0.5d0*d_a(j,i)
1459 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1460 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1461 dacc=dabs(accel(j)-accel_old(j))
1462 c write (iout,*) i,dacc
1463 if (dacc.gt.amax) amax=dacc
1471 accel_old(j)=d_a_old(j,0)
1476 accel_old(j)=accel_old(j)+d_a_old(j,1)
1477 accel(j)=accel(j)+d_a(j,1)
1481 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1483 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1484 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1485 accel(j)=accel(j)+d_a(j,i+nres)
1489 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1490 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1491 dacc=dabs(accel(j)-accel_old(j))
1492 c write (iout,*) "side-chain",i,dacc
1493 if (dacc.gt.amax) amax=dacc
1497 accel_old(j)=accel_old(j)+d_a_old(j,i)
1498 accel(j)=accel(j)+d_a(j,i)
1499 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1504 c---------------------------------------------------------------------
1505 subroutine predict_edrift(epdrift)
1507 c Predict the drift of the potential energy
1510 include 'DIMENSIONS'
1511 include 'COMMON.CONTROL'
1512 include 'COMMON.VAR'
1515 include 'COMMON.LAGRANGE.5diag'
1517 include 'COMMON.LAGRANGE'
1519 include 'COMMON.CHAIN'
1520 include 'COMMON.DERIV'
1521 include 'COMMON.GEO'
1522 include 'COMMON.LOCAL'
1523 include 'COMMON.INTERACT'
1524 include 'COMMON.IOUNITS'
1525 include 'COMMON.MUCA'
1526 double precision epdrift,epdriftij
1528 c Drift of the potential energy
1534 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1535 if (lmuca) epdriftij=epdriftij*factor
1536 c write (iout,*) "back",i,j,epdriftij
1537 if (epdriftij.gt.epdrift) epdrift=epdriftij
1541 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1544 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1545 if (lmuca) epdriftij=epdriftij*factor
1546 c write (iout,*) "side",i,j,epdriftij
1547 if (epdriftij.gt.epdrift) epdrift=epdriftij
1551 epdrift=0.5d0*epdrift*d_time*d_time
1552 c write (iout,*) "epdrift",epdrift
1555 c-----------------------------------------------------------------------
1556 subroutine verlet_bath
1558 c Coupling to the thermostat by using the Berendsen algorithm
1561 include 'DIMENSIONS'
1562 include 'COMMON.CONTROL'
1563 include 'COMMON.VAR'
1566 include 'COMMON.LAGRANGE.5diag'
1568 include 'COMMON.LAGRANGE'
1572 include 'COMMON.LANGEVIN.lang0.5diag'
1574 include 'COMMON.LANGEVIN.lang0'
1577 include 'COMMON.LANGEVIN'
1579 include 'COMMON.CHAIN'
1580 include 'COMMON.DERIV'
1581 include 'COMMON.GEO'
1582 include 'COMMON.LOCAL'
1583 include 'COMMON.INTERACT'
1584 include 'COMMON.IOUNITS'
1585 include 'COMMON.NAMES'
1586 double precision T_half,fact
1588 T_half=2.0d0/(dimen3*Rb)*EK
1589 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1590 c write(iout,*) "T_half", T_half
1591 c write(iout,*) "EK", EK
1592 c write(iout,*) "fact", fact
1594 d_t(j,0)=fact*d_t(j,0)
1598 d_t(j,i)=fact*d_t(j,i)
1602 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1605 d_t(j,inres)=fact*d_t(j,inres)
1611 c---------------------------------------------------------
1613 c Set up the initial conditions of a MD simulation
1615 include 'DIMENSIONS'
1619 integer IERROR,ERRCODE,error_msg,ierr,ierrcode
1621 include 'COMMON.SETUP'
1622 include 'COMMON.CONTROL'
1623 include 'COMMON.VAR'
1626 include 'COMMON.LAGRANGE.5diag'
1628 include 'COMMON.LAGRANGE'
1630 include 'COMMON.QRESTR'
1632 include 'COMMON.LANGEVIN'
1635 include 'COMMON.LANGEVIN.lang0.5diag'
1638 include 'COMMON.LANGEVIN.lang0.5diag'
1640 include 'COMMON.LANGEVIN.lang0'
1644 include 'COMMON.CHAIN'
1645 include 'COMMON.DERIV'
1646 include 'COMMON.GEO'
1647 include 'COMMON.LOCAL'
1648 include 'COMMON.INTERACT'
1649 include 'COMMON.IOUNITS'
1650 include 'COMMON.NAMES'
1651 include 'COMMON.REMD'
1652 include 'COMMON.TIME1'
1656 common /lbfgstat/ status,niter,nfun
1658 integer n_model_try,list_model_try(max_template),k
1659 double precision tt0
1660 real*8 energia_long(0:n_ene),
1661 & energia_short(0:n_ene),vcm(3),incr(3)
1662 double precision cm(3),L(3),xv,sigv,lowb,highb
1663 double precision varia(maxvar),energia(0:n_ene)
1670 integer i,ipos,iq,iw,j,iranmin,nft_sc,iretcode,nfun,itrial,itmp,
1673 double precision etot
1675 integer i_start_models(0:nodes-1)
1676 double precision potEcomp_long(0:Max_Ene)
1677 write (iout,*) "init_MD INDPDB",indpdb
1679 c write(iout,*) "d_time", d_time
1680 c Compute the standard deviations of stochastic forces for Langevin dynamics
1681 c if the friction coefficients do not depend on surface area
1682 if (lang.gt.0 .and. .not.surfarea) then
1684 stdforcp(i)=stdfp*dsqrt(gamp)
1687 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1688 & *dsqrt(gamsc(iabs(itype(i))))
1691 c Open the pdb file for snapshotshots
1694 if (ilen(tmpdir).gt.0)
1695 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1696 & liczba(:ilen(liczba))//".pdb")
1698 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1702 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1703 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1704 & liczba(:ilen(liczba))//".x")
1705 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1708 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1709 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1710 & liczba(:ilen(liczba))//".cx")
1711 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1717 if (ilen(tmpdir).gt.0)
1718 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1719 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1721 if (ilen(tmpdir).gt.0)
1722 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1723 cartname=prefix(:ilen(prefix))//"_MD.cx"
1727 write (qstr,'(256(1h ))')
1730 iq = qinfrag(i,iset)*10
1731 iw = wfrag(i,iset)/100
1733 if(me.eq.king.or..not.out1file)
1734 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1735 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1740 iq = qinpair(i,iset)*10
1741 iw = wpair(i,iset)/100
1743 if(me.eq.king.or..not.out1file)
1744 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1745 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1749 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1751 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1753 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1755 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1758 write (iout,*) "REST ",rest
1760 if (restart1file) then
1762 & inquire(file=mremd_rst_name,exist=file_exist)
1764 c write (*,*) me," Before broadcast: file_exist",file_exist
1765 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1767 c write (*,*) me," After broadcast: file_exist",file_exist
1768 c inquire(file=mremd_rst_name,exist=file_exist)
1770 if(me.eq.king.or..not.out1file)
1771 & write(iout,*) "Initial state read by master and distributed"
1773 if (ilen(tmpdir).gt.0)
1774 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1775 & //liczba(:ilen(liczba))//'.rst')
1776 inquire(file=rest2name,exist=file_exist)
1779 if(.not.restart1file) then
1780 if(me.eq.king.or..not.out1file)
1781 & write(iout,*) "Initial state will be read from file ",
1782 & rest2name(:ilen(rest2name))
1785 call rescale_weights(t_bath)
1788 if(me.eq.king.or..not.out1file)then
1789 if (restart1file) then
1790 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1793 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1796 write(iout,*) "Initial velocities randomly generated"
1803 c Generate initial velocities
1804 if(me.eq.king.or..not.out1file)
1805 & write(iout,*) "Initial velocities randomly generated"
1808 CtotTafm is the variable for AFM time which eclipsed during
1811 c rest2name = prefix(:ilen(prefix))//'.rst'
1812 if(me.eq.king.or..not.out1file)then
1813 write (iout,*) "Initial velocities"
1815 write (iout,'(i7,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1816 & (d_t(j,i+nres),j=1,3)
1820 c Zeroing the total angular momentum of the system
1821 write(iout,*) "Calling the zero-angular
1822 & momentum subroutine"
1824 c Getting the potential energy and forces and velocities and accelerations
1826 write (iout,*) "velocity of the center of the mass:"
1827 write (iout,*) (vcm(j),j=1,3)
1830 d_t(j,0)=d_t(j,0)-vcm(j)
1832 c Removing the velocity of the center of mass
1834 if(me.eq.king.or..not.out1file)then
1835 write (iout,*) "vcm right after adjustment:"
1836 write (iout,*) (vcm(j),j=1,3)
1837 write (iout,*) "Initial velocities after adjustment"
1839 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1840 & (d_t(j,i+nres),j=1,3)
1845 c write (iout,*) "init_MD before initial structure REST ",rest
1846 if(start_from_model .and. (me.eq.king .or. .not. out1file))
1847 & write(iout,*) 'START_FROM_MODELS is ON'
1848 if(start_from_model .and. rest) then
1849 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1851 & 'START_FROM_MODELS is OFF because the run is restarted'
1852 write(iout,*) 'Remove restart keyword from input'
1855 c write (iout,*) "rest ",rest," start_from_model",start_from_model,
1856 c & " nmodel_start",nmodel_start," preminim",preminim
1859 if (iranconf.ne.0) then
1860 c 8/22/17 AL Loop to produce a low-energy random conformation
1864 print *, 'Calling OVERLAP_SC'
1865 call overlap_sc(fail)
1869 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1870 print *,'SC_move',nft_sc,etot
1871 if(me.eq.king.or..not.out1file)
1872 & write(iout,*) 'SC_move',nft_sc,etot
1876 if (me.eq.king.or..not.out1file) write(iout,*)
1877 & 'Minimizing random structure: Calling MINIM_DC'
1878 call minim_dc(etot,iretcode,nfun)
1880 call geom_to_var(nvar,varia)
1881 if (me.eq.king.or..not.out1file) write (iout,*)
1882 & 'Minimizing random structure: Calling MINIMIZE.'
1883 call minimize(etot,varia,iretcode,nfun)
1884 call var_to_geom(nvar,varia)
1886 if (me.eq.king.or..not.out1file)
1887 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1888 if (isnan(etot) .or. etot.gt.1.0d4) then
1889 write (iout,*) "Energy too large",etot,
1890 & " trying another random conformation"
1893 call gen_rand_conf(itmp,*30)
1895 30 write (iout,*) 'Failed to generate random conformation',
1896 & ', itrial=',itrial
1897 write (*,*) 'Processor:',me,
1898 & ' Failed to generate random conformation',
1908 write (iout,'(a,i3,a)') 'Processor:',me,
1909 & ' error in generating random conformation.'
1910 write (*,'(a,i3,a)') 'Processor:',me,
1911 & ' error in generating random conformation.'
1914 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1923 write (iout,'(a,i3,a)') 'Processor:',me,
1924 & ' failed to generate a low-energy random conformation.'
1925 write (*,'(a,i3,a)') 'Processor:',me,
1926 & ' failed to generate a low-energy random conformation.'
1929 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1934 else if (preminim) then
1935 if (start_from_model) then
1939 do while (fail .and. n_model_try.lt.nmodel_start)
1940 write (iout,*) "n_model_try",n_model_try
1942 i_model=iran_num(1,nmodel_start)
1944 if (i_model.eq.list_model_try(k)) exit
1946 if (k.gt.n_model_try) exit
1948 n_model_try=n_model_try+1
1949 list_model_try(n_model_try)=i_model
1950 if (me.eq.king .or. .not. out1file)
1951 & write (iout,*) 'Trying to start from model ',
1952 & pdbfiles_chomo(i_model)(:ilen(pdbfiles_chomo(i_model)))
1955 c(j,i)=chomo(j,i,i_model)
1958 call int_from_cart(.true.,.false.)
1959 call sc_loc_geom(.false.)
1963 dc(j,i)=c(j,i+1)-c(j,i)
1964 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1969 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1970 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1973 if (me.eq.king.or..not.out1file) then
1974 write (iout,*) "Energies before removing overlaps"
1975 call etotal(energia(0))
1976 call enerprint(energia(0))
1978 ! Remove SC overlaps if requested
1980 write (iout,*) 'Calling OVERLAP_SC'
1981 call overlap_sc(fail)
1984 & "Failed to remove overlap from model",i_model
1988 if (me.eq.king.or..not.out1file) then
1989 write (iout,*) "Energies after removing overlaps"
1990 call etotal(energia(0))
1991 call enerprint(energia(0))
1994 ! Search for better SC rotamers if requested
1996 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1997 print *,'SC_move',nft_sc,etot
1998 if (me.eq.king.or..not.out1file)
1999 & write(iout,*) 'SC_move',nft_sc,etot
2001 call etotal(energia(0))
2004 call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),
2005 & 1,MPI_INTEGER,king,CG_COMM,IERROR)
2006 if (n_model_try.gt.nmodel_start .and.
2007 & (me.eq.king .or. out1file)) then
2009 & "All models have irreparable overlaps. Trying randoms starts."
2011 i_model=nmodel_start+1
2015 ! Remove SC overlaps if requested
2017 write (iout,*) 'Calling OVERLAP_SC'
2018 call overlap_sc(fail)
2021 & "Failed to remove overlap"
2024 if (me.eq.king.or..not.out1file) then
2025 write (iout,*) "Energies after removing overlaps"
2026 call etotal(energia(0))
2027 call enerprint(energia(0))
2030 C 8/22/17 AL Minimize initial structure
2032 if (me.eq.king.or..not.out1file) write(iout,*)
2033 & 'Minimizing initial PDB structure: Calling MINIM_DC'
2034 call minim_dc(etot,iretcode,nfun)
2036 call geom_to_var(nvar,varia)
2037 if(me.eq.king.or..not.out1file) write (iout,*)
2038 & 'Minimizing initial PDB structure: Calling MINIMIZE.'
2039 call minimize(etot,varia,iretcode,nfun)
2040 call var_to_geom(nvar,varia)
2042 if (me.eq.king.or..not.out1file)
2043 & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2044 if(me.eq.king.or..not.out1file)
2045 & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2047 if (me.eq.king.or..not.out1file)
2048 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2049 if(me.eq.king.or..not.out1file)
2050 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2054 if (nmodel_start.gt.0 .and. me.eq.king) then
2055 write (iout,'(a)') "Task Starting model"
2057 if (i_start_models(i).gt.nmodel_start) then
2058 write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE"
2060 write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i))
2061 & (:ilen(pdbfiles_chomo(i_start_models(i))))
2066 call chainbuild_cart
2071 kinetic_T=2.0d0/(dimen3*Rb)*EK
2072 write (iout,*) "Initial kinetic energy",EK," kinetic T",kinetic_T
2073 if(me.eq.king.or..not.out1file)then
2083 call etotal(potEcomp)
2084 if (large) call enerprint(potEcomp)
2085 if (RESPA) call etotal_long(potEcomp_long)
2088 t_etotal=t_etotal+MPI_Wtime()-tt0
2090 t_etotal=t_etotal+tcpu()-tt0
2094 c write (iout,*) "PotE-homology",potE-potEcomp(27)
2098 if (amax*d_time .gt. dvmax .and. .not. respa) then
2099 d_time=d_time*dvmax/amax
2100 if(me.eq.king.or..not.out1file) write (iout,*)
2101 & "Time step reduced to",d_time,
2102 & " because of too large initial acceleration."
2104 if(me.eq.king.or..not.out1file)then
2105 write(iout,*) "Potential energy and its components"
2106 call enerprint(potEcomp)
2107 c write(iout,*) (potEcomp(i),i=0,n_ene)
2109 potE=potEcomp(0)-potEcomp(27)
2110 c write (iout,*) "PotE-homology",potE
2114 if (ntwe.ne.0) call statout(itime)
2115 if(me.eq.king.or..not.out1file)
2116 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
2117 & " Kinetic energy",EK," Potential energy",potE,
2118 & " Total energy",totE," Maximum acceleration ",
2121 write (iout,*) "Initial coordinates"
2123 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
2124 & (c(j,i+nres),j=1,3)
2126 write (iout,*) "Initial dC"
2128 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
2129 & (dc(j,i+nres),j=1,3)
2131 write (iout,*) "Initial velocities"
2132 write (iout,"(13x,' backbone ',23x,' side chain')")
2134 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
2135 & (d_t(j,i+nres),j=1,3)
2137 write (iout,*) "Initial accelerations"
2139 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2140 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
2141 & (d_a(j,i+nres),j=1,3)
2147 d_t_old(j,i)=d_t(j,i)
2148 d_a_old(j,i)=d_a(j,i)
2150 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2159 call etotal_short(energia_short)
2160 if (large) call enerprint(potEcomp)
2163 t_eshort=t_eshort+MPI_Wtime()-tt0
2165 t_eshort=t_eshort+tcpu()-tt0
2170 if(.not.out1file .and. large) then
2171 write (iout,*) "energia_long",energia_long(0),
2172 & " energia_short",energia_short(0),
2173 & " total",energia_long(0)+energia_short(0)
2174 write (iout,*) "Initial fast-force accelerations"
2176 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2177 & (d_a(j,i+nres),j=1,3)
2180 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2183 d_a_short(j,i)=d_a(j,i)
2192 call etotal_long(energia_long)
2193 if (large) call enerprint(potEcomp)
2196 t_elong=t_elong+MPI_Wtime()-tt0
2198 t_elong=t_elong+tcpu()-tt0
2203 if(.not.out1file .and. large) then
2204 write (iout,*) "energia_long",energia_long(0)
2205 write (iout,*) "Initial slow-force accelerations"
2207 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2208 & (d_a(j,i+nres),j=1,3)
2212 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2214 t_enegrad=t_enegrad+tcpu()-tt0
2219 c-----------------------------------------------------------
2220 subroutine random_vel
2222 include 'DIMENSIONS'
2223 include 'COMMON.CONTROL'
2224 include 'COMMON.VAR'
2227 include 'COMMON.LAGRANGE.5diag'
2229 include 'COMMON.LAGRANGE'
2232 include 'COMMON.LANGEVIN'
2235 include 'COMMON.LANGEVIN.lang0.5diag'
2237 include 'COMMON.LANGEVIN.lang0'
2240 include 'COMMON.CHAIN'
2241 include 'COMMON.DERIV'
2242 include 'COMMON.GEO'
2243 include 'COMMON.LOCAL'
2244 include 'COMMON.INTERACT'
2245 include 'COMMON.IOUNITS'
2246 include 'COMMON.NAMES'
2247 include 'COMMON.TIME1'
2248 double precision xv,sigv,lowb,highb,vec_afm(3),Ek1,Ek2,Ek3,aux
2249 integer i,ii,j,k,l,ind
2250 double precision anorm_distr
2251 logical lprn /.false./
2253 integer ichain,n,innt,inct,ibeg,ierr
2254 double precision work(8*maxres6)
2255 integer iwork(maxres6)
2256 double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),
2257 & Gvec(maxres2_chain,maxres2_chain)
2258 common /przechowalnia/Ghalf,Geigen,Gvec
2260 double precision inertia(maxres2_chain,maxres2_chain)
2262 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2263 c First generate velocities in the eigenspace of the G matrix
2264 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2267 write (iout,*) "Random_vel, fivediag"
2276 n=dimen_chain(ichain)
2277 innt=iposd_chain(ichain)
2280 write (iout,*) "Chain",ichain," n",n," start",innt
2282 if (i.lt.inct-1) then
2283 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),
2285 else if (i.eq.inct-1) then
2286 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
2288 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
2292 ghalf(ind+1)=dmorig(innt)
2293 ghalf(ind+2)=du1orig(innt)
2294 ghalf(ind+3)=dmorig(innt+1)
2298 c write (iout,*) "i",i," ind",ind," indu2",innt+i-2,
2299 c & " indu1",innt+i-1," indm",innt+i
2300 ghalf(ind+1)=du2orig(innt-1+i-2)
2301 ghalf(ind+2)=du1orig(innt-1+i-1)
2302 ghalf(ind+3)=dmorig(innt-1+i)
2303 c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
2304 c & "DU1",innt-1+i-1,"DM ",innt-1+i
2312 inertia(i,j)=ghalf(ind)
2313 inertia(j,i)=ghalf(ind)
2318 write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
2319 write (iout,*) "Five-diagonal inertia matrix, lower triangle"
2320 call matoutr(n,ghalf)
2322 call gldiag(maxres2_chain,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
2324 write (iout,'(//a,i3)')
2325 & "Eigenvectors and eigenvalues of the G matrix chain",ichain
2326 call eigout(n,n,maxres2_chain,maxres2_chain,Gvec,Geigen)
2329 c check diagonalization
2335 aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
2339 write (iout,*) i,j,aux,geigen(i)
2341 write (iout,*) i,j,aux
2351 sigv=dsqrt((Rb*t_bath)/geigen(i))
2354 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2355 EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
2356 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2357 c & " d_t_work_new",d_t_work_new(ii)
2365 d_t_work(ind)=d_t_work(ind)
2366 & +Gvec(i,j)*d_t_work_new((j-1)*3+k)
2368 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2377 aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
2383 c Transfer to the d_t vector
2384 innt=chain_border(1,ichain)
2385 inct=chain_border(2,ichain)
2387 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
2391 d_t(j,i)=d_t_work(ind)
2393 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2396 d_t(j,i+nres)=d_t_work(ind)
2403 write (iout,*) "Random velocities in the Calpha,SC space"
2405 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2406 & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2409 call kinetic_CASC(Ek1)
2411 ! Transform the velocities to virtual-bond space
2419 if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
2421 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2425 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2426 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2439 c write (iout,*) "Shifting accelerations"
2441 c write (iout,*) "ichain",chain_border1(1,ichain)-1,
2442 c & chain_border1(1,ichain)
2443 d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
2444 d_t(:,chain_border1(1,ichain))=0.0d0
2446 c write (iout,*) "Adding accelerations"
2448 c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2449 c & chain_border(2,ichain-1)
2450 d_t(:,chain_border1(1,ichain)-1)=
2451 & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2452 d_t(:,chain_border(2,ichain-1))=0.0d0
2455 write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2456 & chain_border(2,ichain-1)
2457 d_t(:,chain_border1(1,ichain)-1)=
2458 & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2459 d_t(:,chain_border(2,ichain-1))=0.0d0
2464 c d_t(j,0)=d_t(j,nnt)
2467 innt=chain_border(1,ichain)
2468 inct=chain_border(2,ichain)
2469 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
2470 c write (iout,*) "ibeg",ibeg
2472 d_t(j,ibeg)=d_t(j,innt)
2476 if (iabs(itype(i).eq.10)) then
2477 c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
2479 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2483 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2484 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2493 & "Random velocities in the virtual-bond-vector space"
2494 write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
2496 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2497 & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2500 write (iout,*) "Kinetic energy from inertia matrix eigenvalues",
2503 & "Kinetic temperatures from inertia matrix eigenvalues",
2506 write (iout,*) "Kinetic energy from inertia matrix",Ek3
2507 write (iout,*) "Kinetic temperatures from inertia",
2508 & 2*Ek3/(3*dimen*Rb)
2510 write (iout,*) "Kinetic energy from velocities in CA-SC space",
2513 & "Kinetic temperatures from velovities in CA-SC space",
2514 & 2*Ek1/(3*dimen*Rb)
2517 & "Kinetic energy from virtual-bond-vector velocities",Ek1
2519 & "Kinetic temperature from virtual-bond-vector velocities ",
2528 sigv=dsqrt((Rb*t_bath)/geigen(i))
2531 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2533 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2534 c & " d_t_work_new",d_t_work_new(ii)
2537 C if (SELFGUIDE.gt.0) then
2540 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
2541 C distance=distance+vec_afm(j)**2
2543 C distance=dsqrt(distance)
2545 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
2546 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
2547 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
2548 C & d_t_work_new(j+(afmend-1)*3)
2559 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2562 c write (iout,*) "Ek from eigenvectors",Ek1
2564 c Transform velocities to UNRES coordinate space
2570 d_t_work(ind)=d_t_work(ind)
2571 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2573 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2577 c Transfer to the d_t vector
2579 d_t(j,0)=d_t_work(j)
2585 d_t(j,i)=d_t_work(ind)
2589 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2592 d_t(j,i+nres)=d_t_work(ind)
2597 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2598 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2604 c-----------------------------------------------------------
2605 subroutine sd_verlet_p_setup
2606 c Sets up the parameters of stochastic Verlet algorithm
2608 include 'DIMENSIONS'
2612 include 'COMMON.CONTROL'
2613 include 'COMMON.VAR'
2616 include 'COMMON.LANGEVIN'
2619 include 'COMMON.LANGEVIN.lang0.5diag'
2621 include 'COMMON.LANGEVIN.lang0'
2624 include 'COMMON.CHAIN'
2625 include 'COMMON.DERIV'
2626 include 'COMMON.GEO'
2627 include 'COMMON.LOCAL'
2628 include 'COMMON.INTERACT'
2629 include 'COMMON.IOUNITS'
2630 include 'COMMON.NAMES'
2631 include 'COMMON.TIME1'
2632 double precision emgdt(MAXRES6),
2633 & pterm,vterm,rho,rhoc,vsig,
2634 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2635 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2636 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2637 logical lprn /.false./
2638 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2639 double precision ktm
2646 c AL 8/17/04 Code adapted from tinker
2648 c Get the frictional and random terms for stochastic dynamics in the
2649 c eigenspace of mass-scaled UNRES friction matrix
2652 gdt = fricgam(i) * d_time
2654 c Stochastic dynamics reduces to simple MD for zero friction
2656 if (gdt .le. zero) then
2657 pfric_vec(i) = 1.0d0
2658 vfric_vec(i) = d_time
2659 afric_vec(i) = 0.5d0 * d_time * d_time
2660 prand_vec(i) = 0.0d0
2661 vrand_vec1(i) = 0.0d0
2662 vrand_vec2(i) = 0.0d0
2664 c Analytical expressions when friction coefficient is large
2667 if (gdt .ge. gdt_radius) then
2670 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2671 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2672 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2673 vterm = 1.0d0 - egdt**2
2674 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2676 c Use series expansions when friction coefficient is small
2687 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2688 & - gdt5/120.0d0 + gdt6/720.0d0
2689 & - gdt7/5040.0d0 + gdt8/40320.0d0
2690 & - gdt9/362880.0d0) / fricgam(i)**2
2691 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2692 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2693 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2694 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2695 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2696 & + 127.0d0*gdt9/90720.0d0
2697 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2698 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2699 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2700 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2701 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2702 & - 17.0d0*gdt2/1280.0d0
2703 & + 17.0d0*gdt3/6144.0d0
2704 & + 40967.0d0*gdt4/34406400.0d0
2705 & - 57203.0d0*gdt5/275251200.0d0
2706 & - 1429487.0d0*gdt6/13212057600.0d0)
2709 c Compute the scaling factors of random terms for the nonzero friction case
2711 ktm = 0.5d0*d_time/fricgam(i)
2712 psig = dsqrt(ktm*pterm) / fricgam(i)
2713 vsig = dsqrt(ktm*vterm)
2714 rhoc = dsqrt(1.0d0 - rho*rho)
2716 vrand_vec1(i) = vsig * rho
2717 vrand_vec2(i) = vsig * rhoc
2722 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2725 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2726 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2730 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2732 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2733 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2734 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2735 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2736 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2737 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2739 t_sdsetup=t_sdsetup+MPI_Wtime()
2741 t_sdsetup=t_sdsetup+tcpu()-tt0
2745 c-------------------------------------------------------------
2746 subroutine eigtransf1(n,ndim,ab,d,c)
2749 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2755 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2761 c-------------------------------------------------------------
2762 subroutine eigtransf(n,ndim,a,b,d,c)
2765 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2771 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2777 c-------------------------------------------------------------
2778 subroutine sd_verlet1
2779 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2781 include 'DIMENSIONS'
2782 include 'COMMON.CONTROL'
2783 include 'COMMON.VAR'
2786 include 'COMMON.LAGRANGE.5diag'
2788 include 'COMMON.LAGRANGE'
2791 include 'COMMON.LANGEVIN'
2794 include 'COMMON.LANGEVIN.lang0.5diag'
2796 include 'COMMON.LANGEVIN.lang0'
2799 include 'COMMON.CHAIN'
2800 include 'COMMON.DERIV'
2801 include 'COMMON.GEO'
2802 include 'COMMON.LOCAL'
2803 include 'COMMON.INTERACT'
2804 include 'COMMON.IOUNITS'
2805 include 'COMMON.NAMES'
2806 double precision stochforcvec(MAXRES6)
2807 common /stochcalc/ stochforcvec
2808 logical lprn /.false./
2809 integer i,j,ind,inres
2811 c write (iout,*) "dc_old"
2813 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2814 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2817 dc_work(j)=dc_old(j,0)
2818 d_t_work(j)=d_t_old(j,0)
2819 d_a_work(j)=d_a_old(j,0)
2824 dc_work(ind+j)=dc_old(j,i)
2825 d_t_work(ind+j)=d_t_old(j,i)
2826 d_a_work(ind+j)=d_a_old(j,i)
2831 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2833 dc_work(ind+j)=dc_old(j,i+nres)
2834 d_t_work(ind+j)=d_t_old(j,i+nres)
2835 d_a_work(ind+j)=d_a_old(j,i+nres)
2843 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2847 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2848 & vfric_mat(i,j),afric_mat(i,j),
2849 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2857 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2858 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2859 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2860 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2862 d_t_work_new(i)=ddt1+0.5d0*ddt2
2863 d_t_work(i)=ddt1+ddt2
2868 d_t(j,0)=d_t_work(j)
2873 dc(j,i)=dc_work(ind+j)
2874 d_t(j,i)=d_t_work(ind+j)
2879 if (itype(i).ne.10) then
2882 dc(j,inres)=dc_work(ind+j)
2883 d_t(j,inres)=d_t_work(ind+j)
2890 c--------------------------------------------------------------------------
2891 subroutine sd_verlet2
2892 c Calculating the adjusted velocities for accelerations
2894 include 'DIMENSIONS'
2895 include 'COMMON.CONTROL'
2896 include 'COMMON.VAR'
2899 include 'COMMON.LAGRANGE.5diag'
2901 include 'COMMON.LAGRANGE'
2904 include 'COMMON.LANGEVIN'
2907 include 'COMMON.LANGEVIN.lang0.5diag'
2909 include 'COMMON.LANGEVIN.lang0'
2912 include 'COMMON.CHAIN'
2913 include 'COMMON.DERIV'
2914 include 'COMMON.GEO'
2915 include 'COMMON.LOCAL'
2916 include 'COMMON.INTERACT'
2917 include 'COMMON.IOUNITS'
2918 include 'COMMON.NAMES'
2919 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2920 common /stochcalc/ stochforcvec
2921 integer i,j,ind,inres
2923 c Compute the stochastic forces which contribute to velocity change
2925 call stochastic_force(stochforcvecV)
2932 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2933 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2934 & vrand_mat2(i,j)*stochforcvecV(j)
2936 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2940 d_t(j,0)=d_t_work(j)
2945 d_t(j,i)=d_t_work(ind+j)
2950 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2953 d_t(j,inres)=d_t_work(ind+j)
2960 c-----------------------------------------------------------
2961 subroutine sd_verlet_ciccotti_setup
2962 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2965 include 'DIMENSIONS'
2969 include 'COMMON.CONTROL'
2970 include 'COMMON.VAR'
2973 include 'COMMON.LAGRANGE.5diag'
2975 include 'COMMON.LAGRANGE'
2977 include 'COMMON.LANGEVIN'
2978 include 'COMMON.CHAIN'
2979 include 'COMMON.DERIV'
2980 include 'COMMON.GEO'
2981 include 'COMMON.LOCAL'
2982 include 'COMMON.INTERACT'
2983 include 'COMMON.IOUNITS'
2984 include 'COMMON.NAMES'
2985 include 'COMMON.TIME1'
2986 double precision emgdt(MAXRES6),
2987 & pterm,vterm,rho,rhoc,vsig,
2988 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2989 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2990 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2991 logical lprn /.false./
2992 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2993 double precision ktm
3001 c AL 8/17/04 Code adapted from tinker
3003 c Get the frictional and random terms for stochastic dynamics in the
3004 c eigenspace of mass-scaled UNRES friction matrix
3007 write (iout,*) "i",i," fricgam",fricgam(i)
3008 gdt = fricgam(i) * d_time
3010 c Stochastic dynamics reduces to simple MD for zero friction
3012 if (gdt .le. zero) then
3013 pfric_vec(i) = 1.0d0
3014 vfric_vec(i) = d_time
3015 afric_vec(i) = 0.5d0*d_time*d_time
3016 prand_vec(i) = afric_vec(i)
3017 vrand_vec2(i) = vfric_vec(i)
3019 c Analytical expressions when friction coefficient is large
3024 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3025 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3026 prand_vec(i) = afric_vec(i)
3027 vrand_vec2(i) = vfric_vec(i)
3029 c Compute the scaling factors of random terms for the nonzero friction case
3031 c ktm = 0.5d0*d_time/fricgam(i)
3032 c psig = dsqrt(ktm*pterm) / fricgam(i)
3033 c vsig = dsqrt(ktm*vterm)
3034 c prand_vec(i) = psig*afric_vec(i)
3035 c vrand_vec2(i) = vsig*vfric_vec(i)
3040 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
3043 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
3044 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3048 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3050 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3051 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3052 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3053 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3054 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3056 t_sdsetup=t_sdsetup+MPI_Wtime()
3058 t_sdsetup=t_sdsetup+tcpu()-tt0
3062 c-------------------------------------------------------------
3063 subroutine sd_verlet1_ciccotti
3064 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
3066 include 'DIMENSIONS'
3070 include 'COMMON.CONTROL'
3071 include 'COMMON.VAR'
3074 include 'COMMON.LAGRANGE.5diag'
3076 include 'COMMON.LAGRANGE'
3078 include 'COMMON.LANGEVIN'
3079 include 'COMMON.CHAIN'
3080 include 'COMMON.DERIV'
3081 include 'COMMON.GEO'
3082 include 'COMMON.LOCAL'
3083 include 'COMMON.INTERACT'
3084 include 'COMMON.IOUNITS'
3085 include 'COMMON.NAMES'
3086 double precision stochforcvec(MAXRES6)
3087 common /stochcalc/ stochforcvec
3088 logical lprn /.false./
3091 c write (iout,*) "dc_old"
3093 c write (iout,'(i5,3f10.5,5x,3f10.5)')
3094 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3097 dc_work(j)=dc_old(j,0)
3098 d_t_work(j)=d_t_old(j,0)
3099 d_a_work(j)=d_a_old(j,0)
3104 dc_work(ind+j)=dc_old(j,i)
3105 d_t_work(ind+j)=d_t_old(j,i)
3106 d_a_work(ind+j)=d_a_old(j,i)
3111 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3113 dc_work(ind+j)=dc_old(j,i+nres)
3114 d_t_work(ind+j)=d_t_old(j,i+nres)
3115 d_a_work(ind+j)=d_a_old(j,i+nres)
3124 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
3128 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
3129 & vfric_mat(i,j),afric_mat(i,j),
3130 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3138 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
3139 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3140 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3141 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3143 d_t_work_new(i)=ddt1+0.5d0*ddt2
3144 d_t_work(i)=ddt1+ddt2
3149 d_t(j,0)=d_t_work(j)
3154 dc(j,i)=dc_work(ind+j)
3155 d_t(j,i)=d_t_work(ind+j)
3160 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3163 dc(j,inres)=dc_work(ind+j)
3164 d_t(j,inres)=d_t_work(ind+j)
3171 c--------------------------------------------------------------------------
3172 subroutine sd_verlet2_ciccotti
3173 c Calculating the adjusted velocities for accelerations
3175 include 'DIMENSIONS'
3176 include 'COMMON.CONTROL'
3177 include 'COMMON.VAR'
3180 include 'COMMON.LAGRANGE.5diag'
3182 include 'COMMON.LAGRANGE'
3184 include 'COMMON.LANGEVIN'
3185 include 'COMMON.CHAIN'
3186 include 'COMMON.DERIV'
3187 include 'COMMON.GEO'
3188 include 'COMMON.LOCAL'
3189 include 'COMMON.INTERACT'
3190 include 'COMMON.IOUNITS'
3191 include 'COMMON.NAMES'
3192 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
3193 common /stochcalc/ stochforcvec
3196 c Compute the stochastic forces which contribute to velocity change
3198 call stochastic_force(stochforcvecV)
3205 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3206 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3207 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3209 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3213 d_t(j,0)=d_t_work(j)
3218 d_t(j,i)=d_t_work(ind+j)
3223 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3226 d_t(j,inres)=d_t_work(ind+j)