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 '
263 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
265 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
266 & " time_fricmatmult",time_fricmatmult," time_fsample ",
271 c-------------------------------------------------------------------------------
272 subroutine velverlet_step(itime)
273 c-------------------------------------------------------------------------------
274 c Perform a single velocity Verlet step; the time step can be rescaled if
275 c increments in accelerations exceed the threshold
276 c-------------------------------------------------------------------------------
277 implicit real*8 (a-h,o-z)
281 integer ierror,ierrcode
283 include 'COMMON.SETUP'
284 include 'COMMON.CONTROL'
288 include 'COMMON.LANGEVIN'
290 include 'COMMON.LANGEVIN.lang0'
292 include 'COMMON.CHAIN'
293 include 'COMMON.DERIV'
295 include 'COMMON.LOCAL'
296 include 'COMMON.INTERACT'
297 include 'COMMON.IOUNITS'
298 include 'COMMON.NAMES'
299 include 'COMMON.TIME1'
300 include 'COMMON.MUCA'
301 double precision vcm(3),incr(3)
302 double precision cm(3),L(3)
303 integer ilen,count,rstcount
306 integer maxcount_scale /20/
308 double precision stochforcvec(MAXRES6)
309 common /stochcalc/ stochforcvec
312 double precision HNose1,HNose,HNose_nh,H,vtnp(maxres6)
313 double precision vtnp_(maxres6),vtnp_a(maxres6)
319 else if (lang.eq.2 .or. lang.eq.3) then
321 call stochastic_force(stochforcvec)
324 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
326 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
333 icount_scale=icount_scale+1
334 if (icount_scale.gt.maxcount_scale) then
336 & "ERROR: too many attempts at scaling down the time step. ",
337 & "amax=",amax,"epdrift=",epdrift,
338 & "damax=",damax,"edriftmax=",edriftmax,
342 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
346 c First step of the velocity Verlet algorithm
351 else if (lang.eq.3) then
353 call sd_verlet1_ciccotti
355 else if (lang.eq.1) then
363 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
366 d_t_old(j,i)=d_t_old(j,i)*scale_nh
372 c Build the chain from the newly calculated coordinates
374 if (rattle) call rattle1
376 if (large.and. mod(itime,ntwe).eq.0) then
377 write (iout,*) "Cartesian and internal coordinates: step 1"
382 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
383 & (dc(j,i+nres),j=1,3)
385 write (iout,*) "Accelerations"
387 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
388 & (d_a(j,i+nres),j=1,3)
390 write (iout,*) "Velocities, step 1"
392 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
393 & (d_t(j,i+nres),j=1,3)
402 c Calculate energy and forces
404 call etotal(potEcomp)
407 t_etotal=t_etotal+MPI_Wtime()-tt0
409 t_etotal=t_etotal+tcpu()-tt0
413 potE=potEcomp(0)-potEcomp(20)
415 c Get the new accelerations
418 t_enegrad=t_enegrad+MPI_Wtime()-tt0
420 t_enegrad=t_enegrad+tcpu()-tt0
422 c Determine maximum acceleration and scale down the timestep if needed
424 amax=amax/(itime_scal+1)**2
425 call predict_edrift(epdrift)
426 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
427 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
429 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
431 itime_scal=itime_scal+ifac_time
432 c fac_time=dmin1(damax/amax,0.5d0)
433 fac_time=0.5d0**ifac_time
434 d_time=d_time*fac_time
435 if (lang.eq.2 .or. lang.eq.3) then
437 c write (iout,*) "Calling sd_verlet_setup: 1"
438 c Rescale the stochastic forces and recalculate or restore
439 c the matrices of tinker integrator
440 if (itime_scal.gt.maxflag_stoch) then
441 if (large) write (iout,'(a,i5,a)')
442 & "Calculate matrices for stochastic step;",
443 & " itime_scal ",itime_scal
445 call sd_verlet_p_setup
447 call sd_verlet_ciccotti_setup
449 write (iout,'(2a,i3,a,i3,1h.)')
450 & "Warning: cannot store matrices for stochastic",
451 & " integration because the index",itime_scal,
452 & " is greater than",maxflag_stoch
453 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
454 & " integration Langevin algorithm for better efficiency."
455 else if (flag_stoch(itime_scal)) then
456 if (large) write (iout,'(a,i5,a,l1)')
457 & "Restore matrices for stochastic step; itime_scal ",
458 & itime_scal," flag ",flag_stoch(itime_scal)
461 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
462 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
463 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
464 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
465 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
466 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
470 if (large) write (iout,'(2a,i5,a,l1)')
471 & "Calculate & store matrices for stochastic step;",
472 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
474 call sd_verlet_p_setup
476 call sd_verlet_ciccotti_setup
478 flag_stoch(ifac_time)=.true.
481 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
482 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
483 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
484 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
485 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
486 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
490 fac_time=1.0d0/dsqrt(fac_time)
492 stochforcvec(i)=fac_time*stochforcvec(i)
495 else if (lang.eq.1) then
496 c Rescale the accelerations due to stochastic forces
497 fac_time=1.0d0/dsqrt(fac_time)
499 d_as_work(i)=d_as_work(i)*fac_time
502 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
503 & "itime",itime," Timestep scaled down to ",
504 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
506 c Second step of the velocity Verlet algorithm
511 else if (lang.eq.3) then
513 call sd_verlet2_ciccotti
515 else if (lang.eq.1) then
525 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
528 d_t(j,i)=d_t(j,i)*scale_nh
533 if (rattle) call rattle2
535 if (d_time.ne.d_time0) then
538 if (lang.eq.2 .or. lang.eq.3) then
539 if (large) write (iout,'(a)')
540 & "Restore original matrices for stochastic step"
541 c write (iout,*) "Calling sd_verlet_setup: 2"
542 c Restore the matrices of tinker integrator if the time step has been restored
545 pfric_mat(i,j)=pfric0_mat(i,j,0)
546 afric_mat(i,j)=afric0_mat(i,j,0)
547 vfric_mat(i,j)=vfric0_mat(i,j,0)
548 prand_mat(i,j)=prand0_mat(i,j,0)
549 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
550 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
559 c Calculate the kinetic and the total energy and the kinetic temperature
560 if (tnp .or. tnp1) then
563 d_t_old(j,i)=d_t(j,i)
564 d_t(j,i)=d_t(j,i)/s_np
572 c write (iout,*) "step",itime," EK",EK," EK1",EK1
574 c Couple the system to Berendsen bath if needed
575 if (tbf .and. lang.eq.0) then
578 kinetic_T=2.0d0/(dimen3*Rb)*EK
579 c Backup the coordinates, velocities, and accelerations
583 if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
584 d_a_old(j,i)=d_a(j,i)
588 if (mod(itime,ntwe).eq.0) then
590 if(tnp .or. tnp1) then
591 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
593 cd write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
594 cd & ,EK+potE+pi_np**2/(2*Q_np)+dimen3*0.001986d0*t_bath*log(s_np)
595 write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
599 HNose1=Hnose_nh(EK,potE)
601 write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
617 if (itype(i).ne.10) then
621 vtnp(itnp)=d_t(j,inres)
626 c Transform velocities from UNRES coordinate space to cartesian and Gvec
633 vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
634 vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
636 vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
640 write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
643 write (iout,*) "Velocities, step 2"
645 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
646 & (d_t(j,i+nres),j=1,3)
653 c-------------------------------------------------------------------------------
654 subroutine RESPA_step(itime)
655 c-------------------------------------------------------------------------------
656 c Perform a single RESPA step.
657 c-------------------------------------------------------------------------------
658 implicit real*8 (a-h,o-z)
662 integer IERROR,ERRCODE
664 include 'COMMON.SETUP'
665 include 'COMMON.CONTROL'
669 include 'COMMON.LANGEVIN'
671 include 'COMMON.LANGEVIN.lang0'
673 include 'COMMON.CHAIN'
674 include 'COMMON.DERIV'
676 include 'COMMON.LOCAL'
677 include 'COMMON.INTERACT'
678 include 'COMMON.IOUNITS'
679 include 'COMMON.NAMES'
680 include 'COMMON.TIME1'
681 double precision energia_short(0:n_ene),
682 & energia_long(0:n_ene)
683 double precision cm(3),L(3),vcm(3),incr(3)
684 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
685 & d_a_old0(3,0:maxres2)
686 integer ilen,count,rstcount
689 integer maxcount_scale /10/
690 common /gucio/ cm,energia_short
691 double precision stochforcvec(MAXRES6)
692 common /stochcalc/ stochforcvec
695 double precision vtnp(maxres6), vtnp_(maxres6), vtnp_a(maxres6)
696 common /cipiszcze/ itt
699 if (large.and. mod(itime,ntwe).eq.0) then
700 write (iout,*) "***************** RESPA itime",itime
701 write (iout,*) "Cartesian and internal coordinates: step 0"
703 call pdbout(0.0d0,"cipiszcze",iout)
705 write (iout,*) "Accelerations from long-range forces"
707 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
708 & (d_a(j,i+nres),j=1,3)
710 write (iout,*) "Velocities, step 0"
712 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
713 & (d_t(j,i+nres),j=1,3)
718 c Perform the initial RESPA step (increment velocities)
719 c write (iout,*) "*********************** RESPA ini"
725 if (tnh.and..not.xiresp) then
726 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
729 d_t(j,i)=d_t(j,i)*scale_nh
736 cd if(tnp .or. tnp1) then
737 cd write (iout,'(a,3f)') "EE1 NP S, pi",totT, s_np, pi_np
741 if (mod(itime,ntwe).eq.0 .and. large) then
742 write (iout,*) "Velocities, end"
744 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
745 & (d_t(j,i+nres),j=1,3)
749 c Compute the short-range forces
755 C 7/2/2009 commented out
757 c call etotal_short(energia_short)
758 if (tnp.or.tnp1) potE=energia_short(0)
761 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
764 d_a(j,i)=d_a_short(j,i)
768 if (large.and. mod(itime,ntwe).eq.0) then
769 write (iout,*) "energia_short",energia_short(0)
770 write (iout,*) "Accelerations from short-range forces"
772 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
773 & (d_a(j,i+nres),j=1,3)
778 t_enegrad=t_enegrad+MPI_Wtime()-tt0
780 t_enegrad=t_enegrad+tcpu()-tt0
785 if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
786 d_a_old(j,i)=d_a(j,i)
789 c 6/30/08 A-MTS: attempt at increasing the split number
792 dc_old0(j,i)=dc_old(j,i)
793 d_t_old0(j,i)=d_t_old(j,i)
794 d_a_old0(j,i)=d_a_old(j,i)
797 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
798 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
805 c write (iout,*) "itime",itime," ntime_split",ntime_split
806 c Split the time step
807 d_time=d_time0/ntime_split
808 c Perform the short-range RESPA steps (velocity Verlet increments of
809 c positions and velocities using short-range forces)
810 c write (iout,*) "*********************** RESPA split"
811 do itsplit=1,ntime_split
814 else if (lang.eq.2 .or. lang.eq.3) then
816 call stochastic_force(stochforcvec)
819 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
821 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
826 c First step of the velocity Verlet algorithm
831 else if (lang.eq.3) then
833 call sd_verlet1_ciccotti
835 else if (lang.eq.1) then
838 call tnp1_respa_i_step1
840 call tnp_respa_i_step1
842 if (tnh.and.xiresp) then
844 call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
847 d_t_old(j,i)=d_t_old(j,i)*scale_nh
850 cd write(iout,*) "SSS1",itsplit,EK,scale_nh
854 c Build the chain from the newly calculated coordinates
856 if (rattle) call rattle1
858 if (large.and. mod(itime,ntwe).eq.0) then
859 write (iout,*) "***** ITSPLIT",itsplit
860 write (iout,*) "Cartesian and internal coordinates: step 1"
861 call pdbout(0.0d0,"cipiszcze",iout)
864 write (iout,*) "Velocities, step 1"
866 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
867 & (d_t(j,i+nres),j=1,3)
876 c Calculate energy and forces
878 call etotal_short(energia_short)
880 potE=energia_short(0)
883 t_eshort=t_eshort+MPI_Wtime()-tt0
885 t_eshort=t_eshort+tcpu()-tt0
889 c Get the new accelerations
891 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
894 d_a_short(j,i)=d_a(j,i)
898 if (large.and. mod(itime,ntwe).eq.0) then
899 write (iout,*)"energia_short",energia_short(0)
900 write (iout,*) "Accelerations from short-range forces"
902 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
903 & (d_a(j,i+nres),j=1,3)
908 c Determine maximum acceleration and scale down the timestep if needed
910 amax=amax/ntime_split**2
911 call predict_edrift(epdrift)
912 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
913 & write (iout,*) "amax",amax," damax",damax,
914 & " epdrift",epdrift," epdriftmax",epdriftmax
915 c Exit loop and try with increased split number if the change of
916 c acceleration is too big
917 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
918 if (ntime_split.lt.maxtime_split) then
920 ntime_split=ntime_split*2
923 dc_old(j,i)=dc_old0(j,i)
924 d_t_old(j,i)=d_t_old0(j,i)
925 d_a_old(j,i)=d_a_old0(j,i)
928 write (iout,*) "acceleration/energy drift too large",amax,
929 & epdrift," split increased to ",ntime_split," itime",itime,
934 & "Uh-hu. Bumpy landscape. Maximum splitting number",
936 & " already reached!!! Trying to carry on!"
940 t_enegrad=t_enegrad+MPI_Wtime()-tt0
942 t_enegrad=t_enegrad+tcpu()-tt0
944 c Second step of the velocity Verlet algorithm
949 else if (lang.eq.3) then
951 call sd_verlet2_ciccotti
953 else if (lang.eq.1) then
956 call tnp1_respa_i_step2
958 call tnp_respa_i_step2
961 if (tnh.and.xiresp) then
963 call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
966 d_t(j,i)=d_t(j,i)*scale_nh
969 cd write(iout,*) "SSS2",itsplit,EK,scale_nh
972 if (rattle) call rattle2
973 c Backup the coordinates, velocities, and accelerations
974 if (tnp .or. tnp1) then
977 d_t_old(j,i)=d_t(j,i)
978 if (tnp) d_t(j,i)=d_t(j,i)/s_np
979 if (tnp1) d_t(j,i)=d_t(j,i)/s_np
987 if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
988 d_a_old(j,i)=d_a(j,i)
995 c Restore the time step
997 c Compute long-range forces
1004 call etotal_long(energia_long)
1005 E_long=energia_long(0)
1006 potE=energia_short(0)+energia_long(0)
1009 t_elong=t_elong+MPI_Wtime()-tt0
1011 t_elong=t_elong+tcpu()-tt0
1017 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1019 t_enegrad=t_enegrad+tcpu()-tt0
1021 c Compute accelerations from long-range forces
1023 if (large.and. mod(itime,ntwe).eq.0) then
1024 write (iout,*) "energia_long",energia_long(0)
1025 write (iout,*) "Cartesian and internal coordinates: step 2"
1027 call pdbout(0.0d0,"cipiszcze",iout)
1029 write (iout,*) "Accelerations from long-range forces"
1031 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1032 & (d_a(j,i+nres),j=1,3)
1034 write (iout,*) "Velocities, step 2"
1036 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1037 & (d_t(j,i+nres),j=1,3)
1041 c Compute the final RESPA step (increment velocities)
1042 c write (iout,*) "*********************** RESPA fin"
1044 call tnp_respa_step2
1046 call tnp_respa_step2
1049 if (tnh.and..not.xiresp) then
1051 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
1054 d_t(j,i)=d_t(j,i)*scale_nh
1060 if (tnp .or. tnp1) then
1063 d_t(j,i)=d_t_old(j,i)/s_np
1068 c Compute the complete potential energy
1070 potEcomp(i)=energia_short(i)+energia_long(i)
1072 potE=potEcomp(0)-potEcomp(20)
1073 c potE=energia_short(0)+energia_long(0)
1075 c Calculate the kinetic and the total energy and the kinetic temperature
1078 c Couple the system to Berendsen bath if needed
1079 if (tbf .and. lang.eq.0) then
1082 kinetic_T=2.0d0/(dimen3*Rb)*EK
1083 c Backup the coordinates, velocities, and accelerations
1085 if (mod(itime,ntwe).eq.0 .and. large) then
1086 write (iout,*) "Velocities, end"
1088 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1089 & (d_t(j,i+nres),j=1,3)
1093 if (mod(itime,ntwe).eq.0) then
1095 if(tnp .or. tnp1) then
1097 write (iout,'(a3,7f)') "TTT",EK,s_np,potE,pi_np,Csplit,
1098 & E_long,energia_short(0)
1100 write (iout,'(a3,7f20.10)') "TTT",EK,s_np,potE,pi_np,Csplit,
1101 & E_long,energia_short(0)
1103 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
1105 cd write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
1106 cd & ,EK+potE+pi_np**2/(2*Q_np)+dimen3*0.001986d0*t_bath*log(s_np)
1107 write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
1108 cd write (iout,'(a,3f)') "EE2 NP S, pi",totT, s_np, pi_np
1112 HNose1=Hnose_nh(EK,potE)
1114 write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
1131 if (itype(i).ne.10) then
1135 vtnp(itnp)=d_t(j,inres)
1140 c Transform velocities from UNRES coordinate space to cartesian and Gvec
1147 vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
1148 vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
1150 vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
1154 write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
1164 c---------------------------------------------------------------------
1165 subroutine RESPA_vel
1166 c First and last RESPA step (incrementing velocities using long-range
1168 implicit real*8 (a-h,o-z)
1169 include 'DIMENSIONS'
1170 include 'COMMON.CONTROL'
1171 include 'COMMON.VAR'
1173 include 'COMMON.CHAIN'
1174 include 'COMMON.DERIV'
1175 include 'COMMON.GEO'
1176 include 'COMMON.LOCAL'
1177 include 'COMMON.INTERACT'
1178 include 'COMMON.IOUNITS'
1179 include 'COMMON.NAMES'
1181 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1185 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1189 if (itype(i).ne.10) then
1192 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1198 c-----------------------------------------------------------------
1200 c Applying velocity Verlet algorithm - step 1 to coordinates
1201 implicit real*8 (a-h,o-z)
1202 include 'DIMENSIONS'
1203 include 'COMMON.CONTROL'
1204 include 'COMMON.VAR'
1206 include 'COMMON.CHAIN'
1207 include 'COMMON.DERIV'
1208 include 'COMMON.GEO'
1209 include 'COMMON.LOCAL'
1210 include 'COMMON.INTERACT'
1211 include 'COMMON.IOUNITS'
1212 include 'COMMON.NAMES'
1213 double precision adt,adt2
1216 write (iout,*) "VELVERLET1 START: DC"
1218 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1219 & (dc(j,i+nres),j=1,3)
1223 adt=d_a_old(j,0)*d_time
1225 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1226 d_t_new(j,0)=d_t_old(j,0)+adt2
1227 d_t(j,0)=d_t_old(j,0)+adt
1231 adt=d_a_old(j,i)*d_time
1233 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1234 d_t_new(j,i)=d_t_old(j,i)+adt2
1235 d_t(j,i)=d_t_old(j,i)+adt
1239 if (itype(i).ne.10) then
1242 adt=d_a_old(j,inres)*d_time
1244 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1245 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1246 d_t(j,inres)=d_t_old(j,inres)+adt
1251 write (iout,*) "VELVERLET1 END: DC"
1253 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1254 & (dc(j,i+nres),j=1,3)
1259 c---------------------------------------------------------------------
1261 c Step 2 of the velocity Verlet algorithm: update velocities
1262 implicit real*8 (a-h,o-z)
1263 include 'DIMENSIONS'
1264 include 'COMMON.CONTROL'
1265 include 'COMMON.VAR'
1267 include 'COMMON.CHAIN'
1268 include 'COMMON.DERIV'
1269 include 'COMMON.GEO'
1270 include 'COMMON.LOCAL'
1271 include 'COMMON.INTERACT'
1272 include 'COMMON.IOUNITS'
1273 include 'COMMON.NAMES'
1275 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1279 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1283 if (itype(i).ne.10) then
1286 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1292 c-----------------------------------------------------------------
1293 subroutine sddir_precalc
1294 c Applying velocity Verlet algorithm - step 1 to coordinates
1295 implicit real*8 (a-h,o-z)
1296 include 'DIMENSIONS'
1300 include 'COMMON.CONTROL'
1301 include 'COMMON.VAR'
1304 include 'COMMON.LANGEVIN'
1306 include 'COMMON.LANGEVIN.lang0'
1308 include 'COMMON.CHAIN'
1309 include 'COMMON.DERIV'
1310 include 'COMMON.GEO'
1311 include 'COMMON.LOCAL'
1312 include 'COMMON.INTERACT'
1313 include 'COMMON.IOUNITS'
1314 include 'COMMON.NAMES'
1315 include 'COMMON.TIME1'
1316 double precision stochforcvec(MAXRES6)
1317 common /stochcalc/ stochforcvec
1319 c Compute friction and stochastic forces
1328 time_fric=time_fric+MPI_Wtime()-time00
1331 time_fric=time_fric+tcpu()-time00
1334 call stochastic_force(stochforcvec)
1336 time_stoch=time_stoch+MPI_Wtime()-time00
1338 time_stoch=time_stoch+tcpu()-time00
1341 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1342 c forces (d_as_work)
1344 call ginv_mult(fric_work, d_af_work)
1345 call ginv_mult(stochforcvec, d_as_work)
1348 c---------------------------------------------------------------------
1349 subroutine sddir_verlet1
1350 c Applying velocity Verlet algorithm - step 1 to velocities
1351 implicit real*8 (a-h,o-z)
1352 include 'DIMENSIONS'
1353 include 'COMMON.CONTROL'
1354 include 'COMMON.VAR'
1357 include 'COMMON.LANGEVIN'
1359 include 'COMMON.LANGEVIN.lang0'
1361 include 'COMMON.CHAIN'
1362 include 'COMMON.DERIV'
1363 include 'COMMON.GEO'
1364 include 'COMMON.LOCAL'
1365 include 'COMMON.INTERACT'
1366 include 'COMMON.IOUNITS'
1367 include 'COMMON.NAMES'
1368 c Revised 3/31/05 AL: correlation between random contributions to
1369 c position and velocity increments included.
1370 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1371 double precision adt,adt2
1373 c Add the contribution from BOTH friction and stochastic force to the
1374 c coordinates, but ONLY the contribution from the friction forces to velocities
1377 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1378 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1379 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1380 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1381 d_t(j,0)=d_t_old(j,0)+adt
1386 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1387 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1388 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1389 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1390 d_t(j,i)=d_t_old(j,i)+adt
1395 if (itype(i).ne.10) then
1398 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1399 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1400 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1401 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1402 d_t(j,inres)=d_t_old(j,inres)+adt
1409 c---------------------------------------------------------------------
1410 subroutine sddir_verlet2
1411 c Calculating the adjusted velocities for accelerations
1412 implicit real*8 (a-h,o-z)
1413 include 'DIMENSIONS'
1414 include 'COMMON.CONTROL'
1415 include 'COMMON.VAR'
1418 include 'COMMON.LANGEVIN'
1420 include 'COMMON.LANGEVIN.lang0'
1422 include 'COMMON.CHAIN'
1423 include 'COMMON.DERIV'
1424 include 'COMMON.GEO'
1425 include 'COMMON.LOCAL'
1426 include 'COMMON.INTERACT'
1427 include 'COMMON.IOUNITS'
1428 include 'COMMON.NAMES'
1429 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1430 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1431 c Revised 3/31/05 AL: correlation between random contributions to
1432 c position and velocity increments included.
1433 c The correlation coefficients are calculated at low-friction limit.
1434 c Also, friction forces are now not calculated with new velocities.
1436 c call friction_force
1437 call stochastic_force(stochforcvec)
1439 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1440 c forces (d_as_work)
1442 call ginv_mult(stochforcvec, d_as_work1)
1448 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1449 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1454 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1455 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1460 if (itype(i).ne.10) then
1463 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1464 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1465 & +cos60*d_as_work1(ind+j))*d_time
1472 c---------------------------------------------------------------------
1473 subroutine max_accel
1475 c Find the maximum difference in the accelerations of the the sites
1476 c at the beginning and the end of the time step.
1478 implicit real*8 (a-h,o-z)
1479 include 'DIMENSIONS'
1480 include 'COMMON.CONTROL'
1481 include 'COMMON.VAR'
1483 include 'COMMON.CHAIN'
1484 include 'COMMON.DERIV'
1485 include 'COMMON.GEO'
1486 include 'COMMON.LOCAL'
1487 include 'COMMON.INTERACT'
1488 include 'COMMON.IOUNITS'
1489 double precision aux(3),accel(3),accel_old(3),dacc
1491 c aux(j)=d_a(j,0)-d_a_old(j,0)
1492 accel_old(j)=d_a_old(j,0)
1499 c 7/3/08 changed to asymmetric difference
1501 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1502 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1503 accel(j)=accel(j)+0.5d0*d_a(j,i)
1504 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1505 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1506 dacc=dabs(accel(j)-accel_old(j))
1507 if (dacc.gt.amax) amax=dacc
1515 accel_old(j)=d_a_old(j,0)
1520 accel_old(j)=accel_old(j)+d_a_old(j,1)
1521 accel(j)=accel(j)+d_a(j,1)
1525 if (itype(i).ne.10) then
1527 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1528 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1529 accel(j)=accel(j)+d_a(j,i+nres)
1533 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1534 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1535 dacc=dabs(accel(j)-accel_old(j))
1536 if (dacc.gt.amax) amax=dacc
1540 accel_old(j)=accel_old(j)+d_a_old(j,i)
1541 accel(j)=accel(j)+d_a(j,i)
1542 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1547 c---------------------------------------------------------------------
1548 subroutine predict_edrift(epdrift)
1550 c Predict the drift of the potential energy
1552 implicit real*8 (a-h,o-z)
1553 include 'DIMENSIONS'
1554 include 'COMMON.CONTROL'
1555 include 'COMMON.VAR'
1557 include 'COMMON.CHAIN'
1558 include 'COMMON.DERIV'
1559 include 'COMMON.GEO'
1560 include 'COMMON.LOCAL'
1561 include 'COMMON.INTERACT'
1562 include 'COMMON.IOUNITS'
1563 include 'COMMON.MUCA'
1564 double precision epdrift,epdriftij
1565 c Drift of the potential energy
1571 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1572 if (lmuca) epdriftij=epdriftij*factor
1573 c write (iout,*) "back",i,j,epdriftij
1574 if (epdriftij.gt.epdrift) epdrift=epdriftij
1578 if (itype(i).ne.10) then
1581 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1582 if (lmuca) epdriftij=epdriftij*factor
1583 c write (iout,*) "side",i,j,epdriftij
1584 if (epdriftij.gt.epdrift) epdrift=epdriftij
1588 epdrift=0.5d0*epdrift*d_time*d_time
1589 c write (iout,*) "epdrift",epdrift
1592 c-----------------------------------------------------------------------
1593 subroutine verlet_bath
1595 c Coupling to the thermostat by using the Berendsen algorithm
1597 implicit real*8 (a-h,o-z)
1598 include 'DIMENSIONS'
1599 include 'COMMON.CONTROL'
1600 include 'COMMON.VAR'
1602 include 'COMMON.CHAIN'
1603 include 'COMMON.DERIV'
1604 include 'COMMON.GEO'
1605 include 'COMMON.LOCAL'
1606 include 'COMMON.INTERACT'
1607 include 'COMMON.IOUNITS'
1608 include 'COMMON.NAMES'
1609 double precision T_half,fact
1611 T_half=2.0d0/(dimen3*Rb)*EK
1612 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1613 c write(iout,*) "T_half", T_half
1614 c write(iout,*) "EK", EK
1615 c write(iout,*) "fact", fact
1617 d_t(j,0)=fact*d_t(j,0)
1621 d_t(j,i)=fact*d_t(j,i)
1625 if (itype(i).ne.10) then
1628 d_t(j,inres)=fact*d_t(j,inres)
1634 c---------------------------------------------------------
1636 c Set up the initial conditions of a MD simulation
1637 implicit real*8 (a-h,o-z)
1638 include 'DIMENSIONS'
1642 integer IERROR,ERRCODE
1644 include 'COMMON.SETUP'
1645 include 'COMMON.CONTROL'
1646 include 'COMMON.VAR'
1649 include 'COMMON.LANGEVIN'
1651 include 'COMMON.LANGEVIN.lang0'
1653 include 'COMMON.CHAIN'
1654 include 'COMMON.DERIV'
1655 include 'COMMON.GEO'
1656 include 'COMMON.LOCAL'
1657 include 'COMMON.INTERACT'
1658 include 'COMMON.IOUNITS'
1659 include 'COMMON.NAMES'
1660 include 'COMMON.REMD'
1661 real*8 energia_long(0:n_ene),
1662 & energia_short(0:n_ene),vcm(3),incr(3),E_short
1663 double precision cm(3),L(3),xv,sigv,lowb,highb
1664 double precision varia(maxvar)
1672 c write(iout,*) "d_time", d_time
1673 c Compute the standard deviations of stochastic forces for Langevin dynamics
1674 c if the friction coefficients do not depend on surface area
1675 if (lang.gt.0 .and. .not.surfarea) then
1677 stdforcp(i)=stdfp*dsqrt(gamp)
1680 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
1683 c Open the pdb file for snapshotshots
1686 if (ilen(tmpdir).gt.0)
1687 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1688 & liczba(:ilen(liczba))//".pdb")
1690 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1694 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1695 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1696 & liczba(:ilen(liczba))//".x")
1697 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1700 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1701 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1702 & liczba(:ilen(liczba))//".cx")
1703 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1709 if (ilen(tmpdir).gt.0)
1710 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1711 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1713 if (ilen(tmpdir).gt.0)
1714 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1715 cartname=prefix(:ilen(prefix))//"_MD.cx"
1719 write (qstr,'(256(1h ))')
1722 iq = qinfrag(i,iset)*10
1723 iw = wfrag(i,iset)/100
1725 if(me.eq.king.or..not.out1file)
1726 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1727 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1732 iq = qinpair(i,iset)*10
1733 iw = wpair(i,iset)/100
1735 if(me.eq.king.or..not.out1file)
1736 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1737 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1741 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1743 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1745 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1747 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1751 if (restart1file) then
1753 & inquire(file=mremd_rst_name,exist=file_exist)
1754 write (*,*) me," Before broadcast: file_exist",file_exist
1756 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1758 write (*,*) me," After broadcast: file_exist",file_exist
1760 c inquire(file=mremd_rst_name,exist=file_exist)
1761 if(me.eq.king.or..not.out1file)
1762 & write(iout,*) "Initial state read by master and distributed"
1764 if (ilen(tmpdir).gt.0)
1765 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1766 & //liczba(:ilen(liczba))//'.rst')
1767 inquire(file=rest2name,exist=file_exist)
1770 if(.not.restart1file) then
1771 if(me.eq.king.or..not.out1file)
1772 & write(iout,*) "Initial state will be read from file ",
1773 & rest2name(:ilen(rest2name))
1776 call rescale_weights(t_bath)
1778 if(me.eq.king.or..not.out1file)then
1779 if (restart1file) then
1780 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1783 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1786 write(iout,*) "Initial velocities randomly generated"
1792 c Generate initial velocities
1793 if(me.eq.king.or..not.out1file)
1794 & write(iout,*) "Initial velocities randomly generated"
1798 c rest2name = prefix(:ilen(prefix))//'.rst'
1799 if(me.eq.king.or..not.out1file)then
1800 write (iout,*) "Initial velocities"
1802 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1803 & (d_t(j,i+nres),j=1,3)
1805 c Zeroing the total angular momentum of the system
1806 write(iout,*) "Calling the zero-angular
1807 & momentum subroutine"
1810 c Getting the potential energy and forces and velocities and accelerations
1812 c write (iout,*) "velocity of the center of the mass:"
1813 c write (iout,*) (vcm(j),j=1,3)
1815 d_t(j,0)=d_t(j,0)-vcm(j)
1817 c Removing the velocity of the center of mass
1819 if(me.eq.king.or..not.out1file)then
1820 write (iout,*) "vcm right after adjustment:"
1821 write (iout,*) (vcm(j),j=1,3)
1825 if(iranconf.ne.0) then
1827 print *, 'Calling OVERLAP_SC'
1828 call overlap_sc(fail)
1832 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1833 print *,'SC_move',nft_sc,etot
1834 if(me.eq.king.or..not.out1file)
1835 & write(iout,*) 'SC_move',nft_sc,etot
1839 print *, 'Calling MINIM_DC'
1840 call minim_dc(etot,iretcode,nfun)
1842 call geom_to_var(nvar,varia)
1843 print *,'Calling MINIMIZE.'
1844 call minimize(etot,varia,iretcode,nfun)
1845 call var_to_geom(nvar,varia)
1847 if(me.eq.king.or..not.out1file)
1848 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1851 call chainbuild_cart
1854 call verlet_bath(EK)
1856 kinetic_T=2.0d0/(dimen3*Rb)*EK
1857 if(me.eq.king.or..not.out1file)then
1867 call etotal(potEcomp)
1870 t_etotal=t_etotal+MPI_Wtime()-tt0
1872 t_etotal=t_etotal+tcpu()-tt0
1877 if(tnp .or. tnp1) then
1880 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
1882 write(iout,*) 'H0= ',H0
1886 HNose1=Hnose_nh(EK,potE)
1888 write (iout,*) 'H0= ',H0
1895 if (amax*d_time .gt. dvmax) then
1896 d_time=d_time*dvmax/amax
1897 if(me.eq.king.or..not.out1file) write (iout,*)
1898 & "Time step reduced to",d_time,
1899 & " because of too large initial acceleration."
1901 if(me.eq.king.or..not.out1file)then
1902 write(iout,*) "Potential energy and its components"
1903 call enerprint(potEcomp)
1904 c write(iout,*) (potEcomp(i),i=0,n_ene)
1906 potE=potEcomp(0)-potEcomp(20)
1909 if (ntwe.ne.0) call statout(itime)
1910 if(me.eq.king.or..not.out1file)
1911 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1912 & " Kinetic energy",EK," potential energy",potE,
1913 & " total energy",totE," maximum acceleration ",
1916 write (iout,*) "Initial coordinates"
1918 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1919 & (c(j,i+nres),j=1,3)
1921 write (iout,*) "Initial dC"
1923 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1924 & (dc(j,i+nres),j=1,3)
1926 write (iout,*) "Initial velocities"
1927 write (iout,"(13x,' backbone ',23x,' side chain')")
1929 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1930 & (d_t(j,i+nres),j=1,3)
1932 write (iout,*) "Initial accelerations"
1934 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1935 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1936 & (d_a(j,i+nres),j=1,3)
1942 d_t_old(j,i)=d_t(j,i)
1943 d_a_old(j,i)=d_a(j,i)
1945 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1954 call etotal_short(energia_short)
1957 t_eshort=t_eshort+MPI_Wtime()-tt0
1959 t_eshort=t_eshort+tcpu()-tt0
1963 if(tnp .or. tnp1) then
1964 E_short=energia_short(0)
1965 HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen3)
1968 c_new_var_csplit Csplit=H0-E_long
1969 c Csplit = H0-energia_short(0)
1970 write(iout,*) 'Csplit= ',Csplit
1976 if(.not.out1file .and. large) then
1977 write (iout,*) "energia_long",energia_long(0),
1978 & " energia_short",energia_short(0),
1979 & " total",energia_long(0)+energia_short(0)
1980 write (iout,*) "Initial fast-force accelerations"
1982 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1983 & (d_a(j,i+nres),j=1,3)
1986 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1989 d_a_short(j,i)=d_a(j,i)
1998 call etotal_long(energia_long)
2001 t_elong=t_elong+MPI_Wtime()-tt0
2003 t_elong=t_elong+tcpu()-tt0
2008 if(.not.out1file .and. large) then
2009 write (iout,*) "energia_long",energia_long(0)
2010 write (iout,*) "Initial slow-force accelerations"
2012 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2013 & (d_a(j,i+nres),j=1,3)
2017 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2019 t_enegrad=t_enegrad+tcpu()-tt0
2024 c-----------------------------------------------------------
2025 subroutine random_vel
2026 implicit real*8 (a-h,o-z)
2027 include 'DIMENSIONS'
2028 include 'COMMON.CONTROL'
2029 include 'COMMON.VAR'
2032 include 'COMMON.LANGEVIN'
2034 include 'COMMON.LANGEVIN.lang0'
2036 include 'COMMON.CHAIN'
2037 include 'COMMON.DERIV'
2038 include 'COMMON.GEO'
2039 include 'COMMON.LOCAL'
2040 include 'COMMON.INTERACT'
2041 include 'COMMON.IOUNITS'
2042 include 'COMMON.NAMES'
2043 include 'COMMON.TIME1'
2044 double precision xv,sigv,lowb,highb
2045 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2046 c First generate velocities in the eigenspace of the G matrix
2047 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2049 c write (iout,*) "RANDOM_VEL dimen",dimen
2055 sigv=dsqrt((Rb*t_bath)/geigen(i))
2058 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2059 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2060 c & " d_t_work_new",d_t_work_new(ii)
2069 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2072 c write (iout,*) "Ek from eigenvectors",Ek1
2074 c Transform velocities to UNRES coordinate space
2080 d_t_work(ind)=d_t_work(ind)
2081 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2083 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2087 c Transfer to the d_t vector
2089 d_t(j,0)=d_t_work(j)
2095 d_t(j,i)=d_t_work(ind)
2099 c write (iout,*) "d_t",i,(d_t(j,i),j=1,3)
2102 if (itype(i).ne.10) then
2105 d_t(j,i+nres)=d_t_work(ind)
2110 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2111 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2116 c-----------------------------------------------------------
2117 subroutine sd_verlet_p_setup
2118 c Sets up the parameters of stochastic Verlet algorithm
2119 implicit real*8 (a-h,o-z)
2120 include 'DIMENSIONS'
2124 include 'COMMON.CONTROL'
2125 include 'COMMON.VAR'
2128 include 'COMMON.LANGEVIN'
2130 include 'COMMON.LANGEVIN.lang0'
2132 include 'COMMON.CHAIN'
2133 include 'COMMON.DERIV'
2134 include 'COMMON.GEO'
2135 include 'COMMON.LOCAL'
2136 include 'COMMON.INTERACT'
2137 include 'COMMON.IOUNITS'
2138 include 'COMMON.NAMES'
2139 include 'COMMON.TIME1'
2140 double precision emgdt(MAXRES6),
2141 & pterm,vterm,rho,rhoc,vsig,
2142 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2143 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2144 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2145 logical lprn /.false./
2146 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2147 double precision ktm
2154 c AL 8/17/04 Code adapted from tinker
2156 c Get the frictional and random terms for stochastic dynamics in the
2157 c eigenspace of mass-scaled UNRES friction matrix
2160 gdt = fricgam(i) * d_time
2162 c Stochastic dynamics reduces to simple MD for zero friction
2164 if (gdt .le. zero) then
2165 pfric_vec(i) = 1.0d0
2166 vfric_vec(i) = d_time
2167 afric_vec(i) = 0.5d0 * d_time * d_time
2168 prand_vec(i) = 0.0d0
2169 vrand_vec1(i) = 0.0d0
2170 vrand_vec2(i) = 0.0d0
2172 c Analytical expressions when friction coefficient is large
2175 if (gdt .ge. gdt_radius) then
2178 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2179 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2180 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2181 vterm = 1.0d0 - egdt**2
2182 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2184 c Use series expansions when friction coefficient is small
2195 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2196 & - gdt5/120.0d0 + gdt6/720.0d0
2197 & - gdt7/5040.0d0 + gdt8/40320.0d0
2198 & - gdt9/362880.0d0) / fricgam(i)**2
2199 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2200 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2201 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2202 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2203 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2204 & + 127.0d0*gdt9/90720.0d0
2205 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2206 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2207 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2208 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2209 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2210 & - 17.0d0*gdt2/1280.0d0
2211 & + 17.0d0*gdt3/6144.0d0
2212 & + 40967.0d0*gdt4/34406400.0d0
2213 & - 57203.0d0*gdt5/275251200.0d0
2214 & - 1429487.0d0*gdt6/13212057600.0d0)
2217 c Compute the scaling factors of random terms for the nonzero friction case
2219 ktm = 0.5d0*d_time/fricgam(i)
2220 psig = dsqrt(ktm*pterm) / fricgam(i)
2221 vsig = dsqrt(ktm*vterm)
2222 rhoc = dsqrt(1.0d0 - rho*rho)
2224 vrand_vec1(i) = vsig * rho
2225 vrand_vec2(i) = vsig * rhoc
2230 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2233 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2234 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2238 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2241 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2242 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2243 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2244 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2245 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2246 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2249 t_sdsetup=t_sdsetup+MPI_Wtime()
2251 t_sdsetup=t_sdsetup+tcpu()-tt0
2255 c-------------------------------------------------------------
2256 subroutine eigtransf1(n,ndim,ab,d,c)
2259 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2265 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2271 c-------------------------------------------------------------
2272 subroutine eigtransf(n,ndim,a,b,d,c)
2275 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2281 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2287 c-------------------------------------------------------------
2288 subroutine sd_verlet1
2289 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2290 implicit real*8 (a-h,o-z)
2291 include 'DIMENSIONS'
2292 include 'COMMON.CONTROL'
2293 include 'COMMON.VAR'
2296 include 'COMMON.LANGEVIN'
2298 include 'COMMON.LANGEVIN.lang0'
2300 include 'COMMON.CHAIN'
2301 include 'COMMON.DERIV'
2302 include 'COMMON.GEO'
2303 include 'COMMON.LOCAL'
2304 include 'COMMON.INTERACT'
2305 include 'COMMON.IOUNITS'
2306 include 'COMMON.NAMES'
2307 double precision stochforcvec(MAXRES6)
2308 common /stochcalc/ stochforcvec
2309 logical lprn /.false./
2311 c write (iout,*) "dc_old"
2313 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2314 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2317 dc_work(j)=dc_old(j,0)
2318 d_t_work(j)=d_t_old(j,0)
2319 d_a_work(j)=d_a_old(j,0)
2324 dc_work(ind+j)=dc_old(j,i)
2325 d_t_work(ind+j)=d_t_old(j,i)
2326 d_a_work(ind+j)=d_a_old(j,i)
2331 if (itype(i).ne.10) then
2333 dc_work(ind+j)=dc_old(j,i+nres)
2334 d_t_work(ind+j)=d_t_old(j,i+nres)
2335 d_a_work(ind+j)=d_a_old(j,i+nres)
2343 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2347 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2348 & vfric_mat(i,j),afric_mat(i,j),
2349 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2357 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2358 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2359 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2360 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2362 d_t_work_new(i)=ddt1+0.5d0*ddt2
2363 d_t_work(i)=ddt1+ddt2
2368 d_t(j,0)=d_t_work(j)
2373 dc(j,i)=dc_work(ind+j)
2374 d_t(j,i)=d_t_work(ind+j)
2379 if (itype(i).ne.10) then
2382 dc(j,inres)=dc_work(ind+j)
2383 d_t(j,inres)=d_t_work(ind+j)
2390 c--------------------------------------------------------------------------
2391 subroutine sd_verlet2
2392 c Calculating the adjusted velocities for accelerations
2393 implicit real*8 (a-h,o-z)
2394 include 'DIMENSIONS'
2395 include 'COMMON.CONTROL'
2396 include 'COMMON.VAR'
2399 include 'COMMON.LANGEVIN'
2401 include 'COMMON.LANGEVIN.lang0'
2403 include 'COMMON.CHAIN'
2404 include 'COMMON.DERIV'
2405 include 'COMMON.GEO'
2406 include 'COMMON.LOCAL'
2407 include 'COMMON.INTERACT'
2408 include 'COMMON.IOUNITS'
2409 include 'COMMON.NAMES'
2410 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2411 common /stochcalc/ stochforcvec
2413 c Compute the stochastic forces which contribute to velocity change
2415 call stochastic_force(stochforcvecV)
2422 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2423 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2424 & vrand_mat2(i,j)*stochforcvecV(j)
2426 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2430 d_t(j,0)=d_t_work(j)
2435 d_t(j,i)=d_t_work(ind+j)
2440 if (itype(i).ne.10) then
2443 d_t(j,inres)=d_t_work(ind+j)
2450 c-----------------------------------------------------------
2451 subroutine sd_verlet_ciccotti_setup
2452 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2454 implicit real*8 (a-h,o-z)
2455 include 'DIMENSIONS'
2459 include 'COMMON.CONTROL'
2460 include 'COMMON.VAR'
2463 include 'COMMON.LANGEVIN'
2465 include 'COMMON.LANGEVIN.lang0'
2467 include 'COMMON.CHAIN'
2468 include 'COMMON.DERIV'
2469 include 'COMMON.GEO'
2470 include 'COMMON.LOCAL'
2471 include 'COMMON.INTERACT'
2472 include 'COMMON.IOUNITS'
2473 include 'COMMON.NAMES'
2474 include 'COMMON.TIME1'
2475 double precision emgdt(MAXRES6),
2476 & pterm,vterm,rho,rhoc,vsig,
2477 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2478 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2479 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2480 logical lprn /.false./
2481 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2482 double precision ktm
2489 c AL 8/17/04 Code adapted from tinker
2491 c Get the frictional and random terms for stochastic dynamics in the
2492 c eigenspace of mass-scaled UNRES friction matrix
2495 write (iout,*) "i",i," fricgam",fricgam(i)
2496 gdt = fricgam(i) * d_time
2498 c Stochastic dynamics reduces to simple MD for zero friction
2500 if (gdt .le. zero) then
2501 pfric_vec(i) = 1.0d0
2502 vfric_vec(i) = d_time
2503 afric_vec(i) = 0.5d0*d_time*d_time
2504 prand_vec(i) = afric_vec(i)
2505 vrand_vec2(i) = vfric_vec(i)
2507 c Analytical expressions when friction coefficient is large
2512 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2513 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2514 prand_vec(i) = afric_vec(i)
2515 vrand_vec2(i) = vfric_vec(i)
2517 c Compute the scaling factors of random terms for the nonzero friction case
2519 c ktm = 0.5d0*d_time/fricgam(i)
2520 c psig = dsqrt(ktm*pterm) / fricgam(i)
2521 c vsig = dsqrt(ktm*vterm)
2522 c prand_vec(i) = psig*afric_vec(i)
2523 c vrand_vec2(i) = vsig*vfric_vec(i)
2528 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2531 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2532 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2536 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2538 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2539 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2540 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2541 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2542 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2544 t_sdsetup=t_sdsetup+MPI_Wtime()
2546 t_sdsetup=t_sdsetup+tcpu()-tt0
2550 c-------------------------------------------------------------
2551 subroutine sd_verlet1_ciccotti
2552 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2553 implicit real*8 (a-h,o-z)
2554 include 'DIMENSIONS'
2558 include 'COMMON.CONTROL'
2559 include 'COMMON.VAR'
2562 include 'COMMON.LANGEVIN'
2564 include 'COMMON.LANGEVIN.lang0'
2566 include 'COMMON.CHAIN'
2567 include 'COMMON.DERIV'
2568 include 'COMMON.GEO'
2569 include 'COMMON.LOCAL'
2570 include 'COMMON.INTERACT'
2571 include 'COMMON.IOUNITS'
2572 include 'COMMON.NAMES'
2573 double precision stochforcvec(MAXRES6)
2574 common /stochcalc/ stochforcvec
2575 logical lprn /.false./
2577 c write (iout,*) "dc_old"
2579 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2580 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2583 dc_work(j)=dc_old(j,0)
2584 d_t_work(j)=d_t_old(j,0)
2585 d_a_work(j)=d_a_old(j,0)
2590 dc_work(ind+j)=dc_old(j,i)
2591 d_t_work(ind+j)=d_t_old(j,i)
2592 d_a_work(ind+j)=d_a_old(j,i)
2597 if (itype(i).ne.10) then
2599 dc_work(ind+j)=dc_old(j,i+nres)
2600 d_t_work(ind+j)=d_t_old(j,i+nres)
2601 d_a_work(ind+j)=d_a_old(j,i+nres)
2610 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2614 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2615 & vfric_mat(i,j),afric_mat(i,j),
2616 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2624 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2625 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2626 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2627 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2629 d_t_work_new(i)=ddt1+0.5d0*ddt2
2630 d_t_work(i)=ddt1+ddt2
2635 d_t(j,0)=d_t_work(j)
2640 dc(j,i)=dc_work(ind+j)
2641 d_t(j,i)=d_t_work(ind+j)
2646 if (itype(i).ne.10) then
2649 dc(j,inres)=dc_work(ind+j)
2650 d_t(j,inres)=d_t_work(ind+j)
2657 c--------------------------------------------------------------------------
2658 subroutine sd_verlet2_ciccotti
2659 c Calculating the adjusted velocities for accelerations
2660 implicit real*8 (a-h,o-z)
2661 include 'DIMENSIONS'
2662 include 'COMMON.CONTROL'
2663 include 'COMMON.VAR'
2666 include 'COMMON.LANGEVIN'
2668 include 'COMMON.LANGEVIN.lang0'
2670 include 'COMMON.CHAIN'
2671 include 'COMMON.DERIV'
2672 include 'COMMON.GEO'
2673 include 'COMMON.LOCAL'
2674 include 'COMMON.INTERACT'
2675 include 'COMMON.IOUNITS'
2676 include 'COMMON.NAMES'
2677 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2678 common /stochcalc/ stochforcvec
2680 c Compute the stochastic forces which contribute to velocity change
2682 call stochastic_force(stochforcvecV)
2689 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2690 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2691 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2693 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2697 d_t(j,0)=d_t_work(j)
2702 d_t(j,i)=d_t_work(ind+j)
2707 if (itype(i).ne.10) then
2710 d_t(j,inres)=d_t_work(ind+j)
2718 c------------------------------------------------------
2719 double precision function HNose(ek,s,e,pi,Q,t_bath,dimenl)
2721 double precision ek,s,e,pi,Q,t_bath,Rb
2724 HNose=ek+e+pi**2/(2*Q)+dimenl*Rb*t_bath*log(s)
2725 c print '(6f15.5,i5,a2,2f15.5)',ek,s,e,pi,Q,t_bath,dimenl,"--",
2726 c & pi**2/(2*Q),dimenl*Rb*t_bath*log(s)
2729 c-----------------------------------------------------------------
2730 double precision function HNose_nh(eki,e)
2731 implicit real*8 (a-h,o-z)
2732 include 'DIMENSIONS'
2734 HNose_nh=eki+e+dimen3*Rb*t_bath*xlogs(1)+qmass(1)*vlogs(1)**2/2
2736 HNose_nh=HNose_nh+qmass(i)*vlogs(i)**2/2+Rb*t_bath*xlogs(i)
2738 c write(4,'(5e15.5)')
2739 c & vlogs(1),xlogs(1),HNose,eki,e
2742 c-----------------------------------------------------------------
2743 SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8)
2744 implicit real*8 (a-h,o-z)
2745 include 'DIMENSIONS'
2747 double precision akin,gnkt,dt,aa,gkt,scale
2748 double precision wdti(maxyosh),wdti2(maxyosh),
2749 & wdti4(maxyosh),wdti8(maxyosh)
2750 integer i,iresn,iyosh,inos,nnos1
2759 C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE
2760 C INTEGRATION FROM t=0 TO t=DT/2
2761 C GET THE TOTAL KINETIC ENERGY
2763 c CALL GETKINP(MASS,VX,VY,VZ,AKIN)
2765 GLOGS(1) = (AKIN - GNKT)/QMASS(1)
2766 C START THE MULTIPLE TIME STEP PROCEDURE
2769 C UPDATE THE THERMOSTAT VELOCITIES
2770 VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2772 AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) )
2773 VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA
2774 & + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA
2776 C UPDATE THE PARTICLE VELOCITIES
2777 AA = EXP(-WDTI2(IYOSH)*VLOGS(1) )
2780 GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1)
2781 C UPDATE THE THERMOSTAT POSITIONS
2783 XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH)
2785 C UPDATE THE THERMOSTAT VELOCITIES
2787 AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) )
2788 VLOGS(INOS) = VLOGS(INOS)*AA*AA
2789 & + WDTI4(IYOSH)*GLOGS(INOS)*AA
2790 GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS)
2791 & -GKT)/QMASS(INOS+1)
2793 VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2796 C UPDATE THE PARTICLE VELOCITIES
2797 c outside of this subroutine
2799 c VX(I) = VX(I)*SCALE
2800 c VY(I) = VY(I)*SCALE
2801 c VZ(I) = VZ(I)*SCALE
2805 c-----------------------------------------------------------------
2806 subroutine tnp1_respa_i_step1
2807 c Applying Nose-Poincare algorithm - step 1 to coordinates
2808 c JPSJ 70 75 (2001) S. Nose
2810 c d_t is not updated here
2812 implicit real*8 (a-h,o-z)
2813 include 'DIMENSIONS'
2814 include 'COMMON.CONTROL'
2815 include 'COMMON.VAR'
2817 include 'COMMON.CHAIN'
2818 include 'COMMON.DERIV'
2819 include 'COMMON.GEO'
2820 include 'COMMON.LOCAL'
2821 include 'COMMON.INTERACT'
2822 include 'COMMON.IOUNITS'
2823 include 'COMMON.NAMES'
2824 double precision adt,adt2,tmp
2826 tmp=1+pi_np/(2*Q_np)*0.5*d_time
2829 s12_dt=d_time/s12_np
2830 d_time_s12=d_time*0.5*s12_np
2833 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2834 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2838 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2839 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2843 if (itype(i).ne.10) then
2846 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2847 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2853 c---------------------------------------------------------------------
2854 subroutine tnp1_respa_i_step2
2855 c Step 2 of the velocity Verlet algorithm: update velocities
2856 implicit real*8 (a-h,o-z)
2857 include 'DIMENSIONS'
2858 include 'COMMON.CONTROL'
2859 include 'COMMON.VAR'
2861 include 'COMMON.CHAIN'
2862 include 'COMMON.DERIV'
2863 include 'COMMON.GEO'
2864 include 'COMMON.LOCAL'
2865 include 'COMMON.INTERACT'
2866 include 'COMMON.IOUNITS'
2867 include 'COMMON.NAMES'
2869 double precision d_time_s12
2873 d_t(j,i)=d_t_new(j,i)
2880 d_time_s12=0.5d0*s12_np*d_time
2883 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
2887 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
2891 if (itype(i).ne.10) then
2894 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
2899 pistar=pistar+(EK-0.5*(E_old+potE)
2900 & -dimen3*Rb*t_bath*log(s12_np)+Csplit-dimen3*Rb*t_bath)*d_time
2901 tmp=1+pistar/(2*Q_np)*0.5*d_time
2907 c-------------------------------------------------------
2909 subroutine tnp1_step1
2910 c Applying Nose-Poincare algorithm - step 1 to coordinates
2911 c JPSJ 70 75 (2001) S. Nose
2913 c d_t is not updated here
2915 implicit real*8 (a-h,o-z)
2916 include 'DIMENSIONS'
2917 include 'COMMON.CONTROL'
2918 include 'COMMON.VAR'
2920 include 'COMMON.CHAIN'
2921 include 'COMMON.DERIV'
2922 include 'COMMON.GEO'
2923 include 'COMMON.LOCAL'
2924 include 'COMMON.INTERACT'
2925 include 'COMMON.IOUNITS'
2926 include 'COMMON.NAMES'
2927 double precision adt,adt2,tmp
2929 tmp=1+pi_np/(2*Q_np)*0.5*d_time
2932 s12_dt=d_time/s12_np
2933 d_time_s12=d_time*0.5*s12_np
2936 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2937 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2941 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2942 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2946 if (itype(i).ne.10) then
2949 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2950 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2956 c---------------------------------------------------------------------
2957 subroutine tnp1_step2
2958 c Step 2 of the velocity Verlet algorithm: update velocities
2959 implicit real*8 (a-h,o-z)
2960 include 'DIMENSIONS'
2961 include 'COMMON.CONTROL'
2962 include 'COMMON.VAR'
2964 include 'COMMON.CHAIN'
2965 include 'COMMON.DERIV'
2966 include 'COMMON.GEO'
2967 include 'COMMON.LOCAL'
2968 include 'COMMON.INTERACT'
2969 include 'COMMON.IOUNITS'
2970 include 'COMMON.NAMES'
2972 double precision d_time_s12
2976 d_t(j,i)=d_t_new(j,i)
2983 d_time_s12=0.5d0*s12_np*d_time
2986 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
2990 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
2994 if (itype(i).ne.10) then
2997 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
3002 cd write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np
3003 pistar=pistar+(EK-0.5*(E_old+potE)
3004 & -dimen3*Rb*t_bath*log(s12_np)+H0-dimen3*Rb*t_bath)*d_time
3005 tmp=1+pistar/(2*Q_np)*0.5*d_time
3012 c-----------------------------------------------------------------
3013 subroutine tnp_respa_i_step1
3014 c Applying Nose-Poincare algorithm - step 1 to coordinates
3015 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3017 c d_t is not updated here, it is destroyed
3019 implicit real*8 (a-h,o-z)
3020 include 'DIMENSIONS'
3021 include 'COMMON.CONTROL'
3022 include 'COMMON.VAR'
3024 include 'COMMON.CHAIN'
3025 include 'COMMON.DERIV'
3026 include 'COMMON.GEO'
3027 include 'COMMON.LOCAL'
3028 include 'COMMON.INTERACT'
3029 include 'COMMON.IOUNITS'
3030 include 'COMMON.NAMES'
3031 double precision C_np,d_time_s,tmp,d_time_ss
3033 d_time_s=d_time*0.5*s_np
3034 ct2 d_time_s=d_time*0.5*s12_np
3037 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3041 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3045 if (itype(i).ne.10) then
3048 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3055 d_t(j,i)=d_t_new(j,i)
3062 C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-Csplit)
3065 pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3066 tmp=0.5*d_time*pistar/Q_np
3067 s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3069 d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3070 ct2 d_time_ss=d_time/s12_np
3071 c d_time_ss=0.5*d_time*(1.0/sold_np+1.0/s_np)
3074 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3078 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3082 if (itype(i).ne.10) then
3085 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3092 c---------------------------------------------------------------------
3094 subroutine tnp_respa_i_step2
3095 c Step 2 of the velocity Verlet algorithm: update velocities
3096 implicit real*8 (a-h,o-z)
3097 include 'DIMENSIONS'
3098 include 'COMMON.CONTROL'
3099 include 'COMMON.VAR'
3101 include 'COMMON.CHAIN'
3102 include 'COMMON.DERIV'
3103 include 'COMMON.GEO'
3104 include 'COMMON.LOCAL'
3105 include 'COMMON.INTERACT'
3106 include 'COMMON.IOUNITS'
3107 include 'COMMON.NAMES'
3109 double precision d_time_s
3111 EK=EK*(s_np/s12_np)**2
3112 HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3113 pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath
3116 cr print '(a,5f)','i_step2',EK,potE,HNose1,pi_np,E_long
3117 d_time_s=d_time*0.5*s12_np
3118 c d_time_s=d_time*0.5*s_np
3121 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3125 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3129 if (itype(i).ne.10) then
3132 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
3141 c-----------------------------------------------------------------
3142 subroutine tnp_respa_step1
3143 c Applying Nose-Poincare algorithm - step 1 to vel for RESPA
3144 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3146 c d_t is not updated here, it is destroyed
3148 implicit real*8 (a-h,o-z)
3149 include 'DIMENSIONS'
3150 include 'COMMON.CONTROL'
3151 include 'COMMON.VAR'
3153 include 'COMMON.CHAIN'
3154 include 'COMMON.DERIV'
3155 include 'COMMON.GEO'
3156 include 'COMMON.LOCAL'
3157 include 'COMMON.INTERACT'
3158 include 'COMMON.IOUNITS'
3159 include 'COMMON.NAMES'
3160 double precision C_np,d_time_s,tmp,d_time_ss
3161 double precision energia(0:n_ene)
3163 d_time_s=d_time*0.5*s_np
3166 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3170 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3174 if (itype(i).ne.10) then
3177 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3183 c C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3186 c pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3187 c tmp=0.5*d_time*pistar/Q_np
3188 c s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3189 c write(iout,*) 'tnp_respa_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3195 c-------------------------------------
3196 c test of reviewer's comment
3197 pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3198 cr print '(a,3f)','1 pi_np,s_np',pi_np,s_np,E_long
3199 c-------------------------------------
3203 c---------------------------------------------------------------------
3204 subroutine tnp_respa_step2
3205 c Step 2 of the velocity Verlet algorithm: update velocities for RESPA
3206 implicit real*8 (a-h,o-z)
3207 include 'DIMENSIONS'
3208 include 'COMMON.CONTROL'
3209 include 'COMMON.VAR'
3211 include 'COMMON.CHAIN'
3212 include 'COMMON.DERIV'
3213 include 'COMMON.GEO'
3214 include 'COMMON.LOCAL'
3215 include 'COMMON.INTERACT'
3216 include 'COMMON.IOUNITS'
3217 include 'COMMON.NAMES'
3219 double precision d_time_s
3225 ct HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3226 ct pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3227 ct & -0.5*d_time*(HNose1-H0)
3229 c-------------------------------------
3230 c test of reviewer's comment
3231 pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3232 cr print '(a,3f)','2 pi_np,s_np',pi_np,s_np,E_long
3233 c-------------------------------------
3234 d_time_s=d_time*0.5*s_np
3237 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3241 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3245 if (itype(i).ne.10) then
3248 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3257 c---------------------------------------------------------------------
3258 subroutine tnp_step1
3259 c Applying Nose-Poincare algorithm - step 1 to coordinates
3260 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3262 c d_t is not updated here, it is destroyed
3264 implicit real*8 (a-h,o-z)
3265 include 'DIMENSIONS'
3266 include 'COMMON.CONTROL'
3267 include 'COMMON.VAR'
3269 include 'COMMON.CHAIN'
3270 include 'COMMON.DERIV'
3271 include 'COMMON.GEO'
3272 include 'COMMON.LOCAL'
3273 include 'COMMON.INTERACT'
3274 include 'COMMON.IOUNITS'
3275 include 'COMMON.NAMES'
3276 double precision C_np,d_time_s,tmp,d_time_ss
3278 d_time_s=d_time*0.5*s_np
3281 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3285 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3289 if (itype(i).ne.10) then
3292 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3299 d_t(j,i)=d_t_new(j,i)
3306 C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3309 pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3310 tmp=0.5*d_time*pistar/Q_np
3311 s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3312 c write(iout,*) 'tnp_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3314 d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3317 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3321 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3325 if (itype(i).ne.10) then
3328 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3335 c-----------------------------------------------------------------
3336 subroutine tnp_step2
3337 c Step 2 of the velocity Verlet algorithm: update velocities
3338 implicit real*8 (a-h,o-z)
3339 include 'DIMENSIONS'
3340 include 'COMMON.CONTROL'
3341 include 'COMMON.VAR'
3343 include 'COMMON.CHAIN'
3344 include 'COMMON.DERIV'
3345 include 'COMMON.GEO'
3346 include 'COMMON.LOCAL'
3347 include 'COMMON.INTERACT'
3348 include 'COMMON.IOUNITS'
3349 include 'COMMON.NAMES'
3351 double precision d_time_s
3353 EK=EK*(s_np/s12_np)**2
3354 HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3355 pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3356 & -0.5*d_time*(HNose1-H0)
3358 cd write(iout,'(a,4f)') 'mmm',EK,potE,HNose1,pi_np
3359 d_time_s=d_time*0.5*s12_np
3362 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3366 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3370 if (itype(i).ne.10) then
3373 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s