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
630 C 7/2/2009 commented out
632 c call etotal_short(energia_short)
635 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
638 d_a(j,i)=d_a_short(j,i)
642 if (large.and. mod(itime,ntwe).eq.0) then
643 write (iout,*) "energia_short",energia_short(0)
644 write (iout,*) "Accelerations from short-range forces"
646 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
647 & (d_a(j,i+nres),j=1,3)
652 t_enegrad=t_enegrad+MPI_Wtime()-tt0
654 t_enegrad=t_enegrad+tcpu()-tt0
659 d_t_old(j,i)=d_t(j,i)
660 d_a_old(j,i)=d_a(j,i)
663 c 6/30/08 A-MTS: attempt at increasing the split number
666 dc_old0(j,i)=dc_old(j,i)
667 d_t_old0(j,i)=d_t_old(j,i)
668 d_a_old0(j,i)=d_a_old(j,i)
671 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
672 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
679 c write (iout,*) "itime",itime," ntime_split",ntime_split
680 c Split the time step
681 d_time=d_time0/ntime_split
682 c Perform the short-range RESPA steps (velocity Verlet increments of
683 c positions and velocities using short-range forces)
684 c write (iout,*) "*********************** RESPA split"
685 do itsplit=1,ntime_split
688 else if (lang.eq.2 .or. lang.eq.3) then
690 call stochastic_force(stochforcvec)
693 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
695 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
700 c First step of the velocity Verlet algorithm
705 else if (lang.eq.3) then
707 call sd_verlet1_ciccotti
709 else if (lang.eq.1) then
714 c Build the chain from the newly calculated coordinates
716 if (rattle) call rattle1
718 if (large.and. mod(itime,ntwe).eq.0) then
719 write (iout,*) "***** ITSPLIT",itsplit
720 write (iout,*) "Cartesian and internal coordinates: step 1"
721 call pdbout(0.0d0,"cipiszcze",iout)
724 write (iout,*) "Velocities, step 1"
726 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
727 & (d_t(j,i+nres),j=1,3)
736 c Calculate energy and forces
738 call etotal_short(energia_short)
740 c Get the new accelerations
742 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
745 d_a_short(j,i)=d_a(j,i)
749 if (large.and. mod(itime,ntwe).eq.0) then
750 write (iout,*)"energia_short",energia_short(0)
751 write (iout,*) "Accelerations from short-range forces"
753 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
754 & (d_a(j,i+nres),j=1,3)
759 c Determine maximum acceleration and scale down the timestep if needed
761 amax=amax/ntime_split**2
762 call predict_edrift(epdrift)
763 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
764 & write (iout,*) "amax",amax," damax",damax,
765 & " epdrift",epdrift," epdriftmax",epdriftmax
766 c Exit loop and try with increased split number if the change of
767 c acceleration is too big
768 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
769 if (ntime_split.lt.maxtime_split) then
771 ntime_split=ntime_split*2
774 dc_old(j,i)=dc_old0(j,i)
775 d_t_old(j,i)=d_t_old0(j,i)
776 d_a_old(j,i)=d_a_old0(j,i)
779 write (iout,*) "acceleration/energy drift too large",amax,
780 & epdrift," split increased to ",ntime_split," itime",itime,
785 & "Uh-hu. Bumpy landscape. Maximum splitting number",
787 & " already reached!!! Trying to carry on!"
791 t_enegrad=t_enegrad+MPI_Wtime()-tt0
793 t_enegrad=t_enegrad+tcpu()-tt0
795 c Second step of the velocity Verlet algorithm
800 else if (lang.eq.3) then
802 call sd_verlet2_ciccotti
804 else if (lang.eq.1) then
809 if (rattle) call rattle2
810 c Backup the coordinates, velocities, and accelerations
814 d_t_old(j,i)=d_t(j,i)
815 d_a_old(j,i)=d_a(j,i)
822 c Restore the time step
824 c Compute long-range forces
831 call etotal_long(energia_long)
835 t_enegrad=t_enegrad+MPI_Wtime()-tt0
837 t_enegrad=t_enegrad+tcpu()-tt0
839 c Compute accelerations from long-range forces
841 if (large.and. mod(itime,ntwe).eq.0) then
842 write (iout,*) "energia_long",energia_long(0)
843 write (iout,*) "Cartesian and internal coordinates: step 2"
845 call pdbout(0.0d0,"cipiszcze",iout)
847 write (iout,*) "Accelerations from long-range forces"
849 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
850 & (d_a(j,i+nres),j=1,3)
852 write (iout,*) "Velocities, step 2"
854 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
855 & (d_t(j,i+nres),j=1,3)
859 c Compute the final RESPA step (increment velocities)
860 c write (iout,*) "*********************** RESPA fin"
862 c Compute the complete potential energy
864 potEcomp(i)=energia_short(i)+energia_long(i)
866 potE=potEcomp(0)-potEcomp(20)
867 c potE=energia_short(0)+energia_long(0)
869 c Calculate the kinetic and the total energy and the kinetic temperature
872 c Couple the system to Berendsen bath if needed
873 if (tbf .and. lang.eq.0) then
876 kinetic_T=2.0d0/(dimen3*Rb)*EK
877 c Backup the coordinates, velocities, and accelerations
879 if (mod(itime,ntwe).eq.0 .and. large) then
880 write (iout,*) "Velocities, end"
882 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
883 & (d_t(j,i+nres),j=1,3)
889 c---------------------------------------------------------------------
891 c First and last RESPA step (incrementing velocities using long-range
893 implicit real*8 (a-h,o-z)
895 include 'COMMON.CONTROL'
898 include 'COMMON.CHAIN'
899 include 'COMMON.DERIV'
901 include 'COMMON.LOCAL'
902 include 'COMMON.INTERACT'
903 include 'COMMON.IOUNITS'
904 include 'COMMON.NAMES'
906 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
910 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
914 if (itype(i).ne.10) then
917 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
923 c-----------------------------------------------------------------
925 c Applying velocity Verlet algorithm - step 1 to coordinates
926 implicit real*8 (a-h,o-z)
928 include 'COMMON.CONTROL'
931 include 'COMMON.CHAIN'
932 include 'COMMON.DERIV'
934 include 'COMMON.LOCAL'
935 include 'COMMON.INTERACT'
936 include 'COMMON.IOUNITS'
937 include 'COMMON.NAMES'
938 double precision adt,adt2
941 write (iout,*) "VELVERLET1 START: DC"
943 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
944 & (dc(j,i+nres),j=1,3)
948 adt=d_a_old(j,0)*d_time
950 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
951 d_t_new(j,0)=d_t_old(j,0)+adt2
952 d_t(j,0)=d_t_old(j,0)+adt
956 adt=d_a_old(j,i)*d_time
958 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
959 d_t_new(j,i)=d_t_old(j,i)+adt2
960 d_t(j,i)=d_t_old(j,i)+adt
964 if (itype(i).ne.10) then
967 adt=d_a_old(j,inres)*d_time
969 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
970 d_t_new(j,inres)=d_t_old(j,inres)+adt2
971 d_t(j,inres)=d_t_old(j,inres)+adt
976 write (iout,*) "VELVERLET1 END: DC"
978 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
979 & (dc(j,i+nres),j=1,3)
984 c---------------------------------------------------------------------
986 c Step 2 of the velocity Verlet algorithm: update velocities
987 implicit real*8 (a-h,o-z)
989 include 'COMMON.CONTROL'
992 include 'COMMON.CHAIN'
993 include 'COMMON.DERIV'
995 include 'COMMON.LOCAL'
996 include 'COMMON.INTERACT'
997 include 'COMMON.IOUNITS'
998 include 'COMMON.NAMES'
1000 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1004 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1008 if (itype(i).ne.10) then
1011 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1017 c-----------------------------------------------------------------
1018 subroutine sddir_precalc
1019 c Applying velocity Verlet algorithm - step 1 to coordinates
1020 implicit real*8 (a-h,o-z)
1021 include 'DIMENSIONS'
1022 include 'COMMON.CONTROL'
1023 include 'COMMON.VAR'
1026 include 'COMMON.LANGEVIN'
1028 include 'COMMON.LANGEVIN.lang0'
1030 include 'COMMON.CHAIN'
1031 include 'COMMON.DERIV'
1032 include 'COMMON.GEO'
1033 include 'COMMON.LOCAL'
1034 include 'COMMON.INTERACT'
1035 include 'COMMON.IOUNITS'
1036 include 'COMMON.NAMES'
1037 double precision stochforcvec(MAXRES6)
1038 common /stochcalc/ stochforcvec
1040 c Compute friction and stochastic forces
1043 call stochastic_force(stochforcvec)
1045 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1046 c forces (d_as_work)
1048 call ginv_mult(fric_work, d_af_work)
1049 call ginv_mult(stochforcvec, d_as_work)
1052 c---------------------------------------------------------------------
1053 subroutine sddir_verlet1
1054 c Applying velocity Verlet algorithm - step 1 to velocities
1055 implicit real*8 (a-h,o-z)
1056 include 'DIMENSIONS'
1057 include 'COMMON.CONTROL'
1058 include 'COMMON.VAR'
1061 include 'COMMON.LANGEVIN'
1063 include 'COMMON.LANGEVIN.lang0'
1065 include 'COMMON.CHAIN'
1066 include 'COMMON.DERIV'
1067 include 'COMMON.GEO'
1068 include 'COMMON.LOCAL'
1069 include 'COMMON.INTERACT'
1070 include 'COMMON.IOUNITS'
1071 include 'COMMON.NAMES'
1072 c Revised 3/31/05 AL: correlation between random contributions to
1073 c position and velocity increments included.
1074 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1075 double precision adt,adt2
1077 c Add the contribution from BOTH friction and stochastic force to the
1078 c coordinates, but ONLY the contribution from the friction forces to velocities
1081 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1082 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1083 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1084 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1085 d_t(j,0)=d_t_old(j,0)+adt
1090 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1091 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1092 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1093 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1094 d_t(j,i)=d_t_old(j,i)+adt
1099 if (itype(i).ne.10) then
1102 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1103 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1104 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1105 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1106 d_t(j,inres)=d_t_old(j,inres)+adt
1113 c---------------------------------------------------------------------
1114 subroutine sddir_verlet2
1115 c Calculating the adjusted velocities for accelerations
1116 implicit real*8 (a-h,o-z)
1117 include 'DIMENSIONS'
1118 include 'COMMON.CONTROL'
1119 include 'COMMON.VAR'
1122 include 'COMMON.LANGEVIN'
1124 include 'COMMON.LANGEVIN.lang0'
1126 include 'COMMON.CHAIN'
1127 include 'COMMON.DERIV'
1128 include 'COMMON.GEO'
1129 include 'COMMON.LOCAL'
1130 include 'COMMON.INTERACT'
1131 include 'COMMON.IOUNITS'
1132 include 'COMMON.NAMES'
1133 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1134 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1135 c Revised 3/31/05 AL: correlation between random contributions to
1136 c position and velocity increments included.
1137 c The correlation coefficients are calculated at low-friction limit.
1138 c Also, friction forces are now not calculated with new velocities.
1140 c call friction_force
1141 call stochastic_force(stochforcvec)
1143 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1144 c forces (d_as_work)
1146 call ginv_mult(stochforcvec, d_as_work1)
1152 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1153 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1158 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1159 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1164 if (itype(i).ne.10) then
1167 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1168 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1169 & +cos60*d_as_work1(ind+j))*d_time
1176 c---------------------------------------------------------------------
1177 subroutine max_accel
1179 c Find the maximum difference in the accelerations of the the sites
1180 c at the beginning and the end of the time step.
1182 implicit real*8 (a-h,o-z)
1183 include 'DIMENSIONS'
1184 include 'COMMON.CONTROL'
1185 include 'COMMON.VAR'
1187 include 'COMMON.CHAIN'
1188 include 'COMMON.DERIV'
1189 include 'COMMON.GEO'
1190 include 'COMMON.LOCAL'
1191 include 'COMMON.INTERACT'
1192 include 'COMMON.IOUNITS'
1193 double precision aux(3),accel(3),accel_old(3),dacc
1195 c aux(j)=d_a(j,0)-d_a_old(j,0)
1196 accel_old(j)=d_a_old(j,0)
1203 c 7/3/08 changed to asymmetric difference
1205 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1206 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1207 accel(j)=accel(j)+0.5d0*d_a(j,i)
1208 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1209 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1210 dacc=dabs(accel(j)-accel_old(j))
1211 if (dacc.gt.amax) amax=dacc
1219 accel_old(j)=d_a_old(j,0)
1224 accel_old(j)=accel_old(j)+d_a_old(j,1)
1225 accel(j)=accel(j)+d_a(j,1)
1229 if (itype(i).ne.10) then
1231 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1232 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1233 accel(j)=accel(j)+d_a(j,i+nres)
1237 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1238 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1239 dacc=dabs(accel(j)-accel_old(j))
1240 if (dacc.gt.amax) amax=dacc
1244 accel_old(j)=accel_old(j)+d_a_old(j,i)
1245 accel(j)=accel(j)+d_a(j,i)
1246 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1251 c---------------------------------------------------------------------
1252 subroutine predict_edrift(epdrift)
1254 c Predict the drift of the potential energy
1256 implicit real*8 (a-h,o-z)
1257 include 'DIMENSIONS'
1258 include 'COMMON.CONTROL'
1259 include 'COMMON.VAR'
1261 include 'COMMON.CHAIN'
1262 include 'COMMON.DERIV'
1263 include 'COMMON.GEO'
1264 include 'COMMON.LOCAL'
1265 include 'COMMON.INTERACT'
1266 include 'COMMON.IOUNITS'
1267 include 'COMMON.MUCA'
1268 double precision epdrift,epdriftij
1269 c Drift of the potential energy
1275 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1276 if (lmuca) epdriftij=epdriftij*factor
1277 c write (iout,*) "back",i,j,epdriftij
1278 if (epdriftij.gt.epdrift) epdrift=epdriftij
1282 if (itype(i).ne.10) then
1285 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1286 if (lmuca) epdriftij=epdriftij*factor
1287 c write (iout,*) "side",i,j,epdriftij
1288 if (epdriftij.gt.epdrift) epdrift=epdriftij
1292 epdrift=0.5d0*epdrift*d_time*d_time
1293 c write (iout,*) "epdrift",epdrift
1296 c-----------------------------------------------------------------------
1297 subroutine verlet_bath
1299 c Coupling to the thermostat by using the Berendsen algorithm
1301 implicit real*8 (a-h,o-z)
1302 include 'DIMENSIONS'
1303 include 'COMMON.CONTROL'
1304 include 'COMMON.VAR'
1306 include 'COMMON.CHAIN'
1307 include 'COMMON.DERIV'
1308 include 'COMMON.GEO'
1309 include 'COMMON.LOCAL'
1310 include 'COMMON.INTERACT'
1311 include 'COMMON.IOUNITS'
1312 include 'COMMON.NAMES'
1313 double precision T_half,fact
1315 T_half=2.0d0/(dimen3*Rb)*EK
1316 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1317 c write(iout,*) "T_half", T_half
1318 c write(iout,*) "EK", EK
1319 c write(iout,*) "fact", fact
1321 d_t(j,0)=fact*d_t(j,0)
1325 d_t(j,i)=fact*d_t(j,i)
1329 if (itype(i).ne.10) then
1332 d_t(j,inres)=fact*d_t(j,inres)
1338 c---------------------------------------------------------
1340 c Set up the initial conditions of a MD simulation
1341 implicit real*8 (a-h,o-z)
1342 include 'DIMENSIONS'
1346 integer IERROR,ERRCODE
1348 include 'COMMON.SETUP'
1349 include 'COMMON.CONTROL'
1350 include 'COMMON.VAR'
1353 include 'COMMON.LANGEVIN'
1355 include 'COMMON.LANGEVIN.lang0'
1357 include 'COMMON.CHAIN'
1358 include 'COMMON.DERIV'
1359 include 'COMMON.GEO'
1360 include 'COMMON.LOCAL'
1361 include 'COMMON.INTERACT'
1362 include 'COMMON.IOUNITS'
1363 include 'COMMON.NAMES'
1364 include 'COMMON.REMD'
1365 real*8 energia_long(0:n_ene),
1366 & energia_short(0:n_ene),vcm(3),incr(3)
1367 double precision cm(3),L(3),xv,sigv,lowb,highb
1368 double precision varia(maxvar)
1376 c write(iout,*) "d_time", d_time
1377 c Compute the standard deviations of stochastic forces for Langevin dynamics
1378 c if the friction coefficients do not depend on surface area
1379 if (lang.gt.0 .and. .not.surfarea) then
1381 stdforcp(i)=stdfp*dsqrt(gamp)
1384 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
1387 c Open the pdb file for snapshotshots
1390 if (ilen(tmpdir).gt.0)
1391 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1392 & liczba(:ilen(liczba))//".pdb")
1394 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1398 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1399 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1400 & liczba(:ilen(liczba))//".x")
1401 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1404 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1405 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1406 & liczba(:ilen(liczba))//".cx")
1407 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1413 if (ilen(tmpdir).gt.0)
1414 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1415 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1417 if (ilen(tmpdir).gt.0)
1418 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1419 cartname=prefix(:ilen(prefix))//"_MD.cx"
1423 write (qstr,'(256(1h ))')
1426 iq = qinfrag(i,iset)*10
1427 iw = wfrag(i,iset)/100
1429 if(me.eq.king.or..not.out1file)
1430 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1431 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1436 iq = qinpair(i,iset)*10
1437 iw = wpair(i,iset)/100
1439 if(me.eq.king.or..not.out1file)
1440 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1441 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1445 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1447 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1449 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1451 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1455 if (restart1file) then
1457 & inquire(file=mremd_rst_name,exist=file_exist)
1458 write (*,*) me," Before broadcast: file_exist",file_exist
1459 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1461 write (*,*) me," After broadcast: file_exist",file_exist
1462 c inquire(file=mremd_rst_name,exist=file_exist)
1463 if(me.eq.king.or..not.out1file)
1464 & write(iout,*) "Initial state read by master and distributed"
1466 if (ilen(tmpdir).gt.0)
1467 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1468 & //liczba(:ilen(liczba))//'.rst')
1469 inquire(file=rest2name,exist=file_exist)
1472 if(.not.restart1file) then
1473 if(me.eq.king.or..not.out1file)
1474 & write(iout,*) "Initial state will be read from file ",
1475 & rest2name(:ilen(rest2name))
1478 call rescale_weights(t_bath)
1480 if(me.eq.king.or..not.out1file)then
1481 if (restart1file) then
1482 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1485 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1488 write(iout,*) "Initial velocities randomly generated"
1494 c Generate initial velocities
1495 if(me.eq.king.or..not.out1file)
1496 & write(iout,*) "Initial velocities randomly generated"
1500 c rest2name = prefix(:ilen(prefix))//'.rst'
1501 if(me.eq.king.or..not.out1file)then
1502 write(iout,*) "Initial backbone velocities"
1504 write(iout,*) (d_t(j,i),j=1,3)
1506 write(iout,*) "Initial side-chain velocities"
1508 write(iout,*) (d_t(j,i+nres),j=1,3)
1510 c Zeroing the total angular momentum of the system
1511 write(iout,*) "Calling the zero-angular
1512 & momentum subroutine"
1515 c Getting the potential energy and forces and velocities and accelerations
1517 c write (iout,*) "velocity of the center of the mass:"
1518 c write (iout,*) (vcm(j),j=1,3)
1520 d_t(j,0)=d_t(j,0)-vcm(j)
1522 c Removing the velocity of the center of mass
1524 if(me.eq.king.or..not.out1file)then
1525 write (iout,*) "vcm right after adjustment:"
1526 write (iout,*) (vcm(j),j=1,3)
1530 if(iranconf.ne.0) then
1532 print *, 'Calling OVERLAP_SC'
1533 call overlap_sc(fail)
1537 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1538 print *,'SC_move',nft_sc,etot
1539 if(me.eq.king.or..not.out1file)
1540 & write(iout,*) 'SC_move',nft_sc,etot
1544 print *, 'Calling MINIM_DC'
1545 call minim_dc(etot,iretcode,nfun)
1547 call geom_to_var(nvar,varia)
1548 print *,'Calling MINIMIZE.'
1549 call minimize(etot,varia,iretcode,nfun)
1550 call var_to_geom(nvar,varia)
1552 if(me.eq.king.or..not.out1file)
1553 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1556 call chainbuild_cart
1559 call verlet_bath(EK)
1561 kinetic_T=2.0d0/(dimen3*Rb)*EK
1562 if(me.eq.king.or..not.out1file)then
1567 call etotal(potEcomp)
1572 if (amax*d_time .gt. dvmax) then
1573 d_time=d_time*dvmax/amax
1574 if(me.eq.king.or..not.out1file) write (iout,*)
1575 & "Time step reduced to",d_time,
1576 & " because of too large initial acceleration."
1578 if(me.eq.king.or..not.out1file)then
1579 write(iout,*) "Potential energy and its components"
1580 write(iout,*) (potEcomp(i),i=0,n_ene)
1582 potE=potEcomp(0)-potEcomp(20)
1585 if (ntwe.ne.0) call statout(itime)
1586 if(me.eq.king.or..not.out1file)
1587 & write (iout,*) "Initial:",
1588 & " Kinetic energy",EK," potential energy",potE,
1589 & " total energy",totE," maximum acceleration ",
1592 write (iout,*) "Initial coordinates"
1594 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1595 & (c(j,i+nres),j=1,3)
1597 write (iout,*) "Initial dC"
1599 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1600 & (dc(j,i+nres),j=1,3)
1602 write (iout,*) "Initial velocities"
1604 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1605 & (d_t(j,i+nres),j=1,3)
1607 write (iout,*) "Initial accelerations"
1609 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1610 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1611 & (d_a(j,i+nres),j=1,3)
1617 d_t_old(j,i)=d_t(j,i)
1618 d_a_old(j,i)=d_a(j,i)
1620 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1629 call etotal_short(energia_short)
1632 if(.not.out1file .and. large) then
1633 write (iout,*) "energia_long",energia_long(0),
1634 & " energia_short",energia_short(0),
1635 & " total",energia_long(0)+energia_short(0)
1636 write (iout,*) "Initial fast-force accelerations"
1638 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1639 & (d_a(j,i+nres),j=1,3)
1642 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1645 d_a_short(j,i)=d_a(j,i)
1649 call etotal_long(energia_long)
1652 if(.not.out1file .and. large) then
1653 write (iout,*) "energia_long",energia_long(0)
1654 write (iout,*) "Initial slow-force accelerations"
1656 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1657 & (d_a(j,i+nres),j=1,3)
1661 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1663 t_enegrad=t_enegrad+tcpu()-tt0
1668 c-----------------------------------------------------------
1669 subroutine random_vel
1670 implicit real*8 (a-h,o-z)
1671 include 'DIMENSIONS'
1672 include 'COMMON.CONTROL'
1673 include 'COMMON.VAR'
1676 include 'COMMON.LANGEVIN'
1678 include 'COMMON.LANGEVIN.lang0'
1680 include 'COMMON.CHAIN'
1681 include 'COMMON.DERIV'
1682 include 'COMMON.GEO'
1683 include 'COMMON.LOCAL'
1684 include 'COMMON.INTERACT'
1685 include 'COMMON.IOUNITS'
1686 include 'COMMON.NAMES'
1687 include 'COMMON.TIME1'
1688 double precision xv,sigv,lowb,highb
1689 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1690 c First generate velocities in the eigenspace of the G matrix
1691 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1698 sigv=dsqrt((Rb*t_bath)/geigen(i))
1701 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1702 c write (iout,*) "ii",ii," d_t_work_new",d_t_work_new(ii)
1711 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1714 c write (iout,*) "Ek from eigenvectors",Ek1
1716 c Transform velocities to UNRES coordinate space
1722 d_t_work(ind)=d_t_work(ind)
1723 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1725 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1729 c Transfer to the d_t vector
1731 d_t(j,0)=d_t_work(j)
1737 d_t(j,i)=d_t_work(ind)
1741 if (itype(i).ne.10) then
1744 d_t(j,i+nres)=d_t_work(ind)
1749 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1750 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1755 c-----------------------------------------------------------
1756 subroutine sd_verlet_p_setup
1757 c Sets up the parameters of stochastic Verlet algorithm
1758 implicit real*8 (a-h,o-z)
1759 include 'DIMENSIONS'
1763 include 'COMMON.CONTROL'
1764 include 'COMMON.VAR'
1767 include 'COMMON.LANGEVIN'
1769 include 'COMMON.LANGEVIN.lang0'
1771 include 'COMMON.CHAIN'
1772 include 'COMMON.DERIV'
1773 include 'COMMON.GEO'
1774 include 'COMMON.LOCAL'
1775 include 'COMMON.INTERACT'
1776 include 'COMMON.IOUNITS'
1777 include 'COMMON.NAMES'
1778 include 'COMMON.TIME1'
1779 double precision emgdt(MAXRES6),
1780 & pterm,vterm,rho,rhoc,vsig,
1781 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1782 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1783 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1784 logical lprn /.false./
1785 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1786 double precision ktm
1793 c AL 8/17/04 Code adapted from tinker
1795 c Get the frictional and random terms for stochastic dynamics in the
1796 c eigenspace of mass-scaled UNRES friction matrix
1799 gdt = fricgam(i) * d_time
1801 c Stochastic dynamics reduces to simple MD for zero friction
1803 if (gdt .le. zero) then
1804 pfric_vec(i) = 1.0d0
1805 vfric_vec(i) = d_time
1806 afric_vec(i) = 0.5d0 * d_time * d_time
1807 prand_vec(i) = 0.0d0
1808 vrand_vec1(i) = 0.0d0
1809 vrand_vec2(i) = 0.0d0
1811 c Analytical expressions when friction coefficient is large
1814 if (gdt .ge. gdt_radius) then
1817 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1818 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1819 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1820 vterm = 1.0d0 - egdt**2
1821 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1823 c Use series expansions when friction coefficient is small
1834 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1835 & - gdt5/120.0d0 + gdt6/720.0d0
1836 & - gdt7/5040.0d0 + gdt8/40320.0d0
1837 & - gdt9/362880.0d0) / fricgam(i)**2
1838 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1839 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1840 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1841 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1842 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1843 & + 127.0d0*gdt9/90720.0d0
1844 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1845 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1846 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1847 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1848 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1849 & - 17.0d0*gdt2/1280.0d0
1850 & + 17.0d0*gdt3/6144.0d0
1851 & + 40967.0d0*gdt4/34406400.0d0
1852 & - 57203.0d0*gdt5/275251200.0d0
1853 & - 1429487.0d0*gdt6/13212057600.0d0)
1856 c Compute the scaling factors of random terms for the nonzero friction case
1858 ktm = 0.5d0*d_time/fricgam(i)
1859 psig = dsqrt(ktm*pterm) / fricgam(i)
1860 vsig = dsqrt(ktm*vterm)
1861 rhoc = dsqrt(1.0d0 - rho*rho)
1863 vrand_vec1(i) = vsig * rho
1864 vrand_vec2(i) = vsig * rhoc
1869 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1872 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1873 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1877 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1880 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1881 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1882 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1883 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1884 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1885 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1888 t_sdsetup=t_sdsetup+MPI_Wtime()
1890 t_sdsetup=t_sdsetup+tcpu()-tt0
1894 c-------------------------------------------------------------
1895 subroutine eigtransf1(n,ndim,ab,d,c)
1898 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
1904 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
1910 c-------------------------------------------------------------
1911 subroutine eigtransf(n,ndim,a,b,d,c)
1914 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
1920 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
1926 c-------------------------------------------------------------
1927 subroutine sd_verlet1
1928 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
1929 implicit real*8 (a-h,o-z)
1930 include 'DIMENSIONS'
1931 include 'COMMON.CONTROL'
1932 include 'COMMON.VAR'
1935 include 'COMMON.LANGEVIN'
1937 include 'COMMON.LANGEVIN.lang0'
1939 include 'COMMON.CHAIN'
1940 include 'COMMON.DERIV'
1941 include 'COMMON.GEO'
1942 include 'COMMON.LOCAL'
1943 include 'COMMON.INTERACT'
1944 include 'COMMON.IOUNITS'
1945 include 'COMMON.NAMES'
1946 double precision stochforcvec(MAXRES6)
1947 common /stochcalc/ stochforcvec
1948 logical lprn /.false./
1950 c write (iout,*) "dc_old"
1952 c write (iout,'(i5,3f10.5,5x,3f10.5)')
1953 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
1956 dc_work(j)=dc_old(j,0)
1957 d_t_work(j)=d_t_old(j,0)
1958 d_a_work(j)=d_a_old(j,0)
1963 dc_work(ind+j)=dc_old(j,i)
1964 d_t_work(ind+j)=d_t_old(j,i)
1965 d_a_work(ind+j)=d_a_old(j,i)
1970 if (itype(i).ne.10) then
1972 dc_work(ind+j)=dc_old(j,i+nres)
1973 d_t_work(ind+j)=d_t_old(j,i+nres)
1974 d_a_work(ind+j)=d_a_old(j,i+nres)
1982 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
1986 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
1987 & vfric_mat(i,j),afric_mat(i,j),
1988 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
1996 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
1997 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
1998 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
1999 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2001 d_t_work_new(i)=ddt1+0.5d0*ddt2
2002 d_t_work(i)=ddt1+ddt2
2007 d_t(j,0)=d_t_work(j)
2012 dc(j,i)=dc_work(ind+j)
2013 d_t(j,i)=d_t_work(ind+j)
2018 if (itype(i).ne.10) then
2021 dc(j,inres)=dc_work(ind+j)
2022 d_t(j,inres)=d_t_work(ind+j)
2029 c--------------------------------------------------------------------------
2030 subroutine sd_verlet2
2031 c Calculating the adjusted velocities for accelerations
2032 implicit real*8 (a-h,o-z)
2033 include 'DIMENSIONS'
2034 include 'COMMON.CONTROL'
2035 include 'COMMON.VAR'
2038 include 'COMMON.LANGEVIN'
2040 include 'COMMON.LANGEVIN.lang0'
2042 include 'COMMON.CHAIN'
2043 include 'COMMON.DERIV'
2044 include 'COMMON.GEO'
2045 include 'COMMON.LOCAL'
2046 include 'COMMON.INTERACT'
2047 include 'COMMON.IOUNITS'
2048 include 'COMMON.NAMES'
2049 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2050 common /stochcalc/ stochforcvec
2052 c Compute the stochastic forces which contribute to velocity change
2054 call stochastic_force(stochforcvecV)
2061 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2062 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2063 & vrand_mat2(i,j)*stochforcvecV(j)
2065 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2069 d_t(j,0)=d_t_work(j)
2074 d_t(j,i)=d_t_work(ind+j)
2079 if (itype(i).ne.10) then
2082 d_t(j,inres)=d_t_work(ind+j)
2089 c-----------------------------------------------------------
2090 subroutine sd_verlet_ciccotti_setup
2091 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2093 implicit real*8 (a-h,o-z)
2094 include 'DIMENSIONS'
2098 include 'COMMON.CONTROL'
2099 include 'COMMON.VAR'
2102 include 'COMMON.LANGEVIN'
2104 include 'COMMON.LANGEVIN.lang0'
2106 include 'COMMON.CHAIN'
2107 include 'COMMON.DERIV'
2108 include 'COMMON.GEO'
2109 include 'COMMON.LOCAL'
2110 include 'COMMON.INTERACT'
2111 include 'COMMON.IOUNITS'
2112 include 'COMMON.NAMES'
2113 include 'COMMON.TIME1'
2114 double precision emgdt(MAXRES6),
2115 & pterm,vterm,rho,rhoc,vsig,
2116 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2117 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2118 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2119 logical lprn /.false./
2120 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2121 double precision ktm
2128 c AL 8/17/04 Code adapted from tinker
2130 c Get the frictional and random terms for stochastic dynamics in the
2131 c eigenspace of mass-scaled UNRES friction matrix
2134 write (iout,*) "i",i," fricgam",fricgam(i)
2135 gdt = fricgam(i) * d_time
2137 c Stochastic dynamics reduces to simple MD for zero friction
2139 if (gdt .le. zero) then
2140 pfric_vec(i) = 1.0d0
2141 vfric_vec(i) = d_time
2142 afric_vec(i) = 0.5d0*d_time*d_time
2143 prand_vec(i) = afric_vec(i)
2144 vrand_vec2(i) = vfric_vec(i)
2146 c Analytical expressions when friction coefficient is large
2151 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2152 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2153 prand_vec(i) = afric_vec(i)
2154 vrand_vec2(i) = vfric_vec(i)
2156 c Compute the scaling factors of random terms for the nonzero friction case
2158 c ktm = 0.5d0*d_time/fricgam(i)
2159 c psig = dsqrt(ktm*pterm) / fricgam(i)
2160 c vsig = dsqrt(ktm*vterm)
2161 c prand_vec(i) = psig*afric_vec(i)
2162 c vrand_vec2(i) = vsig*vfric_vec(i)
2167 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2170 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2171 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2175 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2177 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2178 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2179 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2180 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2181 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2183 t_sdsetup=t_sdsetup+MPI_Wtime()
2185 t_sdsetup=t_sdsetup+tcpu()-tt0
2189 c-------------------------------------------------------------
2190 subroutine sd_verlet1_ciccotti
2191 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2192 implicit real*8 (a-h,o-z)
2193 include 'DIMENSIONS'
2197 include 'COMMON.CONTROL'
2198 include 'COMMON.VAR'
2201 include 'COMMON.LANGEVIN'
2203 include 'COMMON.LANGEVIN.lang0'
2205 include 'COMMON.CHAIN'
2206 include 'COMMON.DERIV'
2207 include 'COMMON.GEO'
2208 include 'COMMON.LOCAL'
2209 include 'COMMON.INTERACT'
2210 include 'COMMON.IOUNITS'
2211 include 'COMMON.NAMES'
2212 double precision stochforcvec(MAXRES6)
2213 common /stochcalc/ stochforcvec
2214 logical lprn /.false./
2216 c write (iout,*) "dc_old"
2218 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2219 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2222 dc_work(j)=dc_old(j,0)
2223 d_t_work(j)=d_t_old(j,0)
2224 d_a_work(j)=d_a_old(j,0)
2229 dc_work(ind+j)=dc_old(j,i)
2230 d_t_work(ind+j)=d_t_old(j,i)
2231 d_a_work(ind+j)=d_a_old(j,i)
2236 if (itype(i).ne.10) then
2238 dc_work(ind+j)=dc_old(j,i+nres)
2239 d_t_work(ind+j)=d_t_old(j,i+nres)
2240 d_a_work(ind+j)=d_a_old(j,i+nres)
2249 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2253 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2254 & vfric_mat(i,j),afric_mat(i,j),
2255 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2263 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2264 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2265 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2266 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2268 d_t_work_new(i)=ddt1+0.5d0*ddt2
2269 d_t_work(i)=ddt1+ddt2
2274 d_t(j,0)=d_t_work(j)
2279 dc(j,i)=dc_work(ind+j)
2280 d_t(j,i)=d_t_work(ind+j)
2285 if (itype(i).ne.10) then
2288 dc(j,inres)=dc_work(ind+j)
2289 d_t(j,inres)=d_t_work(ind+j)
2296 c--------------------------------------------------------------------------
2297 subroutine sd_verlet2_ciccotti
2298 c Calculating the adjusted velocities for accelerations
2299 implicit real*8 (a-h,o-z)
2300 include 'DIMENSIONS'
2301 include 'COMMON.CONTROL'
2302 include 'COMMON.VAR'
2305 include 'COMMON.LANGEVIN'
2307 include 'COMMON.LANGEVIN.lang0'
2309 include 'COMMON.CHAIN'
2310 include 'COMMON.DERIV'
2311 include 'COMMON.GEO'
2312 include 'COMMON.LOCAL'
2313 include 'COMMON.INTERACT'
2314 include 'COMMON.IOUNITS'
2315 include 'COMMON.NAMES'
2316 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2317 common /stochcalc/ stochforcvec
2319 c Compute the stochastic forces which contribute to velocity change
2321 call stochastic_force(stochforcvecV)
2328 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2329 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2330 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2332 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2336 d_t(j,0)=d_t_work(j)
2341 d_t(j,i)=d_t_work(ind+j)
2346 if (itype(i).ne.10) then
2349 d_t(j,inres)=d_t_work(ind+j)