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) call statout(itime)
212 if (itype(i).ne.10 .and. itype(i).ne.21) then
215 v_work(ind)=d_t(j,i+nres)
220 write (66,'(80f10.5)')
221 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
225 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
227 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
229 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
232 if (mod(itime,ntwx).eq.0) then
233 write (tytul,'("time",f8.2)') totT
235 call hairpin(.true.,nharp,iharp)
236 call secondary2(.true.)
237 call pdbout(potE,tytul,ipdb)
242 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
243 open(irest2,file=rest2name,status='unknown')
244 write(irest2,*) totT,EK,potE,totE,t_bath
246 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
249 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
260 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
262 & 'MD calculations setup:',t_MDsetup,
263 & 'Energy & gradient evaluation:',t_enegrad,
264 & 'Stochastic MD setup:',t_langsetup,
265 & 'Stochastic MD step setup:',t_sdsetup,
267 write (iout,'(/28(1h=),a25,27(1h=))')
268 & ' End of MD calculation '
270 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
272 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
273 & " time_fricmatmult",time_fricmatmult," time_fsample ",
278 c-------------------------------------------------------------------------------
279 subroutine velverlet_step(itime)
280 c-------------------------------------------------------------------------------
281 c Perform a single velocity Verlet step; the time step can be rescaled if
282 c increments in accelerations exceed the threshold
283 c-------------------------------------------------------------------------------
284 implicit real*8 (a-h,o-z)
288 integer ierror,ierrcode
290 include 'COMMON.SETUP'
291 include 'COMMON.CONTROL'
295 include 'COMMON.LANGEVIN'
297 include 'COMMON.LANGEVIN.lang0'
299 include 'COMMON.CHAIN'
300 include 'COMMON.DERIV'
302 include 'COMMON.LOCAL'
303 include 'COMMON.INTERACT'
304 include 'COMMON.IOUNITS'
305 include 'COMMON.NAMES'
306 include 'COMMON.TIME1'
307 include 'COMMON.MUCA'
308 double precision vcm(3),incr(3)
309 double precision cm(3),L(3)
310 integer ilen,count,rstcount
313 integer maxcount_scale /20/
315 double precision stochforcvec(MAXRES6)
316 common /stochcalc/ stochforcvec
324 else if (lang.eq.2 .or. lang.eq.3) then
326 call stochastic_force(stochforcvec)
329 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
331 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
338 icount_scale=icount_scale+1
339 if (icount_scale.gt.maxcount_scale) then
341 & "ERROR: too many attempts at scaling down the time step. ",
342 & "amax=",amax,"epdrift=",epdrift,
343 & "damax=",damax,"edriftmax=",edriftmax,
347 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
351 c First step of the velocity Verlet algorithm
356 else if (lang.eq.3) then
358 call sd_verlet1_ciccotti
360 else if (lang.eq.1) then
365 c Build the chain from the newly calculated coordinates
367 if (rattle) call rattle1
369 if (large.and. mod(itime,ntwe).eq.0) then
370 write (iout,*) "Cartesian and internal coordinates: step 1"
375 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
376 & (dc(j,i+nres),j=1,3)
378 write (iout,*) "Accelerations"
380 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
381 & (d_a(j,i+nres),j=1,3)
383 write (iout,*) "Velocities, step 1"
385 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
386 & (d_t(j,i+nres),j=1,3)
395 c Calculate energy and forces
397 call etotal(potEcomp)
398 if (large.and. mod(itime,ntwe).eq.0)
399 & call enerprint(potEcomp)
402 t_etotal=t_etotal+MPI_Wtime()-tt0
404 t_etotal=t_etotal+tcpu()-tt0
407 potE=potEcomp(0)-potEcomp(20)
409 c Get the new accelerations
412 t_enegrad=t_enegrad+MPI_Wtime()-tt0
414 t_enegrad=t_enegrad+tcpu()-tt0
416 c Determine maximum acceleration and scale down the timestep if needed
418 amax=amax/(itime_scal+1)**2
419 call predict_edrift(epdrift)
420 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
421 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
423 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
425 itime_scal=itime_scal+ifac_time
426 c fac_time=dmin1(damax/amax,0.5d0)
427 fac_time=0.5d0**ifac_time
428 d_time=d_time*fac_time
429 if (lang.eq.2 .or. lang.eq.3) then
431 c write (iout,*) "Calling sd_verlet_setup: 1"
432 c Rescale the stochastic forces and recalculate or restore
433 c the matrices of tinker integrator
434 if (itime_scal.gt.maxflag_stoch) then
435 if (large) write (iout,'(a,i5,a)')
436 & "Calculate matrices for stochastic step;",
437 & " itime_scal ",itime_scal
439 call sd_verlet_p_setup
441 call sd_verlet_ciccotti_setup
443 write (iout,'(2a,i3,a,i3,1h.)')
444 & "Warning: cannot store matrices for stochastic",
445 & " integration because the index",itime_scal,
446 & " is greater than",maxflag_stoch
447 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
448 & " integration Langevin algorithm for better efficiency."
449 else if (flag_stoch(itime_scal)) then
450 if (large) write (iout,'(a,i5,a,l1)')
451 & "Restore matrices for stochastic step; itime_scal ",
452 & itime_scal," flag ",flag_stoch(itime_scal)
455 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
456 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
457 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
458 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
459 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
460 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
464 if (large) write (iout,'(2a,i5,a,l1)')
465 & "Calculate & store matrices for stochastic step;",
466 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
468 call sd_verlet_p_setup
470 call sd_verlet_ciccotti_setup
472 flag_stoch(ifac_time)=.true.
475 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
476 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
477 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
478 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
479 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
480 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
484 fac_time=1.0d0/dsqrt(fac_time)
486 stochforcvec(i)=fac_time*stochforcvec(i)
489 else if (lang.eq.1) then
490 c Rescale the accelerations due to stochastic forces
491 fac_time=1.0d0/dsqrt(fac_time)
493 d_as_work(i)=d_as_work(i)*fac_time
496 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
497 & "itime",itime," Timestep scaled down to ",
498 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
500 c Second step of the velocity Verlet algorithm
505 else if (lang.eq.3) then
507 call sd_verlet2_ciccotti
509 else if (lang.eq.1) then
514 if (rattle) call rattle2
516 if (d_time.ne.d_time0) then
519 if (lang.eq.2 .or. lang.eq.3) then
520 if (large) write (iout,'(a)')
521 & "Restore original matrices for stochastic step"
522 c write (iout,*) "Calling sd_verlet_setup: 2"
523 c Restore the matrices of tinker integrator if the time step has been restored
526 pfric_mat(i,j)=pfric0_mat(i,j,0)
527 afric_mat(i,j)=afric0_mat(i,j,0)
528 vfric_mat(i,j)=vfric0_mat(i,j,0)
529 prand_mat(i,j)=prand0_mat(i,j,0)
530 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
531 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
540 c Calculate the kinetic and the total energy and the kinetic temperature
545 c write (iout,*) "step",itime," EK",EK," EK1",EK1
547 c Couple the system to Berendsen bath if needed
548 if (tbf .and. lang.eq.0) then
551 kinetic_T=2.0d0/(dimen3*Rb)*EK
552 c Backup the coordinates, velocities, and accelerations
556 d_t_old(j,i)=d_t(j,i)
557 d_a_old(j,i)=d_a(j,i)
561 if (mod(itime,ntwe).eq.0 .and. large) then
562 write (iout,*) "Velocities, step 2"
564 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
565 & (d_t(j,i+nres),j=1,3)
571 c-------------------------------------------------------------------------------
572 subroutine RESPA_step(itime)
573 c-------------------------------------------------------------------------------
574 c Perform a single RESPA step.
575 c-------------------------------------------------------------------------------
576 implicit real*8 (a-h,o-z)
580 integer IERROR,ERRCODE
582 include 'COMMON.SETUP'
583 include 'COMMON.CONTROL'
587 include 'COMMON.LANGEVIN'
589 include 'COMMON.LANGEVIN.lang0'
591 include 'COMMON.CHAIN'
592 include 'COMMON.DERIV'
594 include 'COMMON.LOCAL'
595 include 'COMMON.INTERACT'
596 include 'COMMON.IOUNITS'
597 include 'COMMON.NAMES'
598 include 'COMMON.TIME1'
599 double precision energia_short(0:n_ene),
600 & energia_long(0:n_ene)
601 double precision cm(3),L(3),vcm(3),incr(3)
602 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
603 & d_a_old0(3,0:maxres2)
604 logical PRINT_AMTS_MSG /.false./
605 integer ilen,count,rstcount
608 integer maxcount_scale /10/
610 double precision stochforcvec(MAXRES6)
611 common /stochcalc/ stochforcvec
614 common /cipiszcze/ itt
617 if (large.and. mod(itime,ntwe).eq.0) then
618 write (iout,*) "***************** RESPA itime",itime
619 write (iout,*) "Cartesian and internal coordinates: step 0"
621 call pdbout(0.0d0,"cipiszcze",iout)
623 write (iout,*) "Accelerations from long-range forces"
625 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
626 & (d_a(j,i+nres),j=1,3)
628 write (iout,*) "Velocities, step 0"
630 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
631 & (d_t(j,i+nres),j=1,3)
636 c Perform the initial RESPA step (increment velocities)
637 c write (iout,*) "*********************** RESPA ini"
640 if (mod(itime,ntwe).eq.0 .and. large) then
641 write (iout,*) "Velocities, end"
643 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
644 & (d_t(j,i+nres),j=1,3)
648 c Compute the short-range forces
654 C 7/2/2009 commented out
656 c call etotal_short(energia_short)
659 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
662 d_a(j,i)=d_a_short(j,i)
666 if (large.and. mod(itime,ntwe).eq.0) then
667 write (iout,*) "energia_short",energia_short(0)
668 write (iout,*) "Accelerations from short-range forces"
670 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
671 & (d_a(j,i+nres),j=1,3)
676 t_enegrad=t_enegrad+MPI_Wtime()-tt0
678 t_enegrad=t_enegrad+tcpu()-tt0
683 d_t_old(j,i)=d_t(j,i)
684 d_a_old(j,i)=d_a(j,i)
687 c 6/30/08 A-MTS: attempt at increasing the split number
690 dc_old0(j,i)=dc_old(j,i)
691 d_t_old0(j,i)=d_t_old(j,i)
692 d_a_old0(j,i)=d_a_old(j,i)
695 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
696 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
703 c write (iout,*) "itime",itime," ntime_split",ntime_split
704 c Split the time step
705 d_time=d_time0/ntime_split
706 c Perform the short-range RESPA steps (velocity Verlet increments of
707 c positions and velocities using short-range forces)
708 c write (iout,*) "*********************** RESPA split"
709 do itsplit=1,ntime_split
712 else if (lang.eq.2 .or. lang.eq.3) then
714 call stochastic_force(stochforcvec)
717 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
719 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
724 c First step of the velocity Verlet algorithm
729 else if (lang.eq.3) then
731 call sd_verlet1_ciccotti
733 else if (lang.eq.1) then
738 c Build the chain from the newly calculated coordinates
740 if (rattle) call rattle1
742 if (large.and. mod(itime,ntwe).eq.0) then
743 write (iout,*) "***** ITSPLIT",itsplit
744 write (iout,*) "Cartesian and internal coordinates: step 1"
745 call pdbout(0.0d0,"cipiszcze",iout)
748 write (iout,*) "Velocities, step 1"
750 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
751 & (d_t(j,i+nres),j=1,3)
760 c Calculate energy and forces
762 call etotal_short(energia_short)
763 if (large.and. mod(itime,ntwe).eq.0)
764 & call enerprint(energia_short)
767 t_eshort=t_eshort+MPI_Wtime()-tt0
769 t_eshort=t_eshort+tcpu()-tt0
773 c Get the new accelerations
775 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
778 d_a_short(j,i)=d_a(j,i)
782 if (large.and. mod(itime,ntwe).eq.0) then
783 write (iout,*)"energia_short",energia_short(0)
784 write (iout,*) "Accelerations from short-range forces"
786 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
787 & (d_a(j,i+nres),j=1,3)
792 c Determine maximum acceleration and scale down the timestep if needed
794 amax=amax/ntime_split**2
795 call predict_edrift(epdrift)
796 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
797 & write (iout,*) "amax",amax," damax",damax,
798 & " epdrift",epdrift," epdriftmax",epdriftmax
799 c Exit loop and try with increased split number if the change of
800 c acceleration is too big
801 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
802 if (ntime_split.lt.maxtime_split) then
804 ntime_split=ntime_split*2
807 dc_old(j,i)=dc_old0(j,i)
808 d_t_old(j,i)=d_t_old0(j,i)
809 d_a_old(j,i)=d_a_old0(j,i)
812 if (PRINT_AMTS_MSG) then
813 write (iout,*) "acceleration/energy drift too large",amax,
814 & epdrift," split increased to ",ntime_split," itime",itime,
820 & "Uh-hu. Bumpy landscape. Maximum splitting number",
822 & " already reached!!! Trying to carry on!"
826 t_enegrad=t_enegrad+MPI_Wtime()-tt0
828 t_enegrad=t_enegrad+tcpu()-tt0
830 c Second step of the velocity Verlet algorithm
835 else if (lang.eq.3) then
837 call sd_verlet2_ciccotti
839 else if (lang.eq.1) then
844 if (rattle) call rattle2
845 c Backup the coordinates, velocities, and accelerations
849 d_t_old(j,i)=d_t(j,i)
850 d_a_old(j,i)=d_a(j,i)
857 c Restore the time step
859 c Compute long-range forces
866 call etotal_long(energia_long)
867 if (large.and. mod(itime,ntwe).eq.0)
868 & call enerprint(energia_long)
871 t_elong=t_elong+MPI_Wtime()-tt0
873 t_elong=t_elong+tcpu()-tt0
879 t_enegrad=t_enegrad+MPI_Wtime()-tt0
881 t_enegrad=t_enegrad+tcpu()-tt0
883 c Compute accelerations from long-range forces
885 if (large.and. mod(itime,ntwe).eq.0) then
886 write (iout,*) "energia_long",energia_long(0)
887 write (iout,*) "Cartesian and internal coordinates: step 2"
889 call pdbout(0.0d0,"cipiszcze",iout)
891 write (iout,*) "Accelerations from long-range forces"
893 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
894 & (d_a(j,i+nres),j=1,3)
896 write (iout,*) "Velocities, step 2"
898 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
899 & (d_t(j,i+nres),j=1,3)
903 c Compute the final RESPA step (increment velocities)
904 c write (iout,*) "*********************** RESPA fin"
906 c Compute the complete potential energy
908 potEcomp(i)=energia_short(i)+energia_long(i)
910 potE=potEcomp(0)-potEcomp(20)
911 c potE=energia_short(0)+energia_long(0)
913 c Calculate the kinetic and the total energy and the kinetic temperature
916 c Couple the system to Berendsen bath if needed
917 if (tbf .and. lang.eq.0) then
920 kinetic_T=2.0d0/(dimen3*Rb)*EK
921 c Backup the coordinates, velocities, and accelerations
923 if (mod(itime,ntwe).eq.0 .and. large) then
924 write (iout,*) "Velocities, end"
926 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
927 & (d_t(j,i+nres),j=1,3)
933 c---------------------------------------------------------------------
935 c First and last RESPA step (incrementing velocities using long-range
937 implicit real*8 (a-h,o-z)
939 include 'COMMON.CONTROL'
942 include 'COMMON.CHAIN'
943 include 'COMMON.DERIV'
945 include 'COMMON.LOCAL'
946 include 'COMMON.INTERACT'
947 include 'COMMON.IOUNITS'
948 include 'COMMON.NAMES'
950 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
954 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
958 if (itype(i).ne.10 .and. itype(i).ne.21) then
961 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
967 c-----------------------------------------------------------------
969 c Applying velocity Verlet algorithm - step 1 to coordinates
970 implicit real*8 (a-h,o-z)
972 include 'COMMON.CONTROL'
975 include 'COMMON.CHAIN'
976 include 'COMMON.DERIV'
978 include 'COMMON.LOCAL'
979 include 'COMMON.INTERACT'
980 include 'COMMON.IOUNITS'
981 include 'COMMON.NAMES'
982 double precision adt,adt2
985 write (iout,*) "VELVERLET1 START: DC"
987 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
988 & (dc(j,i+nres),j=1,3)
992 adt=d_a_old(j,0)*d_time
994 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
995 d_t_new(j,0)=d_t_old(j,0)+adt2
996 d_t(j,0)=d_t_old(j,0)+adt
1000 adt=d_a_old(j,i)*d_time
1002 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1003 d_t_new(j,i)=d_t_old(j,i)+adt2
1004 d_t(j,i)=d_t_old(j,i)+adt
1008 if (itype(i).ne.10 .and. itype(i).ne.21) then
1011 adt=d_a_old(j,inres)*d_time
1013 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1014 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1015 d_t(j,inres)=d_t_old(j,inres)+adt
1020 write (iout,*) "VELVERLET1 END: DC"
1022 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1023 & (dc(j,i+nres),j=1,3)
1028 c---------------------------------------------------------------------
1030 c Step 2 of the velocity Verlet algorithm: update velocities
1031 implicit real*8 (a-h,o-z)
1032 include 'DIMENSIONS'
1033 include 'COMMON.CONTROL'
1034 include 'COMMON.VAR'
1036 include 'COMMON.CHAIN'
1037 include 'COMMON.DERIV'
1038 include 'COMMON.GEO'
1039 include 'COMMON.LOCAL'
1040 include 'COMMON.INTERACT'
1041 include 'COMMON.IOUNITS'
1042 include 'COMMON.NAMES'
1044 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1048 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1052 if (itype(i).ne.10 .and. itype(i).ne.21) then
1055 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1061 c-----------------------------------------------------------------
1062 subroutine sddir_precalc
1063 c Applying velocity Verlet algorithm - step 1 to coordinates
1064 implicit real*8 (a-h,o-z)
1065 include 'DIMENSIONS'
1069 include 'COMMON.CONTROL'
1070 include 'COMMON.VAR'
1073 include 'COMMON.LANGEVIN'
1075 include 'COMMON.LANGEVIN.lang0'
1077 include 'COMMON.CHAIN'
1078 include 'COMMON.DERIV'
1079 include 'COMMON.GEO'
1080 include 'COMMON.LOCAL'
1081 include 'COMMON.INTERACT'
1082 include 'COMMON.IOUNITS'
1083 include 'COMMON.NAMES'
1084 include 'COMMON.TIME1'
1085 double precision stochforcvec(MAXRES6)
1086 common /stochcalc/ stochforcvec
1088 c Compute friction and stochastic forces
1092 time_fric=time_fric+MPI_Wtime()-time00
1094 call stochastic_force(stochforcvec)
1095 time_stoch=time_stoch+MPI_Wtime()-time00
1097 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1098 c forces (d_as_work)
1100 call ginv_mult(fric_work, d_af_work)
1101 call ginv_mult(stochforcvec, d_as_work)
1104 c---------------------------------------------------------------------
1105 subroutine sddir_verlet1
1106 c Applying velocity Verlet algorithm - step 1 to velocities
1107 implicit real*8 (a-h,o-z)
1108 include 'DIMENSIONS'
1109 include 'COMMON.CONTROL'
1110 include 'COMMON.VAR'
1113 include 'COMMON.LANGEVIN'
1115 include 'COMMON.LANGEVIN.lang0'
1117 include 'COMMON.CHAIN'
1118 include 'COMMON.DERIV'
1119 include 'COMMON.GEO'
1120 include 'COMMON.LOCAL'
1121 include 'COMMON.INTERACT'
1122 include 'COMMON.IOUNITS'
1123 include 'COMMON.NAMES'
1124 c Revised 3/31/05 AL: correlation between random contributions to
1125 c position and velocity increments included.
1126 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1127 double precision adt,adt2
1129 c Add the contribution from BOTH friction and stochastic force to the
1130 c coordinates, but ONLY the contribution from the friction forces to velocities
1133 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1134 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1135 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1136 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1137 d_t(j,0)=d_t_old(j,0)+adt
1142 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1143 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1144 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1145 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1146 d_t(j,i)=d_t_old(j,i)+adt
1151 if (itype(i).ne.10 .and. itype(i).ne.21) then
1154 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1155 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1156 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1157 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1158 d_t(j,inres)=d_t_old(j,inres)+adt
1165 c---------------------------------------------------------------------
1166 subroutine sddir_verlet2
1167 c Calculating the adjusted velocities for accelerations
1168 implicit real*8 (a-h,o-z)
1169 include 'DIMENSIONS'
1170 include 'COMMON.CONTROL'
1171 include 'COMMON.VAR'
1174 include 'COMMON.LANGEVIN'
1176 include 'COMMON.LANGEVIN.lang0'
1178 include 'COMMON.CHAIN'
1179 include 'COMMON.DERIV'
1180 include 'COMMON.GEO'
1181 include 'COMMON.LOCAL'
1182 include 'COMMON.INTERACT'
1183 include 'COMMON.IOUNITS'
1184 include 'COMMON.NAMES'
1185 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1186 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1187 c Revised 3/31/05 AL: correlation between random contributions to
1188 c position and velocity increments included.
1189 c The correlation coefficients are calculated at low-friction limit.
1190 c Also, friction forces are now not calculated with new velocities.
1192 c call friction_force
1193 call stochastic_force(stochforcvec)
1195 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1196 c forces (d_as_work)
1198 call ginv_mult(stochforcvec, d_as_work1)
1204 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1205 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1210 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1211 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1216 if (itype(i).ne.10 .and. itype(i).ne.21) then
1219 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1220 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1221 & +cos60*d_as_work1(ind+j))*d_time
1228 c---------------------------------------------------------------------
1229 subroutine max_accel
1231 c Find the maximum difference in the accelerations of the the sites
1232 c at the beginning and the end of the time step.
1234 implicit real*8 (a-h,o-z)
1235 include 'DIMENSIONS'
1236 include 'COMMON.CONTROL'
1237 include 'COMMON.VAR'
1239 include 'COMMON.CHAIN'
1240 include 'COMMON.DERIV'
1241 include 'COMMON.GEO'
1242 include 'COMMON.LOCAL'
1243 include 'COMMON.INTERACT'
1244 include 'COMMON.IOUNITS'
1245 double precision aux(3),accel(3),accel_old(3),dacc
1247 c aux(j)=d_a(j,0)-d_a_old(j,0)
1248 accel_old(j)=d_a_old(j,0)
1255 c 7/3/08 changed to asymmetric difference
1257 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1258 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1259 accel(j)=accel(j)+0.5d0*d_a(j,i)
1260 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1261 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1262 dacc=dabs(accel(j)-accel_old(j))
1263 if (dacc.gt.amax) amax=dacc
1271 accel_old(j)=d_a_old(j,0)
1276 accel_old(j)=accel_old(j)+d_a_old(j,1)
1277 accel(j)=accel(j)+d_a(j,1)
1281 if (itype(i).ne.10 .and. itype(i).ne.21) then
1283 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1284 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1285 accel(j)=accel(j)+d_a(j,i+nres)
1289 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1290 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1291 dacc=dabs(accel(j)-accel_old(j))
1292 if (dacc.gt.amax) amax=dacc
1296 accel_old(j)=accel_old(j)+d_a_old(j,i)
1297 accel(j)=accel(j)+d_a(j,i)
1298 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1303 c---------------------------------------------------------------------
1304 subroutine predict_edrift(epdrift)
1306 c Predict the drift of the potential energy
1308 implicit real*8 (a-h,o-z)
1309 include 'DIMENSIONS'
1310 include 'COMMON.CONTROL'
1311 include 'COMMON.VAR'
1313 include 'COMMON.CHAIN'
1314 include 'COMMON.DERIV'
1315 include 'COMMON.GEO'
1316 include 'COMMON.LOCAL'
1317 include 'COMMON.INTERACT'
1318 include 'COMMON.IOUNITS'
1319 include 'COMMON.MUCA'
1320 double precision epdrift,epdriftij
1321 c Drift of the potential energy
1327 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1328 if (lmuca) epdriftij=epdriftij*factor
1329 c write (iout,*) "back",i,j,epdriftij
1330 if (epdriftij.gt.epdrift) epdrift=epdriftij
1334 if (itype(i).ne.10 .and. itype(i).ne.21) then
1337 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1338 if (lmuca) epdriftij=epdriftij*factor
1339 c write (iout,*) "side",i,j,epdriftij
1340 if (epdriftij.gt.epdrift) epdrift=epdriftij
1344 epdrift=0.5d0*epdrift*d_time*d_time
1345 c write (iout,*) "epdrift",epdrift
1348 c-----------------------------------------------------------------------
1349 subroutine verlet_bath
1351 c Coupling to the thermostat by using the Berendsen algorithm
1353 implicit real*8 (a-h,o-z)
1354 include 'DIMENSIONS'
1355 include 'COMMON.CONTROL'
1356 include 'COMMON.VAR'
1358 include 'COMMON.CHAIN'
1359 include 'COMMON.DERIV'
1360 include 'COMMON.GEO'
1361 include 'COMMON.LOCAL'
1362 include 'COMMON.INTERACT'
1363 include 'COMMON.IOUNITS'
1364 include 'COMMON.NAMES'
1365 double precision T_half,fact
1367 T_half=2.0d0/(dimen3*Rb)*EK
1368 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1369 c write(iout,*) "T_half", T_half
1370 c write(iout,*) "EK", EK
1371 c write(iout,*) "fact", fact
1373 d_t(j,0)=fact*d_t(j,0)
1377 d_t(j,i)=fact*d_t(j,i)
1381 if (itype(i).ne.10 .and. itype(i).ne.21) then
1384 d_t(j,inres)=fact*d_t(j,inres)
1390 c---------------------------------------------------------
1392 c Set up the initial conditions of a MD simulation
1393 implicit real*8 (a-h,o-z)
1394 include 'DIMENSIONS'
1398 integer IERROR,ERRCODE
1400 include 'COMMON.SETUP'
1401 include 'COMMON.CONTROL'
1402 include 'COMMON.VAR'
1405 include 'COMMON.LANGEVIN'
1407 include 'COMMON.LANGEVIN.lang0'
1409 include 'COMMON.CHAIN'
1410 include 'COMMON.DERIV'
1411 include 'COMMON.GEO'
1412 include 'COMMON.LOCAL'
1413 include 'COMMON.INTERACT'
1414 include 'COMMON.IOUNITS'
1415 include 'COMMON.NAMES'
1416 include 'COMMON.REMD'
1417 real*8 energia_long(0:n_ene),
1418 & energia_short(0:n_ene),vcm(3),incr(3)
1419 double precision cm(3),L(3),xv,sigv,lowb,highb
1420 double precision varia(maxvar)
1428 c write(iout,*) "d_time", d_time
1429 c Compute the standard deviations of stochastic forces for Langevin dynamics
1430 c if the friction coefficients do not depend on surface area
1431 if (lang.gt.0 .and. .not.surfarea) then
1433 stdforcp(i)=stdfp*dsqrt(gamp)
1436 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
1439 c Open the pdb file for snapshotshots
1442 if (ilen(tmpdir).gt.0)
1443 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1444 & liczba(:ilen(liczba))//".pdb")
1446 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1450 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1451 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1452 & liczba(:ilen(liczba))//".x")
1453 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1456 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1457 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1458 & liczba(:ilen(liczba))//".cx")
1459 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1465 if (ilen(tmpdir).gt.0)
1466 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1467 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1469 if (ilen(tmpdir).gt.0)
1470 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1471 cartname=prefix(:ilen(prefix))//"_MD.cx"
1475 write (qstr,'(256(1h ))')
1478 iq = qinfrag(i,iset)*10
1479 iw = wfrag(i,iset)/100
1481 if(me.eq.king.or..not.out1file)
1482 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1483 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1488 iq = qinpair(i,iset)*10
1489 iw = wpair(i,iset)/100
1491 if(me.eq.king.or..not.out1file)
1492 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1493 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1497 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1499 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1501 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1503 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1507 if (restart1file) then
1509 & inquire(file=mremd_rst_name,exist=file_exist)
1510 write (*,*) me," Before broadcast: file_exist",file_exist
1511 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1513 write (*,*) me," After broadcast: file_exist",file_exist
1514 c inquire(file=mremd_rst_name,exist=file_exist)
1515 if(me.eq.king.or..not.out1file)
1516 & write(iout,*) "Initial state read by master and distributed"
1518 if (ilen(tmpdir).gt.0)
1519 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1520 & //liczba(:ilen(liczba))//'.rst')
1521 inquire(file=rest2name,exist=file_exist)
1524 if(.not.restart1file) then
1525 if(me.eq.king.or..not.out1file)
1526 & write(iout,*) "Initial state will be read from file ",
1527 & rest2name(:ilen(rest2name))
1530 call rescale_weights(t_bath)
1532 if(me.eq.king.or..not.out1file)then
1533 if (restart1file) then
1534 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1537 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1540 write(iout,*) "Initial velocities randomly generated"
1546 c Generate initial velocities
1547 if(me.eq.king.or..not.out1file)
1548 & write(iout,*) "Initial velocities randomly generated"
1552 c rest2name = prefix(:ilen(prefix))//'.rst'
1553 if(me.eq.king.or..not.out1file)then
1554 write (iout,*) "Initial velocities"
1556 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1557 & (d_t(j,i+nres),j=1,3)
1559 c Zeroing the total angular momentum of the system
1560 write(iout,*) "Calling the zero-angular
1561 & momentum subroutine"
1564 c Getting the potential energy and forces and velocities and accelerations
1566 c write (iout,*) "velocity of the center of the mass:"
1567 c write (iout,*) (vcm(j),j=1,3)
1569 d_t(j,0)=d_t(j,0)-vcm(j)
1571 c Removing the velocity of the center of mass
1573 if(me.eq.king.or..not.out1file)then
1574 write (iout,*) "vcm right after adjustment:"
1575 write (iout,*) (vcm(j),j=1,3)
1579 if(iranconf.ne.0) then
1581 print *, 'Calling OVERLAP_SC'
1582 call overlap_sc(fail)
1586 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1587 print *,'SC_move',nft_sc,etot
1588 if(me.eq.king.or..not.out1file)
1589 & write(iout,*) 'SC_move',nft_sc,etot
1593 print *, 'Calling MINIM_DC'
1594 call minim_dc(etot,iretcode,nfun)
1596 call geom_to_var(nvar,varia)
1597 print *,'Calling MINIMIZE.'
1598 call minimize(etot,varia,iretcode,nfun)
1599 call var_to_geom(nvar,varia)
1601 if(me.eq.king.or..not.out1file)
1602 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1605 call chainbuild_cart
1608 call verlet_bath(EK)
1610 kinetic_T=2.0d0/(dimen3*Rb)*EK
1611 if(me.eq.king.or..not.out1file)then
1621 call etotal(potEcomp)
1622 if (large) call enerprint(potEcomp)
1625 t_etotal=t_etotal+MPI_Wtime()-tt0
1627 t_etotal=t_etotal+tcpu()-tt0
1634 if (amax*d_time .gt. dvmax) then
1635 d_time=d_time*dvmax/amax
1636 if(me.eq.king.or..not.out1file) write (iout,*)
1637 & "Time step reduced to",d_time,
1638 & " because of too large initial acceleration."
1640 if(me.eq.king.or..not.out1file)then
1641 write(iout,*) "Potential energy and its components"
1642 call enerprint(potEcomp)
1643 c write(iout,*) (potEcomp(i),i=0,n_ene)
1645 potE=potEcomp(0)-potEcomp(20)
1648 if (ntwe.ne.0) call statout(itime)
1649 if(me.eq.king.or..not.out1file)
1650 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1651 & " Kinetic energy",EK," Potential energy",potE,
1652 & " Total energy",totE," Maximum acceleration ",
1655 write (iout,*) "Initial coordinates"
1657 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1658 & (c(j,i+nres),j=1,3)
1660 write (iout,*) "Initial dC"
1662 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1663 & (dc(j,i+nres),j=1,3)
1665 write (iout,*) "Initial velocities"
1666 write (iout,"(13x,' backbone ',23x,' side chain')")
1668 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1669 & (d_t(j,i+nres),j=1,3)
1671 write (iout,*) "Initial accelerations"
1673 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1674 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1675 & (d_a(j,i+nres),j=1,3)
1681 d_t_old(j,i)=d_t(j,i)
1682 d_a_old(j,i)=d_a(j,i)
1684 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1693 call etotal_short(energia_short)
1694 if (large) call enerprint(potEcomp)
1697 t_eshort=t_eshort+MPI_Wtime()-tt0
1699 t_eshort=t_eshort+tcpu()-tt0
1704 if(.not.out1file .and. large) then
1705 write (iout,*) "energia_long",energia_long(0),
1706 & " energia_short",energia_short(0),
1707 & " total",energia_long(0)+energia_short(0)
1708 write (iout,*) "Initial fast-force accelerations"
1710 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1711 & (d_a(j,i+nres),j=1,3)
1714 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1717 d_a_short(j,i)=d_a(j,i)
1726 call etotal_long(energia_long)
1727 if (large) call enerprint(potEcomp)
1730 t_elong=t_elong+MPI_Wtime()-tt0
1732 t_elong=t_elong+tcpu()-tt0
1737 if(.not.out1file .and. large) then
1738 write (iout,*) "energia_long",energia_long(0)
1739 write (iout,*) "Initial slow-force accelerations"
1741 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1742 & (d_a(j,i+nres),j=1,3)
1746 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1748 t_enegrad=t_enegrad+tcpu()-tt0
1753 c-----------------------------------------------------------
1754 subroutine random_vel
1755 implicit real*8 (a-h,o-z)
1756 include 'DIMENSIONS'
1757 include 'COMMON.CONTROL'
1758 include 'COMMON.VAR'
1761 include 'COMMON.LANGEVIN'
1763 include 'COMMON.LANGEVIN.lang0'
1765 include 'COMMON.CHAIN'
1766 include 'COMMON.DERIV'
1767 include 'COMMON.GEO'
1768 include 'COMMON.LOCAL'
1769 include 'COMMON.INTERACT'
1770 include 'COMMON.IOUNITS'
1771 include 'COMMON.NAMES'
1772 include 'COMMON.TIME1'
1773 double precision xv,sigv,lowb,highb
1774 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1775 c First generate velocities in the eigenspace of the G matrix
1776 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1783 sigv=dsqrt((Rb*t_bath)/geigen(i))
1786 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1787 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1788 c & " d_t_work_new",d_t_work_new(ii)
1797 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1800 c write (iout,*) "Ek from eigenvectors",Ek1
1802 c Transform velocities to UNRES coordinate space
1808 d_t_work(ind)=d_t_work(ind)
1809 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1811 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1815 c Transfer to the d_t vector
1817 d_t(j,0)=d_t_work(j)
1823 d_t(j,i)=d_t_work(ind)
1827 if (itype(i).ne.10 .and. itype(i).ne.21) then
1830 d_t(j,i+nres)=d_t_work(ind)
1835 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1836 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1841 c-----------------------------------------------------------
1842 subroutine sd_verlet_p_setup
1843 c Sets up the parameters of stochastic Verlet algorithm
1844 implicit real*8 (a-h,o-z)
1845 include 'DIMENSIONS'
1849 include 'COMMON.CONTROL'
1850 include 'COMMON.VAR'
1853 include 'COMMON.LANGEVIN'
1855 include 'COMMON.LANGEVIN.lang0'
1857 include 'COMMON.CHAIN'
1858 include 'COMMON.DERIV'
1859 include 'COMMON.GEO'
1860 include 'COMMON.LOCAL'
1861 include 'COMMON.INTERACT'
1862 include 'COMMON.IOUNITS'
1863 include 'COMMON.NAMES'
1864 include 'COMMON.TIME1'
1865 double precision emgdt(MAXRES6),
1866 & pterm,vterm,rho,rhoc,vsig,
1867 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1868 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1869 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1870 logical lprn /.false./
1871 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1872 double precision ktm
1879 c AL 8/17/04 Code adapted from tinker
1881 c Get the frictional and random terms for stochastic dynamics in the
1882 c eigenspace of mass-scaled UNRES friction matrix
1885 gdt = fricgam(i) * d_time
1887 c Stochastic dynamics reduces to simple MD for zero friction
1889 if (gdt .le. zero) then
1890 pfric_vec(i) = 1.0d0
1891 vfric_vec(i) = d_time
1892 afric_vec(i) = 0.5d0 * d_time * d_time
1893 prand_vec(i) = 0.0d0
1894 vrand_vec1(i) = 0.0d0
1895 vrand_vec2(i) = 0.0d0
1897 c Analytical expressions when friction coefficient is large
1900 if (gdt .ge. gdt_radius) then
1903 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1904 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1905 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1906 vterm = 1.0d0 - egdt**2
1907 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1909 c Use series expansions when friction coefficient is small
1920 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1921 & - gdt5/120.0d0 + gdt6/720.0d0
1922 & - gdt7/5040.0d0 + gdt8/40320.0d0
1923 & - gdt9/362880.0d0) / fricgam(i)**2
1924 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1925 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1926 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1927 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1928 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1929 & + 127.0d0*gdt9/90720.0d0
1930 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1931 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1932 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1933 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1934 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1935 & - 17.0d0*gdt2/1280.0d0
1936 & + 17.0d0*gdt3/6144.0d0
1937 & + 40967.0d0*gdt4/34406400.0d0
1938 & - 57203.0d0*gdt5/275251200.0d0
1939 & - 1429487.0d0*gdt6/13212057600.0d0)
1942 c Compute the scaling factors of random terms for the nonzero friction case
1944 ktm = 0.5d0*d_time/fricgam(i)
1945 psig = dsqrt(ktm*pterm) / fricgam(i)
1946 vsig = dsqrt(ktm*vterm)
1947 rhoc = dsqrt(1.0d0 - rho*rho)
1949 vrand_vec1(i) = vsig * rho
1950 vrand_vec2(i) = vsig * rhoc
1955 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1958 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1959 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1963 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1966 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1967 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1968 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1969 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1970 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1971 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1974 t_sdsetup=t_sdsetup+MPI_Wtime()
1976 t_sdsetup=t_sdsetup+tcpu()-tt0
1980 c-------------------------------------------------------------
1981 subroutine eigtransf1(n,ndim,ab,d,c)
1984 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
1990 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
1996 c-------------------------------------------------------------
1997 subroutine eigtransf(n,ndim,a,b,d,c)
2000 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2006 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2012 c-------------------------------------------------------------
2013 subroutine sd_verlet1
2014 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2015 implicit real*8 (a-h,o-z)
2016 include 'DIMENSIONS'
2017 include 'COMMON.CONTROL'
2018 include 'COMMON.VAR'
2021 include 'COMMON.LANGEVIN'
2023 include 'COMMON.LANGEVIN.lang0'
2025 include 'COMMON.CHAIN'
2026 include 'COMMON.DERIV'
2027 include 'COMMON.GEO'
2028 include 'COMMON.LOCAL'
2029 include 'COMMON.INTERACT'
2030 include 'COMMON.IOUNITS'
2031 include 'COMMON.NAMES'
2032 double precision stochforcvec(MAXRES6)
2033 common /stochcalc/ stochforcvec
2034 logical lprn /.false./
2036 c write (iout,*) "dc_old"
2038 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2039 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2042 dc_work(j)=dc_old(j,0)
2043 d_t_work(j)=d_t_old(j,0)
2044 d_a_work(j)=d_a_old(j,0)
2049 dc_work(ind+j)=dc_old(j,i)
2050 d_t_work(ind+j)=d_t_old(j,i)
2051 d_a_work(ind+j)=d_a_old(j,i)
2056 if (itype(i).ne.10 .and. itype(i).ne.21) then
2058 dc_work(ind+j)=dc_old(j,i+nres)
2059 d_t_work(ind+j)=d_t_old(j,i+nres)
2060 d_a_work(ind+j)=d_a_old(j,i+nres)
2068 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2072 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2073 & vfric_mat(i,j),afric_mat(i,j),
2074 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2082 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2083 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2084 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2085 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2087 d_t_work_new(i)=ddt1+0.5d0*ddt2
2088 d_t_work(i)=ddt1+ddt2
2093 d_t(j,0)=d_t_work(j)
2098 dc(j,i)=dc_work(ind+j)
2099 d_t(j,i)=d_t_work(ind+j)
2104 if (itype(i).ne.10) then
2107 dc(j,inres)=dc_work(ind+j)
2108 d_t(j,inres)=d_t_work(ind+j)
2115 c--------------------------------------------------------------------------
2116 subroutine sd_verlet2
2117 c Calculating the adjusted velocities for accelerations
2118 implicit real*8 (a-h,o-z)
2119 include 'DIMENSIONS'
2120 include 'COMMON.CONTROL'
2121 include 'COMMON.VAR'
2124 include 'COMMON.LANGEVIN'
2126 include 'COMMON.LANGEVIN.lang0'
2128 include 'COMMON.CHAIN'
2129 include 'COMMON.DERIV'
2130 include 'COMMON.GEO'
2131 include 'COMMON.LOCAL'
2132 include 'COMMON.INTERACT'
2133 include 'COMMON.IOUNITS'
2134 include 'COMMON.NAMES'
2135 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2136 common /stochcalc/ stochforcvec
2138 c Compute the stochastic forces which contribute to velocity change
2140 call stochastic_force(stochforcvecV)
2147 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2148 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2149 & vrand_mat2(i,j)*stochforcvecV(j)
2151 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2155 d_t(j,0)=d_t_work(j)
2160 d_t(j,i)=d_t_work(ind+j)
2165 if (itype(i).ne.10 .and. itype(i).ne.21) then
2168 d_t(j,inres)=d_t_work(ind+j)
2175 c-----------------------------------------------------------
2176 subroutine sd_verlet_ciccotti_setup
2177 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2179 implicit real*8 (a-h,o-z)
2180 include 'DIMENSIONS'
2184 include 'COMMON.CONTROL'
2185 include 'COMMON.VAR'
2188 include 'COMMON.LANGEVIN'
2190 include 'COMMON.LANGEVIN.lang0'
2192 include 'COMMON.CHAIN'
2193 include 'COMMON.DERIV'
2194 include 'COMMON.GEO'
2195 include 'COMMON.LOCAL'
2196 include 'COMMON.INTERACT'
2197 include 'COMMON.IOUNITS'
2198 include 'COMMON.NAMES'
2199 include 'COMMON.TIME1'
2200 double precision emgdt(MAXRES6),
2201 & pterm,vterm,rho,rhoc,vsig,
2202 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2203 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2204 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2205 logical lprn /.false./
2206 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2207 double precision ktm
2214 c AL 8/17/04 Code adapted from tinker
2216 c Get the frictional and random terms for stochastic dynamics in the
2217 c eigenspace of mass-scaled UNRES friction matrix
2220 write (iout,*) "i",i," fricgam",fricgam(i)
2221 gdt = fricgam(i) * d_time
2223 c Stochastic dynamics reduces to simple MD for zero friction
2225 if (gdt .le. zero) then
2226 pfric_vec(i) = 1.0d0
2227 vfric_vec(i) = d_time
2228 afric_vec(i) = 0.5d0*d_time*d_time
2229 prand_vec(i) = afric_vec(i)
2230 vrand_vec2(i) = vfric_vec(i)
2232 c Analytical expressions when friction coefficient is large
2237 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2238 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2239 prand_vec(i) = afric_vec(i)
2240 vrand_vec2(i) = vfric_vec(i)
2242 c Compute the scaling factors of random terms for the nonzero friction case
2244 c ktm = 0.5d0*d_time/fricgam(i)
2245 c psig = dsqrt(ktm*pterm) / fricgam(i)
2246 c vsig = dsqrt(ktm*vterm)
2247 c prand_vec(i) = psig*afric_vec(i)
2248 c vrand_vec2(i) = vsig*vfric_vec(i)
2253 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2256 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2257 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2261 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2263 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2264 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2265 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2266 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2267 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2269 t_sdsetup=t_sdsetup+MPI_Wtime()
2271 t_sdsetup=t_sdsetup+tcpu()-tt0
2275 c-------------------------------------------------------------
2276 subroutine sd_verlet1_ciccotti
2277 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2278 implicit real*8 (a-h,o-z)
2279 include 'DIMENSIONS'
2283 include 'COMMON.CONTROL'
2284 include 'COMMON.VAR'
2287 include 'COMMON.LANGEVIN'
2289 include 'COMMON.LANGEVIN.lang0'
2291 include 'COMMON.CHAIN'
2292 include 'COMMON.DERIV'
2293 include 'COMMON.GEO'
2294 include 'COMMON.LOCAL'
2295 include 'COMMON.INTERACT'
2296 include 'COMMON.IOUNITS'
2297 include 'COMMON.NAMES'
2298 double precision stochforcvec(MAXRES6)
2299 common /stochcalc/ stochforcvec
2300 logical lprn /.false./
2302 c write (iout,*) "dc_old"
2304 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2305 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2308 dc_work(j)=dc_old(j,0)
2309 d_t_work(j)=d_t_old(j,0)
2310 d_a_work(j)=d_a_old(j,0)
2315 dc_work(ind+j)=dc_old(j,i)
2316 d_t_work(ind+j)=d_t_old(j,i)
2317 d_a_work(ind+j)=d_a_old(j,i)
2322 if (itype(i).ne.10 .and. itype(i).ne.21) then
2324 dc_work(ind+j)=dc_old(j,i+nres)
2325 d_t_work(ind+j)=d_t_old(j,i+nres)
2326 d_a_work(ind+j)=d_a_old(j,i+nres)
2335 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2339 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2340 & vfric_mat(i,j),afric_mat(i,j),
2341 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2349 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2350 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2351 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2352 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2354 d_t_work_new(i)=ddt1+0.5d0*ddt2
2355 d_t_work(i)=ddt1+ddt2
2360 d_t(j,0)=d_t_work(j)
2365 dc(j,i)=dc_work(ind+j)
2366 d_t(j,i)=d_t_work(ind+j)
2371 if (itype(i).ne.10 .and. itype(i).ne.21) then
2374 dc(j,inres)=dc_work(ind+j)
2375 d_t(j,inres)=d_t_work(ind+j)
2382 c--------------------------------------------------------------------------
2383 subroutine sd_verlet2_ciccotti
2384 c Calculating the adjusted velocities for accelerations
2385 implicit real*8 (a-h,o-z)
2386 include 'DIMENSIONS'
2387 include 'COMMON.CONTROL'
2388 include 'COMMON.VAR'
2391 include 'COMMON.LANGEVIN'
2393 include 'COMMON.LANGEVIN.lang0'
2395 include 'COMMON.CHAIN'
2396 include 'COMMON.DERIV'
2397 include 'COMMON.GEO'
2398 include 'COMMON.LOCAL'
2399 include 'COMMON.INTERACT'
2400 include 'COMMON.IOUNITS'
2401 include 'COMMON.NAMES'
2402 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2403 common /stochcalc/ stochforcvec
2405 c Compute the stochastic forces which contribute to velocity change
2407 call stochastic_force(stochforcvecV)
2414 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2415 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2416 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2418 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2422 d_t(j,0)=d_t_work(j)
2427 d_t(j,i)=d_t_work(ind+j)
2432 if (itype(i).ne.10 .and. itype(i).ne.21) then
2435 d_t(j,inres)=d_t_work(ind+j)