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 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.ntyp1) 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.ntyp1) 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.ntyp1) 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(iabs(itype(i)))
1437 & *dsqrt(gamsc(iabs(itype(i))))
1440 c Open the pdb file for snapshotshots
1443 if (ilen(tmpdir).gt.0)
1444 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1445 & liczba(:ilen(liczba))//".pdb")
1447 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1451 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1452 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1453 & liczba(:ilen(liczba))//".x")
1454 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1457 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1458 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1459 & liczba(:ilen(liczba))//".cx")
1460 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1466 if (ilen(tmpdir).gt.0)
1467 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1468 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1470 if (ilen(tmpdir).gt.0)
1471 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1472 cartname=prefix(:ilen(prefix))//"_MD.cx"
1476 write (qstr,'(256(1h ))')
1479 iq = qinfrag(i,iset)*10
1480 iw = wfrag(i,iset)/100
1482 if(me.eq.king.or..not.out1file)
1483 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1484 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1489 iq = qinpair(i,iset)*10
1490 iw = wpair(i,iset)/100
1492 if(me.eq.king.or..not.out1file)
1493 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1494 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1498 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1500 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1502 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1504 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1508 if (restart1file) then
1510 & inquire(file=mremd_rst_name,exist=file_exist)
1511 write (*,*) me," Before broadcast: file_exist",file_exist
1512 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1514 write (*,*) me," After broadcast: file_exist",file_exist
1515 c inquire(file=mremd_rst_name,exist=file_exist)
1516 if(me.eq.king.or..not.out1file)
1517 & write(iout,*) "Initial state read by master and distributed"
1519 if (ilen(tmpdir).gt.0)
1520 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1521 & //liczba(:ilen(liczba))//'.rst')
1522 inquire(file=rest2name,exist=file_exist)
1525 if(.not.restart1file) then
1526 if(me.eq.king.or..not.out1file)
1527 & write(iout,*) "Initial state will be read from file ",
1528 & rest2name(:ilen(rest2name))
1531 call rescale_weights(t_bath)
1533 if(me.eq.king.or..not.out1file)then
1534 if (restart1file) then
1535 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1538 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1541 write(iout,*) "Initial velocities randomly generated"
1547 c Generate initial velocities
1548 if(me.eq.king.or..not.out1file)
1549 & write(iout,*) "Initial velocities randomly generated"
1553 c rest2name = prefix(:ilen(prefix))//'.rst'
1554 if(me.eq.king.or..not.out1file)then
1555 write (iout,*) "Initial velocities"
1557 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1558 & (d_t(j,i+nres),j=1,3)
1560 c Zeroing the total angular momentum of the system
1561 write(iout,*) "Calling the zero-angular
1562 & momentum subroutine"
1565 c Getting the potential energy and forces and velocities and accelerations
1567 c write (iout,*) "velocity of the center of the mass:"
1568 c write (iout,*) (vcm(j),j=1,3)
1570 d_t(j,0)=d_t(j,0)-vcm(j)
1572 c Removing the velocity of the center of mass
1574 if(me.eq.king.or..not.out1file)then
1575 write (iout,*) "vcm right after adjustment:"
1576 write (iout,*) (vcm(j),j=1,3)
1580 if(iranconf.ne.0) then
1582 print *, 'Calling OVERLAP_SC'
1583 call overlap_sc(fail)
1587 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1588 print *,'SC_move',nft_sc,etot
1589 if(me.eq.king.or..not.out1file)
1590 & write(iout,*) 'SC_move',nft_sc,etot
1594 print *, 'Calling MINIM_DC'
1595 call minim_dc(etot,iretcode,nfun)
1597 call geom_to_var(nvar,varia)
1598 print *,'Calling MINIMIZE.'
1599 call minimize(etot,varia,iretcode,nfun)
1600 call var_to_geom(nvar,varia)
1602 if(me.eq.king.or..not.out1file)
1603 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1606 call chainbuild_cart
1609 call verlet_bath(EK)
1611 kinetic_T=2.0d0/(dimen3*Rb)*EK
1612 if(me.eq.king.or..not.out1file)then
1622 call etotal(potEcomp)
1623 if (large) call enerprint(potEcomp)
1626 t_etotal=t_etotal+MPI_Wtime()-tt0
1628 t_etotal=t_etotal+tcpu()-tt0
1635 if (amax*d_time .gt. dvmax) then
1636 d_time=d_time*dvmax/amax
1637 if(me.eq.king.or..not.out1file) write (iout,*)
1638 & "Time step reduced to",d_time,
1639 & " because of too large initial acceleration."
1641 if(me.eq.king.or..not.out1file)then
1642 write(iout,*) "Potential energy and its components"
1643 call enerprint(potEcomp)
1644 c write(iout,*) (potEcomp(i),i=0,n_ene)
1646 potE=potEcomp(0)-potEcomp(20)
1649 if (ntwe.ne.0) call statout(itime)
1650 if(me.eq.king.or..not.out1file)
1651 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1652 & " Kinetic energy",EK," Potential energy",potE,
1653 & " Total energy",totE," Maximum acceleration ",
1656 write (iout,*) "Initial coordinates"
1658 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1659 & (c(j,i+nres),j=1,3)
1661 write (iout,*) "Initial dC"
1663 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1664 & (dc(j,i+nres),j=1,3)
1666 write (iout,*) "Initial velocities"
1667 write (iout,"(13x,' backbone ',23x,' side chain')")
1669 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1670 & (d_t(j,i+nres),j=1,3)
1672 write (iout,*) "Initial accelerations"
1674 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1675 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1676 & (d_a(j,i+nres),j=1,3)
1682 d_t_old(j,i)=d_t(j,i)
1683 d_a_old(j,i)=d_a(j,i)
1685 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1694 call etotal_short(energia_short)
1695 if (large) call enerprint(potEcomp)
1698 t_eshort=t_eshort+MPI_Wtime()-tt0
1700 t_eshort=t_eshort+tcpu()-tt0
1705 if(.not.out1file .and. large) then
1706 write (iout,*) "energia_long",energia_long(0),
1707 & " energia_short",energia_short(0),
1708 & " total",energia_long(0)+energia_short(0)
1709 write (iout,*) "Initial fast-force accelerations"
1711 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1712 & (d_a(j,i+nres),j=1,3)
1715 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1718 d_a_short(j,i)=d_a(j,i)
1727 call etotal_long(energia_long)
1728 if (large) call enerprint(potEcomp)
1731 t_elong=t_elong+MPI_Wtime()-tt0
1733 t_elong=t_elong+tcpu()-tt0
1738 if(.not.out1file .and. large) then
1739 write (iout,*) "energia_long",energia_long(0)
1740 write (iout,*) "Initial slow-force accelerations"
1742 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1743 & (d_a(j,i+nres),j=1,3)
1747 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1749 t_enegrad=t_enegrad+tcpu()-tt0
1754 c-----------------------------------------------------------
1755 subroutine random_vel
1756 implicit real*8 (a-h,o-z)
1757 include 'DIMENSIONS'
1758 include 'COMMON.CONTROL'
1759 include 'COMMON.VAR'
1762 include 'COMMON.LANGEVIN'
1764 include 'COMMON.LANGEVIN.lang0'
1766 include 'COMMON.CHAIN'
1767 include 'COMMON.DERIV'
1768 include 'COMMON.GEO'
1769 include 'COMMON.LOCAL'
1770 include 'COMMON.INTERACT'
1771 include 'COMMON.IOUNITS'
1772 include 'COMMON.NAMES'
1773 include 'COMMON.TIME1'
1774 double precision xv,sigv,lowb,highb
1775 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1776 c First generate velocities in the eigenspace of the G matrix
1777 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1784 sigv=dsqrt((Rb*t_bath)/geigen(i))
1787 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1788 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1789 c & " d_t_work_new",d_t_work_new(ii)
1798 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1801 c write (iout,*) "Ek from eigenvectors",Ek1
1803 c Transform velocities to UNRES coordinate space
1809 d_t_work(ind)=d_t_work(ind)
1810 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1812 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1816 c Transfer to the d_t vector
1818 d_t(j,0)=d_t_work(j)
1824 d_t(j,i)=d_t_work(ind)
1828 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1831 d_t(j,i+nres)=d_t_work(ind)
1836 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1837 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1842 c-----------------------------------------------------------
1843 subroutine sd_verlet_p_setup
1844 c Sets up the parameters of stochastic Verlet algorithm
1845 implicit real*8 (a-h,o-z)
1846 include 'DIMENSIONS'
1850 include 'COMMON.CONTROL'
1851 include 'COMMON.VAR'
1854 include 'COMMON.LANGEVIN'
1856 include 'COMMON.LANGEVIN.lang0'
1858 include 'COMMON.CHAIN'
1859 include 'COMMON.DERIV'
1860 include 'COMMON.GEO'
1861 include 'COMMON.LOCAL'
1862 include 'COMMON.INTERACT'
1863 include 'COMMON.IOUNITS'
1864 include 'COMMON.NAMES'
1865 include 'COMMON.TIME1'
1866 double precision emgdt(MAXRES6),
1867 & pterm,vterm,rho,rhoc,vsig,
1868 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1869 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1870 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1871 logical lprn /.false./
1872 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1873 double precision ktm
1880 c AL 8/17/04 Code adapted from tinker
1882 c Get the frictional and random terms for stochastic dynamics in the
1883 c eigenspace of mass-scaled UNRES friction matrix
1886 gdt = fricgam(i) * d_time
1888 c Stochastic dynamics reduces to simple MD for zero friction
1890 if (gdt .le. zero) then
1891 pfric_vec(i) = 1.0d0
1892 vfric_vec(i) = d_time
1893 afric_vec(i) = 0.5d0 * d_time * d_time
1894 prand_vec(i) = 0.0d0
1895 vrand_vec1(i) = 0.0d0
1896 vrand_vec2(i) = 0.0d0
1898 c Analytical expressions when friction coefficient is large
1901 if (gdt .ge. gdt_radius) then
1904 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1905 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1906 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1907 vterm = 1.0d0 - egdt**2
1908 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1910 c Use series expansions when friction coefficient is small
1921 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1922 & - gdt5/120.0d0 + gdt6/720.0d0
1923 & - gdt7/5040.0d0 + gdt8/40320.0d0
1924 & - gdt9/362880.0d0) / fricgam(i)**2
1925 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1926 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1927 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1928 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1929 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1930 & + 127.0d0*gdt9/90720.0d0
1931 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1932 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1933 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1934 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1935 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1936 & - 17.0d0*gdt2/1280.0d0
1937 & + 17.0d0*gdt3/6144.0d0
1938 & + 40967.0d0*gdt4/34406400.0d0
1939 & - 57203.0d0*gdt5/275251200.0d0
1940 & - 1429487.0d0*gdt6/13212057600.0d0)
1943 c Compute the scaling factors of random terms for the nonzero friction case
1945 ktm = 0.5d0*d_time/fricgam(i)
1946 psig = dsqrt(ktm*pterm) / fricgam(i)
1947 vsig = dsqrt(ktm*vterm)
1948 rhoc = dsqrt(1.0d0 - rho*rho)
1950 vrand_vec1(i) = vsig * rho
1951 vrand_vec2(i) = vsig * rhoc
1956 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1959 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1960 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1964 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1967 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1968 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1969 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1970 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1971 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1972 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1975 t_sdsetup=t_sdsetup+MPI_Wtime()
1977 t_sdsetup=t_sdsetup+tcpu()-tt0
1981 c-------------------------------------------------------------
1982 subroutine eigtransf1(n,ndim,ab,d,c)
1985 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
1991 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
1997 c-------------------------------------------------------------
1998 subroutine eigtransf(n,ndim,a,b,d,c)
2001 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2007 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2013 c-------------------------------------------------------------
2014 subroutine sd_verlet1
2015 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2016 implicit real*8 (a-h,o-z)
2017 include 'DIMENSIONS'
2018 include 'COMMON.CONTROL'
2019 include 'COMMON.VAR'
2022 include 'COMMON.LANGEVIN'
2024 include 'COMMON.LANGEVIN.lang0'
2026 include 'COMMON.CHAIN'
2027 include 'COMMON.DERIV'
2028 include 'COMMON.GEO'
2029 include 'COMMON.LOCAL'
2030 include 'COMMON.INTERACT'
2031 include 'COMMON.IOUNITS'
2032 include 'COMMON.NAMES'
2033 double precision stochforcvec(MAXRES6)
2034 common /stochcalc/ stochforcvec
2035 logical lprn /.false./
2037 c write (iout,*) "dc_old"
2039 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2040 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2043 dc_work(j)=dc_old(j,0)
2044 d_t_work(j)=d_t_old(j,0)
2045 d_a_work(j)=d_a_old(j,0)
2050 dc_work(ind+j)=dc_old(j,i)
2051 d_t_work(ind+j)=d_t_old(j,i)
2052 d_a_work(ind+j)=d_a_old(j,i)
2057 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2059 dc_work(ind+j)=dc_old(j,i+nres)
2060 d_t_work(ind+j)=d_t_old(j,i+nres)
2061 d_a_work(ind+j)=d_a_old(j,i+nres)
2069 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2073 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2074 & vfric_mat(i,j),afric_mat(i,j),
2075 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2083 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2084 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2085 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2086 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2088 d_t_work_new(i)=ddt1+0.5d0*ddt2
2089 d_t_work(i)=ddt1+ddt2
2094 d_t(j,0)=d_t_work(j)
2099 dc(j,i)=dc_work(ind+j)
2100 d_t(j,i)=d_t_work(ind+j)
2105 if (itype(i).ne.10) then
2108 dc(j,inres)=dc_work(ind+j)
2109 d_t(j,inres)=d_t_work(ind+j)
2116 c--------------------------------------------------------------------------
2117 subroutine sd_verlet2
2118 c Calculating the adjusted velocities for accelerations
2119 implicit real*8 (a-h,o-z)
2120 include 'DIMENSIONS'
2121 include 'COMMON.CONTROL'
2122 include 'COMMON.VAR'
2125 include 'COMMON.LANGEVIN'
2127 include 'COMMON.LANGEVIN.lang0'
2129 include 'COMMON.CHAIN'
2130 include 'COMMON.DERIV'
2131 include 'COMMON.GEO'
2132 include 'COMMON.LOCAL'
2133 include 'COMMON.INTERACT'
2134 include 'COMMON.IOUNITS'
2135 include 'COMMON.NAMES'
2136 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2137 common /stochcalc/ stochforcvec
2139 c Compute the stochastic forces which contribute to velocity change
2141 call stochastic_force(stochforcvecV)
2148 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2149 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2150 & vrand_mat2(i,j)*stochforcvecV(j)
2152 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2156 d_t(j,0)=d_t_work(j)
2161 d_t(j,i)=d_t_work(ind+j)
2166 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2169 d_t(j,inres)=d_t_work(ind+j)
2176 c-----------------------------------------------------------
2177 subroutine sd_verlet_ciccotti_setup
2178 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2180 implicit real*8 (a-h,o-z)
2181 include 'DIMENSIONS'
2185 include 'COMMON.CONTROL'
2186 include 'COMMON.VAR'
2189 include 'COMMON.LANGEVIN'
2191 include 'COMMON.LANGEVIN.lang0'
2193 include 'COMMON.CHAIN'
2194 include 'COMMON.DERIV'
2195 include 'COMMON.GEO'
2196 include 'COMMON.LOCAL'
2197 include 'COMMON.INTERACT'
2198 include 'COMMON.IOUNITS'
2199 include 'COMMON.NAMES'
2200 include 'COMMON.TIME1'
2201 double precision emgdt(MAXRES6),
2202 & pterm,vterm,rho,rhoc,vsig,
2203 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2204 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2205 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2206 logical lprn /.false./
2207 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2208 double precision ktm
2215 c AL 8/17/04 Code adapted from tinker
2217 c Get the frictional and random terms for stochastic dynamics in the
2218 c eigenspace of mass-scaled UNRES friction matrix
2221 write (iout,*) "i",i," fricgam",fricgam(i)
2222 gdt = fricgam(i) * d_time
2224 c Stochastic dynamics reduces to simple MD for zero friction
2226 if (gdt .le. zero) then
2227 pfric_vec(i) = 1.0d0
2228 vfric_vec(i) = d_time
2229 afric_vec(i) = 0.5d0*d_time*d_time
2230 prand_vec(i) = afric_vec(i)
2231 vrand_vec2(i) = vfric_vec(i)
2233 c Analytical expressions when friction coefficient is large
2238 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2239 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2240 prand_vec(i) = afric_vec(i)
2241 vrand_vec2(i) = vfric_vec(i)
2243 c Compute the scaling factors of random terms for the nonzero friction case
2245 c ktm = 0.5d0*d_time/fricgam(i)
2246 c psig = dsqrt(ktm*pterm) / fricgam(i)
2247 c vsig = dsqrt(ktm*vterm)
2248 c prand_vec(i) = psig*afric_vec(i)
2249 c vrand_vec2(i) = vsig*vfric_vec(i)
2254 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2257 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2258 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2262 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2264 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2265 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2266 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2267 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2268 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2270 t_sdsetup=t_sdsetup+MPI_Wtime()
2272 t_sdsetup=t_sdsetup+tcpu()-tt0
2276 c-------------------------------------------------------------
2277 subroutine sd_verlet1_ciccotti
2278 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2279 implicit real*8 (a-h,o-z)
2280 include 'DIMENSIONS'
2284 include 'COMMON.CONTROL'
2285 include 'COMMON.VAR'
2288 include 'COMMON.LANGEVIN'
2290 include 'COMMON.LANGEVIN.lang0'
2292 include 'COMMON.CHAIN'
2293 include 'COMMON.DERIV'
2294 include 'COMMON.GEO'
2295 include 'COMMON.LOCAL'
2296 include 'COMMON.INTERACT'
2297 include 'COMMON.IOUNITS'
2298 include 'COMMON.NAMES'
2299 double precision stochforcvec(MAXRES6)
2300 common /stochcalc/ stochforcvec
2301 logical lprn /.false./
2303 c write (iout,*) "dc_old"
2305 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2306 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2309 dc_work(j)=dc_old(j,0)
2310 d_t_work(j)=d_t_old(j,0)
2311 d_a_work(j)=d_a_old(j,0)
2316 dc_work(ind+j)=dc_old(j,i)
2317 d_t_work(ind+j)=d_t_old(j,i)
2318 d_a_work(ind+j)=d_a_old(j,i)
2323 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2325 dc_work(ind+j)=dc_old(j,i+nres)
2326 d_t_work(ind+j)=d_t_old(j,i+nres)
2327 d_a_work(ind+j)=d_a_old(j,i+nres)
2336 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2340 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2341 & vfric_mat(i,j),afric_mat(i,j),
2342 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2350 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2351 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2352 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2353 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2355 d_t_work_new(i)=ddt1+0.5d0*ddt2
2356 d_t_work(i)=ddt1+ddt2
2361 d_t(j,0)=d_t_work(j)
2366 dc(j,i)=dc_work(ind+j)
2367 d_t(j,i)=d_t_work(ind+j)
2372 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2375 dc(j,inres)=dc_work(ind+j)
2376 d_t(j,inres)=d_t_work(ind+j)
2383 c--------------------------------------------------------------------------
2384 subroutine sd_verlet2_ciccotti
2385 c Calculating the adjusted velocities for accelerations
2386 implicit real*8 (a-h,o-z)
2387 include 'DIMENSIONS'
2388 include 'COMMON.CONTROL'
2389 include 'COMMON.VAR'
2392 include 'COMMON.LANGEVIN'
2394 include 'COMMON.LANGEVIN.lang0'
2396 include 'COMMON.CHAIN'
2397 include 'COMMON.DERIV'
2398 include 'COMMON.GEO'
2399 include 'COMMON.LOCAL'
2400 include 'COMMON.INTERACT'
2401 include 'COMMON.IOUNITS'
2402 include 'COMMON.NAMES'
2403 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2404 common /stochcalc/ stochforcvec
2406 c Compute the stochastic forces which contribute to velocity change
2408 call stochastic_force(stochforcvecV)
2415 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2416 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2417 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2419 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2423 d_t(j,0)=d_t_work(j)
2428 d_t(j,i)=d_t_work(ind+j)
2433 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2436 d_t(j,inres)=d_t_work(ind+j)