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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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 c write (iout,*) i,dacc
1264 if (dacc.gt.amax) amax=dacc
1272 accel_old(j)=d_a_old(j,0)
1277 accel_old(j)=accel_old(j)+d_a_old(j,1)
1278 accel(j)=accel(j)+d_a(j,1)
1282 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1284 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1285 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1286 accel(j)=accel(j)+d_a(j,i+nres)
1290 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1291 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1292 dacc=dabs(accel(j)-accel_old(j))
1293 c write (iout,*) "side-chain",i,dacc
1294 if (dacc.gt.amax) amax=dacc
1298 accel_old(j)=accel_old(j)+d_a_old(j,i)
1299 accel(j)=accel(j)+d_a(j,i)
1300 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1305 c---------------------------------------------------------------------
1306 subroutine predict_edrift(epdrift)
1308 c Predict the drift of the potential energy
1310 implicit real*8 (a-h,o-z)
1311 include 'DIMENSIONS'
1312 include 'COMMON.CONTROL'
1313 include 'COMMON.VAR'
1315 include 'COMMON.CHAIN'
1316 include 'COMMON.DERIV'
1317 include 'COMMON.GEO'
1318 include 'COMMON.LOCAL'
1319 include 'COMMON.INTERACT'
1320 include 'COMMON.IOUNITS'
1321 include 'COMMON.MUCA'
1322 double precision epdrift,epdriftij
1323 c Drift of the potential energy
1329 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1330 if (lmuca) epdriftij=epdriftij*factor
1331 c write (iout,*) "back",i,j,epdriftij
1332 if (epdriftij.gt.epdrift) epdrift=epdriftij
1336 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1339 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1340 if (lmuca) epdriftij=epdriftij*factor
1341 c write (iout,*) "side",i,j,epdriftij
1342 if (epdriftij.gt.epdrift) epdrift=epdriftij
1346 epdrift=0.5d0*epdrift*d_time*d_time
1347 c write (iout,*) "epdrift",epdrift
1350 c-----------------------------------------------------------------------
1351 subroutine verlet_bath
1353 c Coupling to the thermostat by using the Berendsen algorithm
1355 implicit real*8 (a-h,o-z)
1356 include 'DIMENSIONS'
1357 include 'COMMON.CONTROL'
1358 include 'COMMON.VAR'
1360 include 'COMMON.CHAIN'
1361 include 'COMMON.DERIV'
1362 include 'COMMON.GEO'
1363 include 'COMMON.LOCAL'
1364 include 'COMMON.INTERACT'
1365 include 'COMMON.IOUNITS'
1366 include 'COMMON.NAMES'
1367 double precision T_half,fact
1369 T_half=2.0d0/(dimen3*Rb)*EK
1370 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1371 c write(iout,*) "T_half", T_half
1372 c write(iout,*) "EK", EK
1373 c write(iout,*) "fact", fact
1375 d_t(j,0)=fact*d_t(j,0)
1379 d_t(j,i)=fact*d_t(j,i)
1383 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1386 d_t(j,inres)=fact*d_t(j,inres)
1392 c---------------------------------------------------------
1394 c Set up the initial conditions of a MD simulation
1395 implicit real*8 (a-h,o-z)
1396 include 'DIMENSIONS'
1400 integer IERROR,ERRCODE
1402 include 'COMMON.SETUP'
1403 include 'COMMON.CONTROL'
1404 include 'COMMON.VAR'
1407 include 'COMMON.LANGEVIN'
1409 include 'COMMON.LANGEVIN.lang0'
1411 include 'COMMON.CHAIN'
1412 include 'COMMON.DERIV'
1413 include 'COMMON.GEO'
1414 include 'COMMON.LOCAL'
1415 include 'COMMON.INTERACT'
1416 include 'COMMON.IOUNITS'
1417 include 'COMMON.NAMES'
1418 include 'COMMON.REMD'
1419 real*8 energia_long(0:n_ene),
1420 & energia_short(0:n_ene),vcm(3),incr(3)
1421 double precision cm(3),L(3),xv,sigv,lowb,highb
1422 double precision varia(maxvar)
1430 c write(iout,*) "d_time", d_time
1431 c Compute the standard deviations of stochastic forces for Langevin dynamics
1432 c if the friction coefficients do not depend on surface area
1433 if (lang.gt.0 .and. .not.surfarea) then
1435 stdforcp(i)=stdfp*dsqrt(gamp)
1438 stdforcsc(i)=stdfsc(iabs(itype(i)))
1439 & *dsqrt(gamsc(iabs(itype(i))))
1442 c Open the pdb file for snapshotshots
1445 if (ilen(tmpdir).gt.0)
1446 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1447 & liczba(:ilen(liczba))//".pdb")
1449 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1453 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1454 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1455 & liczba(:ilen(liczba))//".x")
1456 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1459 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1460 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1461 & liczba(:ilen(liczba))//".cx")
1462 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1468 if (ilen(tmpdir).gt.0)
1469 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1470 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1472 if (ilen(tmpdir).gt.0)
1473 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1474 cartname=prefix(:ilen(prefix))//"_MD.cx"
1478 write (qstr,'(256(1h ))')
1481 iq = qinfrag(i,iset)*10
1482 iw = wfrag(i,iset)/100
1484 if(me.eq.king.or..not.out1file)
1485 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1486 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1491 iq = qinpair(i,iset)*10
1492 iw = wpair(i,iset)/100
1494 if(me.eq.king.or..not.out1file)
1495 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1496 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1500 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1502 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1504 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1506 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1510 if (restart1file) then
1512 & inquire(file=mremd_rst_name,exist=file_exist)
1513 write (*,*) me," Before broadcast: file_exist",file_exist
1514 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1516 write (*,*) me," After broadcast: file_exist",file_exist
1517 c inquire(file=mremd_rst_name,exist=file_exist)
1518 if(me.eq.king.or..not.out1file)
1519 & write(iout,*) "Initial state read by master and distributed"
1521 if (ilen(tmpdir).gt.0)
1522 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1523 & //liczba(:ilen(liczba))//'.rst')
1524 inquire(file=rest2name,exist=file_exist)
1527 if(.not.restart1file) then
1528 if(me.eq.king.or..not.out1file)
1529 & write(iout,*) "Initial state will be read from file ",
1530 & rest2name(:ilen(rest2name))
1533 call rescale_weights(t_bath)
1535 if(me.eq.king.or..not.out1file)then
1536 if (restart1file) then
1537 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1540 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1543 write(iout,*) "Initial velocities randomly generated"
1549 c Generate initial velocities
1550 if(me.eq.king.or..not.out1file)
1551 & write(iout,*) "Initial velocities randomly generated"
1555 c rest2name = prefix(:ilen(prefix))//'.rst'
1556 if(me.eq.king.or..not.out1file)then
1557 write (iout,*) "Initial velocities"
1559 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1560 & (d_t(j,i+nres),j=1,3)
1562 c Zeroing the total angular momentum of the system
1563 write(iout,*) "Calling the zero-angular
1564 & momentum subroutine"
1567 c Getting the potential energy and forces and velocities and accelerations
1569 c write (iout,*) "velocity of the center of the mass:"
1570 c write (iout,*) (vcm(j),j=1,3)
1572 d_t(j,0)=d_t(j,0)-vcm(j)
1574 c Removing the velocity of the center of mass
1576 if(me.eq.king.or..not.out1file)then
1577 write (iout,*) "vcm right after adjustment:"
1578 write (iout,*) (vcm(j),j=1,3)
1582 if(iranconf.ne.0) then
1584 print *, 'Calling OVERLAP_SC'
1585 call overlap_sc(fail)
1589 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1590 print *,'SC_move',nft_sc,etot
1591 if(me.eq.king.or..not.out1file)
1592 & write(iout,*) 'SC_move',nft_sc,etot
1596 print *, 'Calling MINIM_DC'
1597 call minim_dc(etot,iretcode,nfun)
1599 call geom_to_var(nvar,varia)
1600 print *,'Calling MINIMIZE.'
1601 call minimize(etot,varia,iretcode,nfun)
1602 call var_to_geom(nvar,varia)
1604 if(me.eq.king.or..not.out1file)
1605 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1608 call chainbuild_cart
1613 kinetic_T=2.0d0/(dimen3*Rb)*EK
1614 if(me.eq.king.or..not.out1file)then
1624 call etotal(potEcomp)
1625 if (large) call enerprint(potEcomp)
1628 t_etotal=t_etotal+MPI_Wtime()-tt0
1630 t_etotal=t_etotal+tcpu()-tt0
1637 if (amax*d_time .gt. dvmax) then
1638 d_time=d_time*dvmax/amax
1639 if(me.eq.king.or..not.out1file) write (iout,*)
1640 & "Time step reduced to",d_time,
1641 & " because of too large initial acceleration."
1643 if(me.eq.king.or..not.out1file)then
1644 write(iout,*) "Potential energy and its components"
1645 call enerprint(potEcomp)
1646 c write(iout,*) (potEcomp(i),i=0,n_ene)
1648 potE=potEcomp(0)-potEcomp(20)
1651 if (ntwe.ne.0) call statout(itime)
1652 if(me.eq.king.or..not.out1file)
1653 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1654 & " Kinetic energy",EK," Potential energy",potE,
1655 & " Total energy",totE," Maximum acceleration ",
1658 write (iout,*) "Initial coordinates"
1660 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1661 & (c(j,i+nres),j=1,3)
1663 write (iout,*) "Initial dC"
1665 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1666 & (dc(j,i+nres),j=1,3)
1668 write (iout,*) "Initial velocities"
1669 write (iout,"(13x,' backbone ',23x,' side chain')")
1671 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1672 & (d_t(j,i+nres),j=1,3)
1674 write (iout,*) "Initial accelerations"
1676 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1677 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1678 & (d_a(j,i+nres),j=1,3)
1684 d_t_old(j,i)=d_t(j,i)
1685 d_a_old(j,i)=d_a(j,i)
1687 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1696 call etotal_short(energia_short)
1697 if (large) call enerprint(potEcomp)
1700 t_eshort=t_eshort+MPI_Wtime()-tt0
1702 t_eshort=t_eshort+tcpu()-tt0
1707 if(.not.out1file .and. large) then
1708 write (iout,*) "energia_long",energia_long(0),
1709 & " energia_short",energia_short(0),
1710 & " total",energia_long(0)+energia_short(0)
1711 write (iout,*) "Initial fast-force accelerations"
1713 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1714 & (d_a(j,i+nres),j=1,3)
1717 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1720 d_a_short(j,i)=d_a(j,i)
1729 call etotal_long(energia_long)
1730 if (large) call enerprint(potEcomp)
1733 t_elong=t_elong+MPI_Wtime()-tt0
1735 t_elong=t_elong+tcpu()-tt0
1740 if(.not.out1file .and. large) then
1741 write (iout,*) "energia_long",energia_long(0)
1742 write (iout,*) "Initial slow-force accelerations"
1744 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1745 & (d_a(j,i+nres),j=1,3)
1749 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1751 t_enegrad=t_enegrad+tcpu()-tt0
1756 c-----------------------------------------------------------
1757 subroutine random_vel
1758 implicit real*8 (a-h,o-z)
1759 include 'DIMENSIONS'
1760 include 'COMMON.CONTROL'
1761 include 'COMMON.VAR'
1764 include 'COMMON.LANGEVIN'
1766 include 'COMMON.LANGEVIN.lang0'
1768 include 'COMMON.CHAIN'
1769 include 'COMMON.DERIV'
1770 include 'COMMON.GEO'
1771 include 'COMMON.LOCAL'
1772 include 'COMMON.INTERACT'
1773 include 'COMMON.IOUNITS'
1774 include 'COMMON.NAMES'
1775 include 'COMMON.TIME1'
1776 double precision xv,sigv,lowb,highb
1777 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1778 c First generate velocities in the eigenspace of the G matrix
1779 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1786 sigv=dsqrt((Rb*t_bath)/geigen(i))
1789 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1790 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1791 c & " d_t_work_new",d_t_work_new(ii)
1800 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1803 c write (iout,*) "Ek from eigenvectors",Ek1
1805 c Transform velocities to UNRES coordinate space
1811 d_t_work(ind)=d_t_work(ind)
1812 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1814 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1818 c Transfer to the d_t vector
1820 d_t(j,0)=d_t_work(j)
1826 d_t(j,i)=d_t_work(ind)
1830 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1833 d_t(j,i+nres)=d_t_work(ind)
1838 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1839 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1844 c-----------------------------------------------------------
1845 subroutine sd_verlet_p_setup
1846 c Sets up the parameters of stochastic Verlet algorithm
1847 implicit real*8 (a-h,o-z)
1848 include 'DIMENSIONS'
1852 include 'COMMON.CONTROL'
1853 include 'COMMON.VAR'
1856 include 'COMMON.LANGEVIN'
1858 include 'COMMON.LANGEVIN.lang0'
1860 include 'COMMON.CHAIN'
1861 include 'COMMON.DERIV'
1862 include 'COMMON.GEO'
1863 include 'COMMON.LOCAL'
1864 include 'COMMON.INTERACT'
1865 include 'COMMON.IOUNITS'
1866 include 'COMMON.NAMES'
1867 include 'COMMON.TIME1'
1868 double precision emgdt(MAXRES6),
1869 & pterm,vterm,rho,rhoc,vsig,
1870 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1871 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1872 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1873 logical lprn /.false./
1874 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1875 double precision ktm
1882 c AL 8/17/04 Code adapted from tinker
1884 c Get the frictional and random terms for stochastic dynamics in the
1885 c eigenspace of mass-scaled UNRES friction matrix
1888 gdt = fricgam(i) * d_time
1890 c Stochastic dynamics reduces to simple MD for zero friction
1892 if (gdt .le. zero) then
1893 pfric_vec(i) = 1.0d0
1894 vfric_vec(i) = d_time
1895 afric_vec(i) = 0.5d0 * d_time * d_time
1896 prand_vec(i) = 0.0d0
1897 vrand_vec1(i) = 0.0d0
1898 vrand_vec2(i) = 0.0d0
1900 c Analytical expressions when friction coefficient is large
1903 if (gdt .ge. gdt_radius) then
1906 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1907 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1908 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1909 vterm = 1.0d0 - egdt**2
1910 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1912 c Use series expansions when friction coefficient is small
1923 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1924 & - gdt5/120.0d0 + gdt6/720.0d0
1925 & - gdt7/5040.0d0 + gdt8/40320.0d0
1926 & - gdt9/362880.0d0) / fricgam(i)**2
1927 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1928 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1929 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1930 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1931 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1932 & + 127.0d0*gdt9/90720.0d0
1933 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1934 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1935 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1936 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1937 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1938 & - 17.0d0*gdt2/1280.0d0
1939 & + 17.0d0*gdt3/6144.0d0
1940 & + 40967.0d0*gdt4/34406400.0d0
1941 & - 57203.0d0*gdt5/275251200.0d0
1942 & - 1429487.0d0*gdt6/13212057600.0d0)
1945 c Compute the scaling factors of random terms for the nonzero friction case
1947 ktm = 0.5d0*d_time/fricgam(i)
1948 psig = dsqrt(ktm*pterm) / fricgam(i)
1949 vsig = dsqrt(ktm*vterm)
1950 rhoc = dsqrt(1.0d0 - rho*rho)
1952 vrand_vec1(i) = vsig * rho
1953 vrand_vec2(i) = vsig * rhoc
1958 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1961 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1962 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1966 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1969 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1970 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1971 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1972 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1973 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1974 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1977 t_sdsetup=t_sdsetup+MPI_Wtime()
1979 t_sdsetup=t_sdsetup+tcpu()-tt0
1983 c-------------------------------------------------------------
1984 subroutine eigtransf1(n,ndim,ab,d,c)
1987 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
1993 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
1999 c-------------------------------------------------------------
2000 subroutine eigtransf(n,ndim,a,b,d,c)
2003 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2009 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2015 c-------------------------------------------------------------
2016 subroutine sd_verlet1
2017 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2018 implicit real*8 (a-h,o-z)
2019 include 'DIMENSIONS'
2020 include 'COMMON.CONTROL'
2021 include 'COMMON.VAR'
2024 include 'COMMON.LANGEVIN'
2026 include 'COMMON.LANGEVIN.lang0'
2028 include 'COMMON.CHAIN'
2029 include 'COMMON.DERIV'
2030 include 'COMMON.GEO'
2031 include 'COMMON.LOCAL'
2032 include 'COMMON.INTERACT'
2033 include 'COMMON.IOUNITS'
2034 include 'COMMON.NAMES'
2035 double precision stochforcvec(MAXRES6)
2036 common /stochcalc/ stochforcvec
2037 logical lprn /.false./
2039 c write (iout,*) "dc_old"
2041 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2042 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2045 dc_work(j)=dc_old(j,0)
2046 d_t_work(j)=d_t_old(j,0)
2047 d_a_work(j)=d_a_old(j,0)
2052 dc_work(ind+j)=dc_old(j,i)
2053 d_t_work(ind+j)=d_t_old(j,i)
2054 d_a_work(ind+j)=d_a_old(j,i)
2059 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2061 dc_work(ind+j)=dc_old(j,i+nres)
2062 d_t_work(ind+j)=d_t_old(j,i+nres)
2063 d_a_work(ind+j)=d_a_old(j,i+nres)
2071 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2075 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2076 & vfric_mat(i,j),afric_mat(i,j),
2077 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2085 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2086 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2087 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2088 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2090 d_t_work_new(i)=ddt1+0.5d0*ddt2
2091 d_t_work(i)=ddt1+ddt2
2096 d_t(j,0)=d_t_work(j)
2101 dc(j,i)=dc_work(ind+j)
2102 d_t(j,i)=d_t_work(ind+j)
2107 if (itype(i).ne.10) then
2110 dc(j,inres)=dc_work(ind+j)
2111 d_t(j,inres)=d_t_work(ind+j)
2118 c--------------------------------------------------------------------------
2119 subroutine sd_verlet2
2120 c Calculating the adjusted velocities for accelerations
2121 implicit real*8 (a-h,o-z)
2122 include 'DIMENSIONS'
2123 include 'COMMON.CONTROL'
2124 include 'COMMON.VAR'
2127 include 'COMMON.LANGEVIN'
2129 include 'COMMON.LANGEVIN.lang0'
2131 include 'COMMON.CHAIN'
2132 include 'COMMON.DERIV'
2133 include 'COMMON.GEO'
2134 include 'COMMON.LOCAL'
2135 include 'COMMON.INTERACT'
2136 include 'COMMON.IOUNITS'
2137 include 'COMMON.NAMES'
2138 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2139 common /stochcalc/ stochforcvec
2141 c Compute the stochastic forces which contribute to velocity change
2143 call stochastic_force(stochforcvecV)
2150 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2151 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2152 & vrand_mat2(i,j)*stochforcvecV(j)
2154 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2158 d_t(j,0)=d_t_work(j)
2163 d_t(j,i)=d_t_work(ind+j)
2168 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2171 d_t(j,inres)=d_t_work(ind+j)
2178 c-----------------------------------------------------------
2179 subroutine sd_verlet_ciccotti_setup
2180 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2182 implicit real*8 (a-h,o-z)
2183 include 'DIMENSIONS'
2187 include 'COMMON.CONTROL'
2188 include 'COMMON.VAR'
2191 include 'COMMON.LANGEVIN'
2193 include 'COMMON.LANGEVIN.lang0'
2195 include 'COMMON.CHAIN'
2196 include 'COMMON.DERIV'
2197 include 'COMMON.GEO'
2198 include 'COMMON.LOCAL'
2199 include 'COMMON.INTERACT'
2200 include 'COMMON.IOUNITS'
2201 include 'COMMON.NAMES'
2202 include 'COMMON.TIME1'
2203 double precision emgdt(MAXRES6),
2204 & pterm,vterm,rho,rhoc,vsig,
2205 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2206 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2207 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2208 logical lprn /.false./
2209 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2210 double precision ktm
2217 c AL 8/17/04 Code adapted from tinker
2219 c Get the frictional and random terms for stochastic dynamics in the
2220 c eigenspace of mass-scaled UNRES friction matrix
2223 write (iout,*) "i",i," fricgam",fricgam(i)
2224 gdt = fricgam(i) * d_time
2226 c Stochastic dynamics reduces to simple MD for zero friction
2228 if (gdt .le. zero) then
2229 pfric_vec(i) = 1.0d0
2230 vfric_vec(i) = d_time
2231 afric_vec(i) = 0.5d0*d_time*d_time
2232 prand_vec(i) = afric_vec(i)
2233 vrand_vec2(i) = vfric_vec(i)
2235 c Analytical expressions when friction coefficient is large
2240 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2241 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2242 prand_vec(i) = afric_vec(i)
2243 vrand_vec2(i) = vfric_vec(i)
2245 c Compute the scaling factors of random terms for the nonzero friction case
2247 c ktm = 0.5d0*d_time/fricgam(i)
2248 c psig = dsqrt(ktm*pterm) / fricgam(i)
2249 c vsig = dsqrt(ktm*vterm)
2250 c prand_vec(i) = psig*afric_vec(i)
2251 c vrand_vec2(i) = vsig*vfric_vec(i)
2256 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2259 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2260 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2264 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2266 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2267 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2268 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2269 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2270 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2272 t_sdsetup=t_sdsetup+MPI_Wtime()
2274 t_sdsetup=t_sdsetup+tcpu()-tt0
2278 c-------------------------------------------------------------
2279 subroutine sd_verlet1_ciccotti
2280 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2281 implicit real*8 (a-h,o-z)
2282 include 'DIMENSIONS'
2286 include 'COMMON.CONTROL'
2287 include 'COMMON.VAR'
2290 include 'COMMON.LANGEVIN'
2292 include 'COMMON.LANGEVIN.lang0'
2294 include 'COMMON.CHAIN'
2295 include 'COMMON.DERIV'
2296 include 'COMMON.GEO'
2297 include 'COMMON.LOCAL'
2298 include 'COMMON.INTERACT'
2299 include 'COMMON.IOUNITS'
2300 include 'COMMON.NAMES'
2301 double precision stochforcvec(MAXRES6)
2302 common /stochcalc/ stochforcvec
2303 logical lprn /.false./
2305 c write (iout,*) "dc_old"
2307 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2308 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2311 dc_work(j)=dc_old(j,0)
2312 d_t_work(j)=d_t_old(j,0)
2313 d_a_work(j)=d_a_old(j,0)
2318 dc_work(ind+j)=dc_old(j,i)
2319 d_t_work(ind+j)=d_t_old(j,i)
2320 d_a_work(ind+j)=d_a_old(j,i)
2325 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2327 dc_work(ind+j)=dc_old(j,i+nres)
2328 d_t_work(ind+j)=d_t_old(j,i+nres)
2329 d_a_work(ind+j)=d_a_old(j,i+nres)
2338 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2342 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2343 & vfric_mat(i,j),afric_mat(i,j),
2344 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2352 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2353 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2354 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2355 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2357 d_t_work_new(i)=ddt1+0.5d0*ddt2
2358 d_t_work(i)=ddt1+ddt2
2363 d_t(j,0)=d_t_work(j)
2368 dc(j,i)=dc_work(ind+j)
2369 d_t(j,i)=d_t_work(ind+j)
2374 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2377 dc(j,inres)=dc_work(ind+j)
2378 d_t(j,inres)=d_t_work(ind+j)
2385 c--------------------------------------------------------------------------
2386 subroutine sd_verlet2_ciccotti
2387 c Calculating the adjusted velocities for accelerations
2388 implicit real*8 (a-h,o-z)
2389 include 'DIMENSIONS'
2390 include 'COMMON.CONTROL'
2391 include 'COMMON.VAR'
2394 include 'COMMON.LANGEVIN'
2396 include 'COMMON.LANGEVIN.lang0'
2398 include 'COMMON.CHAIN'
2399 include 'COMMON.DERIV'
2400 include 'COMMON.GEO'
2401 include 'COMMON.LOCAL'
2402 include 'COMMON.INTERACT'
2403 include 'COMMON.IOUNITS'
2404 include 'COMMON.NAMES'
2405 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2406 common /stochcalc/ stochforcvec
2408 c Compute the stochastic forces which contribute to velocity change
2410 call stochastic_force(stochforcvecV)
2417 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2418 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2419 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2421 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2425 d_t(j,0)=d_t_work(j)
2430 d_t(j,i)=d_t_work(ind+j)
2435 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2438 d_t(j,inres)=d_t_work(ind+j)