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 double precision cm(3),L(3),vcm(3)
30 double precision v_work(maxres6),v_transf(maxres6)
39 if (ilen(tmpdir).gt.0)
40 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
41 & //liczba(:ilen(liczba))//'.rst')
43 if (ilen(tmpdir).gt.0)
44 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
51 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
57 c Determine the inverse of the inertia matrix.
58 call setup_MD_matrices
62 t_MDsetup = MPI_Wtime()-tt0
64 t_MDsetup = tcpu()-tt0
67 c Entering the MD loop
73 if (lang.eq.2 .or. lang.eq.3) then
77 call sd_verlet_p_setup
79 call sd_verlet_ciccotti_setup
83 pfric0_mat(i,j,0)=pfric_mat(i,j)
84 afric0_mat(i,j,0)=afric_mat(i,j)
85 vfric0_mat(i,j,0)=vfric_mat(i,j)
86 prand0_mat(i,j,0)=prand_mat(i,j)
87 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
88 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
97 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
99 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
103 else if (lang.eq.1 .or. lang.eq.4) then
107 t_langsetup=MPI_Wtime()-tt0
110 t_langsetup=tcpu()-tt0
113 do itime=1,n_timestep
115 if (lang.gt.0 .and. surfarea .and.
116 & mod(itime,reset_fricmat).eq.0) then
117 if (lang.eq.2 .or. lang.eq.3) then
121 call sd_verlet_p_setup
123 call sd_verlet_ciccotti_setup
127 pfric0_mat(i,j,0)=pfric_mat(i,j)
128 afric0_mat(i,j,0)=afric_mat(i,j)
129 vfric0_mat(i,j,0)=vfric_mat(i,j)
130 prand0_mat(i,j,0)=prand_mat(i,j)
131 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
132 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
137 flag_stoch(i)=.false.
140 else if (lang.eq.1 .or. lang.eq.4) then
143 write (iout,'(a,i10)')
144 & "Friction matrix reset based on surface area, itime",itime
146 if (reset_vel .and. tbf .and. lang.eq.0
147 & .and. mod(itime,count_reset_vel).eq.0) then
149 write(iout,'(a,f20.2)')
150 & "Velocities reset to random values, time",totT
153 d_t_old(j,i)=d_t(j,i)
157 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
161 d_t(j,0)=d_t(j,0)-vcm(j)
164 kinetic_T=2.0d0/(dimen3*Rb)*EK
165 scalfac=dsqrt(T_bath/kinetic_T)
166 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
169 d_t_old(j,i)=scalfac*d_t(j,i)
175 c Time-reversible RESPA algorithm
176 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
177 call RESPA_step(itime)
179 c Variable time step algorithm.
180 call velverlet_step(itime)
184 call brown_step(itime)
186 print *,"Brown dynamics not here!"
188 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
194 if (mod(itime,ntwe).eq.0) call statout(itime)
207 if (itype(i).ne.10) then
210 v_work(ind)=d_t(j,i+nres)
215 write (66,'(80f10.5)')
216 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
220 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
222 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
224 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
227 if (mod(itime,ntwx).eq.0) then
228 write (tytul,'("time",f8.2)') totT
230 call pdbout(potE,tytul,ipdb)
235 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
236 open(irest2,file=rest2name,status='unknown')
237 write(irest2,*) totT,EK,potE,totE,t_bath
239 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
242 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
253 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
255 & 'MD calculations setup:',t_MDsetup,
256 & 'Energy & gradient evaluation:',t_enegrad,
257 & 'Stochastic MD setup:',t_langsetup,
258 & 'Stochastic MD step setup:',t_sdsetup,
260 write (iout,'(/28(1h=),a25,27(1h=))')
261 & ' End of MD calculation '
264 c-------------------------------------------------------------------------------
265 subroutine velverlet_step(itime)
266 c-------------------------------------------------------------------------------
267 c Perform a single velocity Verlet step; the time step can be rescaled if
268 c increments in accelerations exceed the threshold
269 c-------------------------------------------------------------------------------
270 implicit real*8 (a-h,o-z)
274 integer ierror,ierrcode
276 include 'COMMON.SETUP'
277 include 'COMMON.CONTROL'
281 include 'COMMON.LANGEVIN'
283 include 'COMMON.LANGEVIN.lang0'
285 include 'COMMON.CHAIN'
286 include 'COMMON.DERIV'
288 include 'COMMON.LOCAL'
289 include 'COMMON.INTERACT'
290 include 'COMMON.IOUNITS'
291 include 'COMMON.NAMES'
292 include 'COMMON.TIME1'
293 include 'COMMON.MUCA'
294 double precision vcm(3),incr(3)
295 double precision cm(3),L(3)
296 integer ilen,count,rstcount
299 integer maxcount_scale /20/
301 double precision stochforcvec(MAXRES6)
302 common /stochcalc/ stochforcvec
310 else if (lang.eq.2 .or. lang.eq.3) then
312 call stochastic_force(stochforcvec)
315 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
317 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
324 icount_scale=icount_scale+1
325 if (icount_scale.gt.maxcount_scale) then
327 & "ERROR: too many attempts at scaling down the time step. ",
328 & "amax=",amax,"epdrift=",epdrift,
329 & "damax=",damax,"edriftmax=",edriftmax,
333 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
337 c First step of the velocity Verlet algorithm
342 else if (lang.eq.3) then
344 call sd_verlet1_ciccotti
346 else if (lang.eq.1) then
351 c Build the chain from the newly calculated coordinates
353 if (rattle) call rattle1
355 if (large.and. mod(itime,ntwe).eq.0) then
356 write (iout,*) "Cartesian and internal coordinates: step 1"
361 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
362 & (dc(j,i+nres),j=1,3)
364 write (iout,*) "Accelerations"
366 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
367 & (d_a(j,i+nres),j=1,3)
369 write (iout,*) "Velocities, step 1"
371 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
372 & (d_t(j,i+nres),j=1,3)
381 c Calculate energy and forces
383 call etotal(potEcomp)
384 potE=potEcomp(0)-potEcomp(20)
386 c Get the new accelerations
389 t_enegrad=t_enegrad+MPI_Wtime()-tt0
391 t_enegrad=t_enegrad+tcpu()-tt0
393 c Determine maximum acceleration and scale down the timestep if needed
395 amax=amax/(itime_scal+1)**2
396 call predict_edrift(epdrift)
397 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
398 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
400 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
402 itime_scal=itime_scal+ifac_time
403 c fac_time=dmin1(damax/amax,0.5d0)
404 fac_time=0.5d0**ifac_time
405 d_time=d_time*fac_time
406 if (lang.eq.2 .or. lang.eq.3) then
408 c write (iout,*) "Calling sd_verlet_setup: 1"
409 c Rescale the stochastic forces and recalculate or restore
410 c the matrices of tinker integrator
411 if (itime_scal.gt.maxflag_stoch) then
412 if (large) write (iout,'(a,i5,a)')
413 & "Calculate matrices for stochastic step;",
414 & " itime_scal ",itime_scal
416 call sd_verlet_p_setup
418 call sd_verlet_ciccotti_setup
420 write (iout,'(2a,i3,a,i3,1h.)')
421 & "Warning: cannot store matrices for stochastic",
422 & " integration because the index",itime_scal,
423 & " is greater than",maxflag_stoch
424 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
425 & " integration Langevin algorithm for better efficiency."
426 else if (flag_stoch(itime_scal)) then
427 if (large) write (iout,'(a,i5,a,l1)')
428 & "Restore matrices for stochastic step; itime_scal ",
429 & itime_scal," flag ",flag_stoch(itime_scal)
432 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
433 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
434 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
435 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
436 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
437 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
441 if (large) write (iout,'(2a,i5,a,l1)')
442 & "Calculate & store matrices for stochastic step;",
443 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
445 call sd_verlet_p_setup
447 call sd_verlet_ciccotti_setup
449 flag_stoch(ifac_time)=.true.
452 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
453 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
454 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
455 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
456 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
457 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
461 fac_time=1.0d0/dsqrt(fac_time)
463 stochforcvec(i)=fac_time*stochforcvec(i)
466 else if (lang.eq.1) then
467 c Rescale the accelerations due to stochastic forces
468 fac_time=1.0d0/dsqrt(fac_time)
470 d_as_work(i)=d_as_work(i)*fac_time
473 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
474 & "itime",itime," Timestep scaled down to ",
475 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
477 c Second step of the velocity Verlet algorithm
482 else if (lang.eq.3) then
484 call sd_verlet2_ciccotti
486 else if (lang.eq.1) then
491 if (rattle) call rattle2
493 if (d_time.ne.d_time0) then
496 if (lang.eq.2 .or. lang.eq.3) then
497 if (large) write (iout,'(a)')
498 & "Restore original matrices for stochastic step"
499 c write (iout,*) "Calling sd_verlet_setup: 2"
500 c Restore the matrices of tinker integrator if the time step has been restored
503 pfric_mat(i,j)=pfric0_mat(i,j,0)
504 afric_mat(i,j)=afric0_mat(i,j,0)
505 vfric_mat(i,j)=vfric0_mat(i,j,0)
506 prand_mat(i,j)=prand0_mat(i,j,0)
507 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
508 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
517 c Calculate the kinetic and the total energy and the kinetic temperature
522 c write (iout,*) "step",itime," EK",EK," EK1",EK1
524 c Couple the system to Berendsen bath if needed
525 if (tbf .and. lang.eq.0) then
528 kinetic_T=2.0d0/(dimen3*Rb)*EK
529 c Backup the coordinates, velocities, and accelerations
533 d_t_old(j,i)=d_t(j,i)
534 d_a_old(j,i)=d_a(j,i)
538 if (mod(itime,ntwe).eq.0 .and. large) then
539 write (iout,*) "Velocities, step 2"
541 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
542 & (d_t(j,i+nres),j=1,3)
548 c-------------------------------------------------------------------------------
549 subroutine RESPA_step(itime)
550 c-------------------------------------------------------------------------------
551 c Perform a single RESPA step.
552 c-------------------------------------------------------------------------------
553 implicit real*8 (a-h,o-z)
557 integer IERROR,ERRCODE
559 include 'COMMON.SETUP'
560 include 'COMMON.CONTROL'
564 include 'COMMON.LANGEVIN'
566 include 'COMMON.LANGEVIN.lang0'
568 include 'COMMON.CHAIN'
569 include 'COMMON.DERIV'
571 include 'COMMON.LOCAL'
572 include 'COMMON.INTERACT'
573 include 'COMMON.IOUNITS'
574 include 'COMMON.NAMES'
575 include 'COMMON.TIME1'
576 double precision energia_short(0:n_ene),
577 & energia_long(0:n_ene)
578 double precision cm(3),L(3),vcm(3),incr(3)
579 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
580 & d_a_old0(3,0:maxres2)
581 integer ilen,count,rstcount
584 integer maxcount_scale /10/
586 double precision stochforcvec(MAXRES6)
587 common /stochcalc/ stochforcvec
590 common /cipiszcze/ itt
593 if (large.and. mod(itime,ntwe).eq.0) then
594 write (iout,*) "***************** RESPA itime",itime
595 write (iout,*) "Cartesian and internal coordinates: step 0"
597 call pdbout(0.0d0,"cipiszcze",iout)
599 write (iout,*) "Accelerations from long-range forces"
601 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
602 & (d_a(j,i+nres),j=1,3)
604 write (iout,*) "Velocities, step 0"
606 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
607 & (d_t(j,i+nres),j=1,3)
612 c Perform the initial RESPA step (increment velocities)
613 c write (iout,*) "*********************** RESPA ini"
616 if (mod(itime,ntwe).eq.0 .and. large) then
617 write (iout,*) "Velocities, end"
619 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
620 & (d_t(j,i+nres),j=1,3)
624 c Compute the short-range forces
631 call etotal_short(energia_short)
635 if (large.and. mod(itime,ntwe).eq.0) then
636 write (iout,*) "energia_short",energia_short(0)
637 write (iout,*) "Accelerations from short-range forces"
639 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
640 & (d_a(j,i+nres),j=1,3)
645 t_enegrad=t_enegrad+MPI_Wtime()-tt0
647 t_enegrad=t_enegrad+tcpu()-tt0
652 d_t_old(j,i)=d_t(j,i)
653 d_a_old(j,i)=d_a(j,i)
656 c 6/30/08 A-MTS: attempt at increasing the split number
659 dc_old0(j,i)=dc_old(j,i)
660 d_t_old0(j,i)=d_t_old(j,i)
661 d_a_old0(j,i)=d_a_old(j,i)
664 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
665 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
672 c write (iout,*) "itime",itime," ntime_split",ntime_split
673 c Split the time step
674 d_time=d_time0/ntime_split
675 c Perform the short-range RESPA steps (velocity Verlet increments of
676 c positions and velocities using short-range forces)
677 c write (iout,*) "*********************** RESPA split"
678 do itsplit=1,ntime_split
681 else if (lang.eq.2 .or. lang.eq.3) then
683 call stochastic_force(stochforcvec)
686 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
688 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
693 c First step of the velocity Verlet algorithm
698 else if (lang.eq.3) then
700 call sd_verlet1_ciccotti
702 else if (lang.eq.1) then
707 c Build the chain from the newly calculated coordinates
709 if (rattle) call rattle1
711 if (large.and. mod(itime,ntwe).eq.0) then
712 write (iout,*) "***** ITSPLIT",itsplit
713 write (iout,*) "Cartesian and internal coordinates: step 1"
714 call pdbout(0.0d0,"cipiszcze",iout)
717 write (iout,*) "Velocities, step 1"
719 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
720 & (d_t(j,i+nres),j=1,3)
729 c Calculate energy and forces
731 call etotal_short(energia_short)
733 c Get the new accelerations
736 if (large.and. mod(itime,ntwe).eq.0) then
737 write (iout,*)"energia_short",energia_short(0)
738 write (iout,*) "Accelerations from short-range forces"
740 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
741 & (d_a(j,i+nres),j=1,3)
746 c Determine maximum acceleration and scale down the timestep if needed
748 amax=amax/ntime_split**2
749 call predict_edrift(epdrift)
750 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
751 & write (iout,*) "amax",amax," damax",damax,
752 & " epdrift",epdrift," epdriftmax",epdriftmax
753 c Exit loop and try with increased split number if the change of
754 c acceleration is too big
755 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
756 if (ntime_split.lt.maxtime_split) then
758 ntime_split=ntime_split*2
761 dc_old(j,i)=dc_old0(j,i)
762 d_t_old(j,i)=d_t_old0(j,i)
763 d_a_old(j,i)=d_a_old0(j,i)
766 write (iout,*) "acceleration/energy drift too large",amax,
767 & epdrift," split increased to ",ntime_split," itime",itime,
772 & "Uh-hu. Bumpy landscape. Maximum splitting number",
774 & " already reached!!! Trying to carry on!"
778 t_enegrad=t_enegrad+MPI_Wtime()-tt0
780 t_enegrad=t_enegrad+tcpu()-tt0
782 c Second step of the velocity Verlet algorithm
787 else if (lang.eq.3) then
789 call sd_verlet2_ciccotti
791 else if (lang.eq.1) then
796 if (rattle) call rattle2
797 c Backup the coordinates, velocities, and accelerations
801 d_t_old(j,i)=d_t(j,i)
802 d_a_old(j,i)=d_a(j,i)
809 c Restore the time step
811 c Compute long-range forces
818 call etotal_long(energia_long)
822 t_enegrad=t_enegrad+MPI_Wtime()-tt0
824 t_enegrad=t_enegrad+tcpu()-tt0
826 c Compute accelerations from long-range forces
828 if (large.and. mod(itime,ntwe).eq.0) then
829 write (iout,*) "energia_long",energia_long(0)
830 write (iout,*) "Cartesian and internal coordinates: step 2"
832 call pdbout(0.0d0,"cipiszcze",iout)
834 write (iout,*) "Accelerations from long-range forces"
836 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
837 & (d_a(j,i+nres),j=1,3)
839 write (iout,*) "Velocities, step 2"
841 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
842 & (d_t(j,i+nres),j=1,3)
846 c Compute the final RESPA step (increment velocities)
847 c write (iout,*) "*********************** RESPA fin"
849 c Compute the complete potential energy
851 potEcomp(i)=energia_short(i)+energia_long(i)
853 potE=potEcomp(0)-potEcomp(20)
854 c potE=energia_short(0)+energia_long(0)
856 c Calculate the kinetic and the total energy and the kinetic temperature
859 c Couple the system to Berendsen bath if needed
860 if (tbf .and. lang.eq.0) then
863 kinetic_T=2.0d0/(dimen3*Rb)*EK
864 c Backup the coordinates, velocities, and accelerations
866 if (mod(itime,ntwe).eq.0 .and. large) then
867 write (iout,*) "Velocities, end"
869 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
870 & (d_t(j,i+nres),j=1,3)
876 c---------------------------------------------------------------------
878 c First and last RESPA step (incrementing velocities using long-range
880 implicit real*8 (a-h,o-z)
882 include 'COMMON.CONTROL'
885 include 'COMMON.CHAIN'
886 include 'COMMON.DERIV'
888 include 'COMMON.LOCAL'
889 include 'COMMON.INTERACT'
890 include 'COMMON.IOUNITS'
891 include 'COMMON.NAMES'
893 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
897 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
901 if (itype(i).ne.10) then
904 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
910 c-----------------------------------------------------------------
912 c Applying velocity Verlet algorithm - step 1 to coordinates
913 implicit real*8 (a-h,o-z)
915 include 'COMMON.CONTROL'
918 include 'COMMON.CHAIN'
919 include 'COMMON.DERIV'
921 include 'COMMON.LOCAL'
922 include 'COMMON.INTERACT'
923 include 'COMMON.IOUNITS'
924 include 'COMMON.NAMES'
925 double precision adt,adt2
928 write (iout,*) "VELVERLET1 START: DC"
930 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
931 & (dc(j,i+nres),j=1,3)
935 adt=d_a_old(j,0)*d_time
937 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
938 d_t_new(j,0)=d_t_old(j,0)+adt2
939 d_t(j,0)=d_t_old(j,0)+adt
943 adt=d_a_old(j,i)*d_time
945 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
946 d_t_new(j,i)=d_t_old(j,i)+adt2
947 d_t(j,i)=d_t_old(j,i)+adt
951 if (itype(i).ne.10) then
954 adt=d_a_old(j,inres)*d_time
956 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
957 d_t_new(j,inres)=d_t_old(j,inres)+adt2
958 d_t(j,inres)=d_t_old(j,inres)+adt
963 write (iout,*) "VELVERLET1 END: DC"
965 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
966 & (dc(j,i+nres),j=1,3)
971 c---------------------------------------------------------------------
973 c Step 2 of the velocity Verlet algorithm: update velocities
974 implicit real*8 (a-h,o-z)
976 include 'COMMON.CONTROL'
979 include 'COMMON.CHAIN'
980 include 'COMMON.DERIV'
982 include 'COMMON.LOCAL'
983 include 'COMMON.INTERACT'
984 include 'COMMON.IOUNITS'
985 include 'COMMON.NAMES'
987 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
991 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
995 if (itype(i).ne.10) then
998 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1004 c-----------------------------------------------------------------
1005 subroutine sddir_precalc
1006 c Applying velocity Verlet algorithm - step 1 to coordinates
1007 implicit real*8 (a-h,o-z)
1008 include 'DIMENSIONS'
1009 include 'COMMON.CONTROL'
1010 include 'COMMON.VAR'
1013 include 'COMMON.LANGEVIN'
1015 include 'COMMON.LANGEVIN.lang0'
1017 include 'COMMON.CHAIN'
1018 include 'COMMON.DERIV'
1019 include 'COMMON.GEO'
1020 include 'COMMON.LOCAL'
1021 include 'COMMON.INTERACT'
1022 include 'COMMON.IOUNITS'
1023 include 'COMMON.NAMES'
1024 double precision stochforcvec(MAXRES6)
1025 common /stochcalc/ stochforcvec
1027 c Compute friction and stochastic forces
1030 call stochastic_force(stochforcvec)
1032 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1033 c forces (d_as_work)
1035 call ginv_mult(fric_work, d_af_work)
1036 call ginv_mult(stochforcvec, d_as_work)
1039 c---------------------------------------------------------------------
1040 subroutine sddir_verlet1
1041 c Applying velocity Verlet algorithm - step 1 to velocities
1042 implicit real*8 (a-h,o-z)
1043 include 'DIMENSIONS'
1044 include 'COMMON.CONTROL'
1045 include 'COMMON.VAR'
1048 include 'COMMON.LANGEVIN'
1050 include 'COMMON.LANGEVIN.lang0'
1052 include 'COMMON.CHAIN'
1053 include 'COMMON.DERIV'
1054 include 'COMMON.GEO'
1055 include 'COMMON.LOCAL'
1056 include 'COMMON.INTERACT'
1057 include 'COMMON.IOUNITS'
1058 include 'COMMON.NAMES'
1059 c Revised 3/31/05 AL: correlation between random contributions to
1060 c position and velocity increments included.
1061 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1062 double precision adt,adt2
1064 c Add the contribution from BOTH friction and stochastic force to the
1065 c coordinates, but ONLY the contribution from the friction forces to velocities
1068 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1069 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1070 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1071 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1072 d_t(j,0)=d_t_old(j,0)+adt
1077 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1078 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1079 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1080 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1081 d_t(j,i)=d_t_old(j,i)+adt
1086 if (itype(i).ne.10) then
1089 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1090 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1091 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1092 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1093 d_t(j,inres)=d_t_old(j,inres)+adt
1100 c---------------------------------------------------------------------
1101 subroutine sddir_verlet2
1102 c Calculating the adjusted velocities for accelerations
1103 implicit real*8 (a-h,o-z)
1104 include 'DIMENSIONS'
1105 include 'COMMON.CONTROL'
1106 include 'COMMON.VAR'
1109 include 'COMMON.LANGEVIN'
1111 include 'COMMON.LANGEVIN.lang0'
1113 include 'COMMON.CHAIN'
1114 include 'COMMON.DERIV'
1115 include 'COMMON.GEO'
1116 include 'COMMON.LOCAL'
1117 include 'COMMON.INTERACT'
1118 include 'COMMON.IOUNITS'
1119 include 'COMMON.NAMES'
1120 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1121 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1122 c Revised 3/31/05 AL: correlation between random contributions to
1123 c position and velocity increments included.
1124 c The correlation coefficients are calculated at low-friction limit.
1125 c Also, friction forces are now not calculated with new velocities.
1127 c call friction_force
1128 call stochastic_force(stochforcvec)
1130 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1131 c forces (d_as_work)
1133 call ginv_mult(stochforcvec, d_as_work1)
1139 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1140 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1145 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1146 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1151 if (itype(i).ne.10) then
1154 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1155 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1156 & +cos60*d_as_work1(ind+j))*d_time
1163 c---------------------------------------------------------------------
1164 subroutine max_accel
1166 c Find the maximum difference in the accelerations of the the sites
1167 c at the beginning and the end of the time step.
1169 implicit real*8 (a-h,o-z)
1170 include 'DIMENSIONS'
1171 include 'COMMON.CONTROL'
1172 include 'COMMON.VAR'
1174 include 'COMMON.CHAIN'
1175 include 'COMMON.DERIV'
1176 include 'COMMON.GEO'
1177 include 'COMMON.LOCAL'
1178 include 'COMMON.INTERACT'
1179 include 'COMMON.IOUNITS'
1180 double precision aux(3),accel(3),accel_old(3),dacc
1182 c aux(j)=d_a(j,0)-d_a_old(j,0)
1183 accel_old(j)=d_a_old(j,0)
1190 c 7/3/08 changed to asymmetric difference
1192 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1193 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1194 accel(j)=accel(j)+0.5d0*d_a(j,i)
1195 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1196 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1197 dacc=dabs(accel(j)-accel_old(j))
1198 if (dacc.gt.amax) amax=dacc
1206 accel_old(j)=d_a_old(j,0)
1211 accel_old(j)=accel_old(j)+d_a_old(j,1)
1212 accel(j)=accel(j)+d_a(j,1)
1216 if (itype(i).ne.10) then
1218 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1219 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1220 accel(j)=accel(j)+d_a(j,i+nres)
1224 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1225 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1226 dacc=dabs(accel(j)-accel_old(j))
1227 if (dacc.gt.amax) amax=dacc
1231 accel_old(j)=accel_old(j)+d_a_old(j,i)
1232 accel(j)=accel(j)+d_a(j,i)
1233 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1238 c---------------------------------------------------------------------
1239 subroutine predict_edrift(epdrift)
1241 c Predict the drift of the potential energy
1243 implicit real*8 (a-h,o-z)
1244 include 'DIMENSIONS'
1245 include 'COMMON.CONTROL'
1246 include 'COMMON.VAR'
1248 include 'COMMON.CHAIN'
1249 include 'COMMON.DERIV'
1250 include 'COMMON.GEO'
1251 include 'COMMON.LOCAL'
1252 include 'COMMON.INTERACT'
1253 include 'COMMON.IOUNITS'
1254 include 'COMMON.MUCA'
1255 double precision epdrift,epdriftij
1256 c Drift of the potential energy
1262 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1263 if (lmuca) epdriftij=epdriftij*factor
1264 c write (iout,*) "back",i,j,epdriftij
1265 if (epdriftij.gt.epdrift) epdrift=epdriftij
1269 if (itype(i).ne.10) then
1272 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1273 if (lmuca) epdriftij=epdriftij*factor
1274 c write (iout,*) "side",i,j,epdriftij
1275 if (epdriftij.gt.epdrift) epdrift=epdriftij
1279 epdrift=0.5d0*epdrift*d_time*d_time
1280 c write (iout,*) "epdrift",epdrift
1283 c-----------------------------------------------------------------------
1284 subroutine verlet_bath
1286 c Coupling to the thermostat by using the Berendsen algorithm
1288 implicit real*8 (a-h,o-z)
1289 include 'DIMENSIONS'
1290 include 'COMMON.CONTROL'
1291 include 'COMMON.VAR'
1293 include 'COMMON.CHAIN'
1294 include 'COMMON.DERIV'
1295 include 'COMMON.GEO'
1296 include 'COMMON.LOCAL'
1297 include 'COMMON.INTERACT'
1298 include 'COMMON.IOUNITS'
1299 include 'COMMON.NAMES'
1300 double precision T_half,fact
1302 T_half=2.0d0/(dimen3*Rb)*EK
1303 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1304 c write(iout,*) "T_half", T_half
1305 c write(iout,*) "EK", EK
1306 c write(iout,*) "fact", fact
1308 d_t(j,0)=fact*d_t(j,0)
1312 d_t(j,i)=fact*d_t(j,i)
1316 if (itype(i).ne.10) then
1319 d_t(j,inres)=fact*d_t(j,inres)
1325 c---------------------------------------------------------
1327 c Set up the initial conditions of a MD simulation
1328 implicit real*8 (a-h,o-z)
1329 include 'DIMENSIONS'
1333 integer IERROR,ERRCODE
1335 include 'COMMON.SETUP'
1336 include 'COMMON.CONTROL'
1337 include 'COMMON.VAR'
1340 include 'COMMON.LANGEVIN'
1342 include 'COMMON.LANGEVIN.lang0'
1344 include 'COMMON.CHAIN'
1345 include 'COMMON.DERIV'
1346 include 'COMMON.GEO'
1347 include 'COMMON.LOCAL'
1348 include 'COMMON.INTERACT'
1349 include 'COMMON.IOUNITS'
1350 include 'COMMON.NAMES'
1351 include 'COMMON.REMD'
1352 real*8 energia_long(0:n_ene),
1353 & energia_short(0:n_ene),vcm(3),incr(3)
1354 double precision cm(3),L(3),xv,sigv,lowb,highb
1355 double precision varia(maxvar)
1363 c write(iout,*) "d_time", d_time
1364 c Compute the standard deviations of stochastic forces for Langevin dynamics
1365 c if the friction coefficients do not depend on surface area
1366 if (lang.gt.0 .and. .not.surfarea) then
1368 stdforcp(i)=stdfp*dsqrt(gamp)
1371 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
1374 c Open the pdb file for snapshotshots
1377 if (ilen(tmpdir).gt.0)
1378 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1379 & liczba(:ilen(liczba))//".pdb")
1381 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1385 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1386 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1387 & liczba(:ilen(liczba))//".x")
1388 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1391 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1392 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1393 & liczba(:ilen(liczba))//".cx")
1394 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1400 if (ilen(tmpdir).gt.0)
1401 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1402 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1404 if (ilen(tmpdir).gt.0)
1405 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1406 cartname=prefix(:ilen(prefix))//"_MD.cx"
1410 write (qstr,'(256(1h ))')
1413 iq = qinfrag(i,iset)*10
1414 iw = wfrag(i,iset)/100
1416 if(me.eq.king.or..not.out1file)
1417 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1418 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1423 iq = qinpair(i,iset)*10
1424 iw = wpair(i,iset)/100
1426 if(me.eq.king.or..not.out1file)
1427 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1428 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1432 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1434 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1436 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1438 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1442 if (restart1file) then
1444 & inquire(file=mremd_rst_name,exist=file_exist)
1445 write (*,*) me," Before broadcast: file_exist",file_exist
1446 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1448 write (*,*) me," After broadcast: file_exist",file_exist
1449 c inquire(file=mremd_rst_name,exist=file_exist)
1450 if(me.eq.king.or..not.out1file)
1451 & write(iout,*) "Initial state read by master and distributed"
1453 if (ilen(tmpdir).gt.0)
1454 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1455 & //liczba(:ilen(liczba))//'.rst')
1456 inquire(file=rest2name,exist=file_exist)
1459 if(.not.restart1file) then
1460 if(me.eq.king.or..not.out1file)
1461 & write(iout,*) "Initial state will be read from file ",
1462 & rest2name(:ilen(rest2name))
1465 call rescale_weights(t_bath)
1467 if(me.eq.king.or..not.out1file)then
1468 if (restart1file) then
1469 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1472 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1475 write(iout,*) "Initial velocities randomly generated"
1481 c Generate initial velocities
1482 if(me.eq.king.or..not.out1file)
1483 & write(iout,*) "Initial velocities randomly generated"
1487 c rest2name = prefix(:ilen(prefix))//'.rst'
1488 if(me.eq.king.or..not.out1file)then
1489 write(iout,*) "Initial backbone velocities"
1491 write(iout,*) (d_t(j,i),j=1,3)
1493 write(iout,*) "Initial side-chain velocities"
1495 write(iout,*) (d_t(j,i+nres),j=1,3)
1497 c Zeroing the total angular momentum of the system
1498 write(iout,*) "Calling the zero-angular
1499 & momentum subroutine"
1502 c Getting the potential energy and forces and velocities and accelerations
1504 c write (iout,*) "velocity of the center of the mass:"
1505 c write (iout,*) (vcm(j),j=1,3)
1507 d_t(j,0)=d_t(j,0)-vcm(j)
1509 c Removing the velocity of the center of mass
1511 if(me.eq.king.or..not.out1file)then
1512 write (iout,*) "vcm right after adjustment:"
1513 write (iout,*) (vcm(j),j=1,3)
1517 if(iranconf.ne.0) then
1519 print *, 'Calling OVERLAP_SC'
1520 call overlap_sc(fail)
1524 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1525 print *,'SC_move',nft_sc,etot
1526 if(me.eq.king.or..not.out1file)
1527 & write(iout,*) 'SC_move',nft_sc,etot
1531 print *, 'Calling MINIM_DC'
1532 call minim_dc(etot,iretcode,nfun)
1534 call geom_to_var(nvar,varia)
1535 print *,'Calling MINIMIZE.'
1536 call minimize(etot,varia,iretcode,nfun)
1537 call var_to_geom(nvar,varia)
1539 if(me.eq.king.or..not.out1file)
1540 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1543 call chainbuild_cart
1546 call verlet_bath(EK)
1548 kinetic_T=2.0d0/(dimen3*Rb)*EK
1549 if(me.eq.king.or..not.out1file)then
1554 call etotal(potEcomp)
1559 if (amax*d_time .gt. dvmax) then
1560 d_time=d_time*dvmax/amax
1561 if(me.eq.king.or..not.out1file) write (iout,*)
1562 & "Time step reduced to",d_time,
1563 & " because of too large initial acceleration."
1565 if(me.eq.king.or..not.out1file)then
1566 write(iout,*) "Potential energy and its components"
1567 write(iout,*) (potEcomp(i),i=0,n_ene)
1569 potE=potEcomp(0)-potEcomp(20)
1572 if (ntwe.ne.0) call statout(itime)
1573 if(me.eq.king.or..not.out1file)
1574 & write (iout,*) "Initial:",
1575 & " Kinetic energy",EK," potential energy",potE,
1576 & " total energy",totE," maximum acceleration ",
1579 write (iout,*) "Initial coordinates"
1581 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1582 & (c(j,i+nres),j=1,3)
1584 write (iout,*) "Initial dC"
1586 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1587 & (dc(j,i+nres),j=1,3)
1589 write (iout,*) "Initial velocities"
1591 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1592 & (d_t(j,i+nres),j=1,3)
1594 write (iout,*) "Initial accelerations"
1596 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1597 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1598 & (d_a(j,i+nres),j=1,3)
1604 d_t_old(j,i)=d_t(j,i)
1605 d_a_old(j,i)=d_a(j,i)
1607 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1616 call etotal_long(energia_long)
1617 if(.not.out1file .and. large) then
1618 write (iout,*) "energia_long",energia_long(0)
1619 write (iout,*) "Initial slow-force accelerations"
1621 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1622 & (d_a(j,i+nres),j=1,3)
1628 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1630 t_enegrad=t_enegrad+tcpu()-tt0
1632 c call etotal_short(energia_short)
1633 c write (iout,*) "energia_long",energia_long(0),
1634 c & " energia_short",energia_short(0),
1635 c & " total",energia_long(0)+energia_short(0)
1639 c-----------------------------------------------------------
1640 subroutine random_vel
1641 implicit real*8 (a-h,o-z)
1642 include 'DIMENSIONS'
1643 include 'COMMON.CONTROL'
1644 include 'COMMON.VAR'
1647 include 'COMMON.LANGEVIN'
1649 include 'COMMON.LANGEVIN.lang0'
1651 include 'COMMON.CHAIN'
1652 include 'COMMON.DERIV'
1653 include 'COMMON.GEO'
1654 include 'COMMON.LOCAL'
1655 include 'COMMON.INTERACT'
1656 include 'COMMON.IOUNITS'
1657 include 'COMMON.NAMES'
1658 include 'COMMON.TIME1'
1659 double precision xv,sigv,lowb,highb
1660 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1661 c First generate velocities in the eigenspace of the G matrix
1662 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1669 sigv=dsqrt((Rb*t_bath)/geigen(i))
1672 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1673 c write (iout,*) "ii",ii," d_t_work_new",d_t_work_new(ii)
1682 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1685 c write (iout,*) "Ek from eigenvectors",Ek1
1687 c Transform velocities to UNRES coordinate space
1693 d_t_work(ind)=d_t_work(ind)
1694 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1696 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1700 c Transfer to the d_t vector
1702 d_t(j,0)=d_t_work(j)
1708 d_t(j,i)=d_t_work(ind)
1712 if (itype(i).ne.10) then
1715 d_t(j,i+nres)=d_t_work(ind)
1720 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1721 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1726 c-----------------------------------------------------------
1727 subroutine sd_verlet_p_setup
1728 c Sets up the parameters of stochastic Verlet algorithm
1729 implicit real*8 (a-h,o-z)
1730 include 'DIMENSIONS'
1734 include 'COMMON.CONTROL'
1735 include 'COMMON.VAR'
1738 include 'COMMON.LANGEVIN'
1740 include 'COMMON.LANGEVIN.lang0'
1742 include 'COMMON.CHAIN'
1743 include 'COMMON.DERIV'
1744 include 'COMMON.GEO'
1745 include 'COMMON.LOCAL'
1746 include 'COMMON.INTERACT'
1747 include 'COMMON.IOUNITS'
1748 include 'COMMON.NAMES'
1749 include 'COMMON.TIME1'
1750 double precision emgdt(MAXRES6),
1751 & pterm,vterm,rho,rhoc,vsig,
1752 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1753 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1754 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1755 logical lprn /.false./
1756 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1757 double precision ktm
1764 c AL 8/17/04 Code adapted from tinker
1766 c Get the frictional and random terms for stochastic dynamics in the
1767 c eigenspace of mass-scaled UNRES friction matrix
1770 gdt = fricgam(i) * d_time
1772 c Stochastic dynamics reduces to simple MD for zero friction
1774 if (gdt .le. zero) then
1775 pfric_vec(i) = 1.0d0
1776 vfric_vec(i) = d_time
1777 afric_vec(i) = 0.5d0 * d_time * d_time
1778 prand_vec(i) = 0.0d0
1779 vrand_vec1(i) = 0.0d0
1780 vrand_vec2(i) = 0.0d0
1782 c Analytical expressions when friction coefficient is large
1785 if (gdt .ge. gdt_radius) then
1788 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1789 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1790 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1791 vterm = 1.0d0 - egdt**2
1792 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1794 c Use series expansions when friction coefficient is small
1805 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1806 & - gdt5/120.0d0 + gdt6/720.0d0
1807 & - gdt7/5040.0d0 + gdt8/40320.0d0
1808 & - gdt9/362880.0d0) / fricgam(i)**2
1809 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1810 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1811 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1812 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1813 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1814 & + 127.0d0*gdt9/90720.0d0
1815 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1816 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1817 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1818 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1819 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1820 & - 17.0d0*gdt2/1280.0d0
1821 & + 17.0d0*gdt3/6144.0d0
1822 & + 40967.0d0*gdt4/34406400.0d0
1823 & - 57203.0d0*gdt5/275251200.0d0
1824 & - 1429487.0d0*gdt6/13212057600.0d0)
1827 c Compute the scaling factors of random terms for the nonzero friction case
1829 ktm = 0.5d0*d_time/fricgam(i)
1830 psig = dsqrt(ktm*pterm) / fricgam(i)
1831 vsig = dsqrt(ktm*vterm)
1832 rhoc = dsqrt(1.0d0 - rho*rho)
1834 vrand_vec1(i) = vsig * rho
1835 vrand_vec2(i) = vsig * rhoc
1840 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1843 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1844 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1848 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1851 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1852 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1853 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1854 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1855 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1856 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1859 t_sdsetup=t_sdsetup+MPI_Wtime()
1861 t_sdsetup=t_sdsetup+tcpu()-tt0
1865 c-------------------------------------------------------------
1866 subroutine eigtransf1(n,ndim,ab,d,c)
1869 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
1875 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
1881 c-------------------------------------------------------------
1882 subroutine eigtransf(n,ndim,a,b,d,c)
1885 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
1891 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
1897 c-------------------------------------------------------------
1898 subroutine sd_verlet1
1899 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
1900 implicit real*8 (a-h,o-z)
1901 include 'DIMENSIONS'
1902 include 'COMMON.CONTROL'
1903 include 'COMMON.VAR'
1906 include 'COMMON.LANGEVIN'
1908 include 'COMMON.LANGEVIN.lang0'
1910 include 'COMMON.CHAIN'
1911 include 'COMMON.DERIV'
1912 include 'COMMON.GEO'
1913 include 'COMMON.LOCAL'
1914 include 'COMMON.INTERACT'
1915 include 'COMMON.IOUNITS'
1916 include 'COMMON.NAMES'
1917 double precision stochforcvec(MAXRES6)
1918 common /stochcalc/ stochforcvec
1919 logical lprn /.false./
1921 c write (iout,*) "dc_old"
1923 c write (iout,'(i5,3f10.5,5x,3f10.5)')
1924 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
1927 dc_work(j)=dc_old(j,0)
1928 d_t_work(j)=d_t_old(j,0)
1929 d_a_work(j)=d_a_old(j,0)
1934 dc_work(ind+j)=dc_old(j,i)
1935 d_t_work(ind+j)=d_t_old(j,i)
1936 d_a_work(ind+j)=d_a_old(j,i)
1941 if (itype(i).ne.10) then
1943 dc_work(ind+j)=dc_old(j,i+nres)
1944 d_t_work(ind+j)=d_t_old(j,i+nres)
1945 d_a_work(ind+j)=d_a_old(j,i+nres)
1953 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
1957 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
1958 & vfric_mat(i,j),afric_mat(i,j),
1959 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
1967 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
1968 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
1969 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
1970 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
1972 d_t_work_new(i)=ddt1+0.5d0*ddt2
1973 d_t_work(i)=ddt1+ddt2
1978 d_t(j,0)=d_t_work(j)
1983 dc(j,i)=dc_work(ind+j)
1984 d_t(j,i)=d_t_work(ind+j)
1989 if (itype(i).ne.10) then
1992 dc(j,inres)=dc_work(ind+j)
1993 d_t(j,inres)=d_t_work(ind+j)
2000 c--------------------------------------------------------------------------
2001 subroutine sd_verlet2
2002 c Calculating the adjusted velocities for accelerations
2003 implicit real*8 (a-h,o-z)
2004 include 'DIMENSIONS'
2005 include 'COMMON.CONTROL'
2006 include 'COMMON.VAR'
2009 include 'COMMON.LANGEVIN'
2011 include 'COMMON.LANGEVIN.lang0'
2013 include 'COMMON.CHAIN'
2014 include 'COMMON.DERIV'
2015 include 'COMMON.GEO'
2016 include 'COMMON.LOCAL'
2017 include 'COMMON.INTERACT'
2018 include 'COMMON.IOUNITS'
2019 include 'COMMON.NAMES'
2020 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2021 common /stochcalc/ stochforcvec
2023 c Compute the stochastic forces which contribute to velocity change
2025 call stochastic_force(stochforcvecV)
2032 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2033 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2034 & vrand_mat2(i,j)*stochforcvecV(j)
2036 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2040 d_t(j,0)=d_t_work(j)
2045 d_t(j,i)=d_t_work(ind+j)
2050 if (itype(i).ne.10) then
2053 d_t(j,inres)=d_t_work(ind+j)
2060 c-----------------------------------------------------------
2061 subroutine sd_verlet_ciccotti_setup
2062 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2064 implicit real*8 (a-h,o-z)
2065 include 'DIMENSIONS'
2069 include 'COMMON.CONTROL'
2070 include 'COMMON.VAR'
2073 include 'COMMON.LANGEVIN'
2075 include 'COMMON.LANGEVIN.lang0'
2077 include 'COMMON.CHAIN'
2078 include 'COMMON.DERIV'
2079 include 'COMMON.GEO'
2080 include 'COMMON.LOCAL'
2081 include 'COMMON.INTERACT'
2082 include 'COMMON.IOUNITS'
2083 include 'COMMON.NAMES'
2084 include 'COMMON.TIME1'
2085 double precision emgdt(MAXRES6),
2086 & pterm,vterm,rho,rhoc,vsig,
2087 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2088 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2089 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2090 logical lprn /.false./
2091 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2092 double precision ktm
2099 c AL 8/17/04 Code adapted from tinker
2101 c Get the frictional and random terms for stochastic dynamics in the
2102 c eigenspace of mass-scaled UNRES friction matrix
2105 write (iout,*) "i",i," fricgam",fricgam(i)
2106 gdt = fricgam(i) * d_time
2108 c Stochastic dynamics reduces to simple MD for zero friction
2110 if (gdt .le. zero) then
2111 pfric_vec(i) = 1.0d0
2112 vfric_vec(i) = d_time
2113 afric_vec(i) = 0.5d0*d_time*d_time
2114 prand_vec(i) = afric_vec(i)
2115 vrand_vec2(i) = vfric_vec(i)
2117 c Analytical expressions when friction coefficient is large
2122 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2123 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2124 prand_vec(i) = afric_vec(i)
2125 vrand_vec2(i) = vfric_vec(i)
2127 c Compute the scaling factors of random terms for the nonzero friction case
2129 c ktm = 0.5d0*d_time/fricgam(i)
2130 c psig = dsqrt(ktm*pterm) / fricgam(i)
2131 c vsig = dsqrt(ktm*vterm)
2132 c prand_vec(i) = psig*afric_vec(i)
2133 c vrand_vec2(i) = vsig*vfric_vec(i)
2138 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2141 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2142 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2146 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2148 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2149 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2150 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2151 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2152 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2154 t_sdsetup=t_sdsetup+MPI_Wtime()
2156 t_sdsetup=t_sdsetup+tcpu()-tt0
2160 c-------------------------------------------------------------
2161 subroutine sd_verlet1_ciccotti
2162 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2163 implicit real*8 (a-h,o-z)
2164 include 'DIMENSIONS'
2168 include 'COMMON.CONTROL'
2169 include 'COMMON.VAR'
2172 include 'COMMON.LANGEVIN'
2174 include 'COMMON.LANGEVIN.lang0'
2176 include 'COMMON.CHAIN'
2177 include 'COMMON.DERIV'
2178 include 'COMMON.GEO'
2179 include 'COMMON.LOCAL'
2180 include 'COMMON.INTERACT'
2181 include 'COMMON.IOUNITS'
2182 include 'COMMON.NAMES'
2183 double precision stochforcvec(MAXRES6)
2184 common /stochcalc/ stochforcvec
2185 logical lprn /.false./
2187 c write (iout,*) "dc_old"
2189 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2190 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2193 dc_work(j)=dc_old(j,0)
2194 d_t_work(j)=d_t_old(j,0)
2195 d_a_work(j)=d_a_old(j,0)
2200 dc_work(ind+j)=dc_old(j,i)
2201 d_t_work(ind+j)=d_t_old(j,i)
2202 d_a_work(ind+j)=d_a_old(j,i)
2207 if (itype(i).ne.10) then
2209 dc_work(ind+j)=dc_old(j,i+nres)
2210 d_t_work(ind+j)=d_t_old(j,i+nres)
2211 d_a_work(ind+j)=d_a_old(j,i+nres)
2220 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2224 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2225 & vfric_mat(i,j),afric_mat(i,j),
2226 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2234 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2235 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2236 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2237 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2239 d_t_work_new(i)=ddt1+0.5d0*ddt2
2240 d_t_work(i)=ddt1+ddt2
2245 d_t(j,0)=d_t_work(j)
2250 dc(j,i)=dc_work(ind+j)
2251 d_t(j,i)=d_t_work(ind+j)
2256 if (itype(i).ne.10) then
2259 dc(j,inres)=dc_work(ind+j)
2260 d_t(j,inres)=d_t_work(ind+j)
2267 c--------------------------------------------------------------------------
2268 subroutine sd_verlet2_ciccotti
2269 c Calculating the adjusted velocities for accelerations
2270 implicit real*8 (a-h,o-z)
2271 include 'DIMENSIONS'
2272 include 'COMMON.CONTROL'
2273 include 'COMMON.VAR'
2276 include 'COMMON.LANGEVIN'
2278 include 'COMMON.LANGEVIN.lang0'
2280 include 'COMMON.CHAIN'
2281 include 'COMMON.DERIV'
2282 include 'COMMON.GEO'
2283 include 'COMMON.LOCAL'
2284 include 'COMMON.INTERACT'
2285 include 'COMMON.IOUNITS'
2286 include 'COMMON.NAMES'
2287 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2288 common /stochcalc/ stochforcvec
2290 c Compute the stochastic forces which contribute to velocity change
2292 call stochastic_force(stochforcvecV)
2299 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2300 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2301 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2303 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2307 d_t(j,0)=d_t_work(j)
2312 d_t(j,i)=d_t_work(ind+j)
2317 if (itype(i).ne.10) then
2320 d_t(j,inres)=d_t_work(ind+j)