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 cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
600 HNose1=Hnose_nh(EK,potE)
603 cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
619 if (itype(i).ne.10) then
623 vtnp(itnp)=d_t(j,inres)
628 c Transform velocities from UNRES coordinate space to cartesian and Gvec
635 vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
636 vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
638 vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
642 write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
645 write (iout,*) "Velocities, step 2"
647 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
648 & (d_t(j,i+nres),j=1,3)
655 c-------------------------------------------------------------------------------
656 subroutine RESPA_step(itime)
657 c-------------------------------------------------------------------------------
658 c Perform a single RESPA step.
659 c-------------------------------------------------------------------------------
660 implicit real*8 (a-h,o-z)
664 integer IERROR,ERRCODE
666 include 'COMMON.SETUP'
667 include 'COMMON.CONTROL'
671 include 'COMMON.LANGEVIN'
673 include 'COMMON.LANGEVIN.lang0'
675 include 'COMMON.CHAIN'
676 include 'COMMON.DERIV'
678 include 'COMMON.LOCAL'
679 include 'COMMON.INTERACT'
680 include 'COMMON.IOUNITS'
681 include 'COMMON.NAMES'
682 include 'COMMON.TIME1'
683 double precision energia_short(0:n_ene),
684 & energia_long(0:n_ene)
685 double precision cm(3),L(3),vcm(3),incr(3)
686 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
687 & d_a_old0(3,0:maxres2)
688 integer ilen,count,rstcount
691 integer maxcount_scale /10/
692 common /gucio/ cm,energia_short
693 double precision stochforcvec(MAXRES6)
694 common /stochcalc/ stochforcvec
697 double precision vtnp(maxres6), vtnp_(maxres6), vtnp_a(maxres6)
698 common /cipiszcze/ itt
701 if (large.and. mod(itime,ntwe).eq.0) then
702 write (iout,*) "***************** RESPA itime",itime
703 write (iout,*) "Cartesian and internal coordinates: step 0"
705 call pdbout(0.0d0,"cipiszcze",iout)
707 write (iout,*) "Accelerations from long-range forces"
709 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
710 & (d_a(j,i+nres),j=1,3)
712 write (iout,*) "Velocities, step 0"
714 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
715 & (d_t(j,i+nres),j=1,3)
720 c Perform the initial RESPA step (increment velocities)
721 c write (iout,*) "*********************** RESPA ini"
727 if (tnh.and..not.xiresp) then
728 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
731 d_t(j,i)=d_t(j,i)*scale_nh
738 cd if(tnp .or. tnp1) then
739 cd write (iout,'(a,3f)') "EE1 NP S, pi",totT, s_np, pi_np
743 if (mod(itime,ntwe).eq.0 .and. large) then
744 write (iout,*) "Velocities, end"
746 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
747 & (d_t(j,i+nres),j=1,3)
751 c Compute the short-range forces
757 C 7/2/2009 commented out
759 c call etotal_short(energia_short)
760 if (tnp.or.tnp1) potE=energia_short(0)
763 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
766 d_a(j,i)=d_a_short(j,i)
770 if (large.and. mod(itime,ntwe).eq.0) then
771 write (iout,*) "energia_short",energia_short(0)
772 write (iout,*) "Accelerations from short-range forces"
774 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
775 & (d_a(j,i+nres),j=1,3)
780 t_enegrad=t_enegrad+MPI_Wtime()-tt0
782 t_enegrad=t_enegrad+tcpu()-tt0
787 if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
788 d_a_old(j,i)=d_a(j,i)
791 c 6/30/08 A-MTS: attempt at increasing the split number
794 dc_old0(j,i)=dc_old(j,i)
795 d_t_old0(j,i)=d_t_old(j,i)
796 d_a_old0(j,i)=d_a_old(j,i)
799 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
800 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
807 c write (iout,*) "itime",itime," ntime_split",ntime_split
808 c Split the time step
809 d_time=d_time0/ntime_split
810 c Perform the short-range RESPA steps (velocity Verlet increments of
811 c positions and velocities using short-range forces)
812 c write (iout,*) "*********************** RESPA split"
813 do itsplit=1,ntime_split
816 else if (lang.eq.2 .or. lang.eq.3) then
818 call stochastic_force(stochforcvec)
821 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
823 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
828 c First step of the velocity Verlet algorithm
833 else if (lang.eq.3) then
835 call sd_verlet1_ciccotti
837 else if (lang.eq.1) then
840 call tnp1_respa_i_step1
842 call tnp_respa_i_step1
844 if (tnh.and.xiresp) then
846 call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
849 d_t_old(j,i)=d_t_old(j,i)*scale_nh
852 cd write(iout,*) "SSS1",itsplit,EK,scale_nh
856 c Build the chain from the newly calculated coordinates
858 if (rattle) call rattle1
860 if (large.and. mod(itime,ntwe).eq.0) then
861 write (iout,*) "***** ITSPLIT",itsplit
862 write (iout,*) "Cartesian and internal coordinates: step 1"
863 call pdbout(0.0d0,"cipiszcze",iout)
866 write (iout,*) "Velocities, step 1"
868 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
869 & (d_t(j,i+nres),j=1,3)
878 c Calculate energy and forces
880 call etotal_short(energia_short)
882 potE=energia_short(0)
885 t_eshort=t_eshort+MPI_Wtime()-tt0
887 t_eshort=t_eshort+tcpu()-tt0
891 c Get the new accelerations
893 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
896 d_a_short(j,i)=d_a(j,i)
900 if (large.and. mod(itime,ntwe).eq.0) then
901 write (iout,*)"energia_short",energia_short(0)
902 write (iout,*) "Accelerations from short-range forces"
904 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
905 & (d_a(j,i+nres),j=1,3)
910 c Determine maximum acceleration and scale down the timestep if needed
912 amax=amax/ntime_split**2
913 call predict_edrift(epdrift)
914 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
915 & write (iout,*) "amax",amax," damax",damax,
916 & " epdrift",epdrift," epdriftmax",epdriftmax
917 c Exit loop and try with increased split number if the change of
918 c acceleration is too big
919 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
920 if (ntime_split.lt.maxtime_split) then
922 ntime_split=ntime_split*2
925 dc_old(j,i)=dc_old0(j,i)
926 d_t_old(j,i)=d_t_old0(j,i)
927 d_a_old(j,i)=d_a_old0(j,i)
930 write (iout,*) "acceleration/energy drift too large",amax,
931 & epdrift," split increased to ",ntime_split," itime",itime,
936 & "Uh-hu. Bumpy landscape. Maximum splitting number",
938 & " already reached!!! Trying to carry on!"
942 t_enegrad=t_enegrad+MPI_Wtime()-tt0
944 t_enegrad=t_enegrad+tcpu()-tt0
946 c Second step of the velocity Verlet algorithm
951 else if (lang.eq.3) then
953 call sd_verlet2_ciccotti
955 else if (lang.eq.1) then
958 call tnp1_respa_i_step2
960 call tnp_respa_i_step2
963 if (tnh.and.xiresp) then
965 call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
968 d_t(j,i)=d_t(j,i)*scale_nh
971 cd write(iout,*) "SSS2",itsplit,EK,scale_nh
974 if (rattle) call rattle2
975 c Backup the coordinates, velocities, and accelerations
976 if (tnp .or. tnp1) then
979 d_t_old(j,i)=d_t(j,i)
980 if (tnp) d_t(j,i)=d_t(j,i)/s_np
981 if (tnp1) d_t(j,i)=d_t(j,i)/s_np
989 if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
990 d_a_old(j,i)=d_a(j,i)
997 c Restore the time step
999 c Compute long-range forces
1006 call etotal_long(energia_long)
1007 E_long=energia_long(0)
1008 potE=energia_short(0)+energia_long(0)
1011 t_elong=t_elong+MPI_Wtime()-tt0
1013 t_elong=t_elong+tcpu()-tt0
1019 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1021 t_enegrad=t_enegrad+tcpu()-tt0
1023 c Compute accelerations from long-range forces
1025 if (large.and. mod(itime,ntwe).eq.0) then
1026 write (iout,*) "energia_long",energia_long(0)
1027 write (iout,*) "Cartesian and internal coordinates: step 2"
1029 call pdbout(0.0d0,"cipiszcze",iout)
1031 write (iout,*) "Accelerations from long-range forces"
1033 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1034 & (d_a(j,i+nres),j=1,3)
1036 write (iout,*) "Velocities, step 2"
1038 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1039 & (d_t(j,i+nres),j=1,3)
1043 c Compute the final RESPA step (increment velocities)
1044 c write (iout,*) "*********************** RESPA fin"
1046 call tnp_respa_step2
1048 call tnp_respa_step2
1051 if (tnh.and..not.xiresp) then
1053 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
1056 d_t(j,i)=d_t(j,i)*scale_nh
1062 if (tnp .or. tnp1) then
1065 d_t(j,i)=d_t_old(j,i)/s_np
1070 c Compute the complete potential energy
1072 potEcomp(i)=energia_short(i)+energia_long(i)
1074 potE=potEcomp(0)-potEcomp(20)
1075 c potE=energia_short(0)+energia_long(0)
1077 c Calculate the kinetic and the total energy and the kinetic temperature
1080 c Couple the system to Berendsen bath if needed
1081 if (tbf .and. lang.eq.0) then
1084 kinetic_T=2.0d0/(dimen3*Rb)*EK
1085 c Backup the coordinates, velocities, and accelerations
1087 if (mod(itime,ntwe).eq.0 .and. large) then
1088 write (iout,*) "Velocities, end"
1090 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1091 & (d_t(j,i+nres),j=1,3)
1095 if (mod(itime,ntwe).eq.0) then
1097 if(tnp .or. tnp1) then
1099 write (iout,'(a3,7f20.10)') "TTT",EK,s_np,potE,pi_np,Csplit,
1100 & E_long,energia_short(0)
1102 write (iout,'(a3,7f20.10)') "TTT",EK,s_np,potE,pi_np,Csplit,
1103 & E_long,energia_short(0)
1105 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
1107 cd write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
1108 cd & ,EK+potE+pi_np**2/(2*Q_np)+dimen3*0.001986d0*t_bath*log(s_np)
1109 cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
1111 cd write (iout,'(a,3f)') "EE2 NP S, pi",totT, s_np, pi_np
1115 HNose1=Hnose_nh(EK,potE)
1117 cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
1135 if (itype(i).ne.10) then
1139 vtnp(itnp)=d_t(j,inres)
1144 c Transform velocities from UNRES coordinate space to cartesian and Gvec
1151 vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
1152 vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
1154 vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
1158 write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
1168 c---------------------------------------------------------------------
1169 subroutine RESPA_vel
1170 c First and last RESPA step (incrementing velocities using long-range
1172 implicit real*8 (a-h,o-z)
1173 include 'DIMENSIONS'
1174 include 'COMMON.CONTROL'
1175 include 'COMMON.VAR'
1177 include 'COMMON.CHAIN'
1178 include 'COMMON.DERIV'
1179 include 'COMMON.GEO'
1180 include 'COMMON.LOCAL'
1181 include 'COMMON.INTERACT'
1182 include 'COMMON.IOUNITS'
1183 include 'COMMON.NAMES'
1185 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1189 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1193 if (itype(i).ne.10) then
1196 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1202 c-----------------------------------------------------------------
1204 c Applying velocity Verlet algorithm - step 1 to coordinates
1205 implicit real*8 (a-h,o-z)
1206 include 'DIMENSIONS'
1207 include 'COMMON.CONTROL'
1208 include 'COMMON.VAR'
1210 include 'COMMON.CHAIN'
1211 include 'COMMON.DERIV'
1212 include 'COMMON.GEO'
1213 include 'COMMON.LOCAL'
1214 include 'COMMON.INTERACT'
1215 include 'COMMON.IOUNITS'
1216 include 'COMMON.NAMES'
1217 double precision adt,adt2
1220 write (iout,*) "VELVERLET1 START: DC"
1222 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1223 & (dc(j,i+nres),j=1,3)
1227 adt=d_a_old(j,0)*d_time
1229 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1230 d_t_new(j,0)=d_t_old(j,0)+adt2
1231 d_t(j,0)=d_t_old(j,0)+adt
1235 adt=d_a_old(j,i)*d_time
1237 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1238 d_t_new(j,i)=d_t_old(j,i)+adt2
1239 d_t(j,i)=d_t_old(j,i)+adt
1243 if (itype(i).ne.10) then
1246 adt=d_a_old(j,inres)*d_time
1248 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1249 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1250 d_t(j,inres)=d_t_old(j,inres)+adt
1255 write (iout,*) "VELVERLET1 END: DC"
1257 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1258 & (dc(j,i+nres),j=1,3)
1263 c---------------------------------------------------------------------
1265 c Step 2 of the velocity Verlet algorithm: update velocities
1266 implicit real*8 (a-h,o-z)
1267 include 'DIMENSIONS'
1268 include 'COMMON.CONTROL'
1269 include 'COMMON.VAR'
1271 include 'COMMON.CHAIN'
1272 include 'COMMON.DERIV'
1273 include 'COMMON.GEO'
1274 include 'COMMON.LOCAL'
1275 include 'COMMON.INTERACT'
1276 include 'COMMON.IOUNITS'
1277 include 'COMMON.NAMES'
1279 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1283 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1287 if (itype(i).ne.10) then
1290 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1296 c-----------------------------------------------------------------
1297 subroutine sddir_precalc
1298 c Applying velocity Verlet algorithm - step 1 to coordinates
1299 implicit real*8 (a-h,o-z)
1300 include 'DIMENSIONS'
1304 include 'COMMON.CONTROL'
1305 include 'COMMON.VAR'
1308 include 'COMMON.LANGEVIN'
1310 include 'COMMON.LANGEVIN.lang0'
1312 include 'COMMON.CHAIN'
1313 include 'COMMON.DERIV'
1314 include 'COMMON.GEO'
1315 include 'COMMON.LOCAL'
1316 include 'COMMON.INTERACT'
1317 include 'COMMON.IOUNITS'
1318 include 'COMMON.NAMES'
1319 include 'COMMON.TIME1'
1320 double precision stochforcvec(MAXRES6)
1321 common /stochcalc/ stochforcvec
1323 c Compute friction and stochastic forces
1332 time_fric=time_fric+MPI_Wtime()-time00
1335 time_fric=time_fric+tcpu()-time00
1338 call stochastic_force(stochforcvec)
1340 time_stoch=time_stoch+MPI_Wtime()-time00
1342 time_stoch=time_stoch+tcpu()-time00
1345 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1346 c forces (d_as_work)
1348 call ginv_mult(fric_work, d_af_work)
1349 call ginv_mult(stochforcvec, d_as_work)
1352 c---------------------------------------------------------------------
1353 subroutine sddir_verlet1
1354 c Applying velocity Verlet algorithm - step 1 to velocities
1355 implicit real*8 (a-h,o-z)
1356 include 'DIMENSIONS'
1357 include 'COMMON.CONTROL'
1358 include 'COMMON.VAR'
1361 include 'COMMON.LANGEVIN'
1363 include 'COMMON.LANGEVIN.lang0'
1365 include 'COMMON.CHAIN'
1366 include 'COMMON.DERIV'
1367 include 'COMMON.GEO'
1368 include 'COMMON.LOCAL'
1369 include 'COMMON.INTERACT'
1370 include 'COMMON.IOUNITS'
1371 include 'COMMON.NAMES'
1372 c Revised 3/31/05 AL: correlation between random contributions to
1373 c position and velocity increments included.
1374 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1375 double precision adt,adt2
1377 c Add the contribution from BOTH friction and stochastic force to the
1378 c coordinates, but ONLY the contribution from the friction forces to velocities
1381 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1382 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1383 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1384 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1385 d_t(j,0)=d_t_old(j,0)+adt
1390 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1391 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1392 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1393 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1394 d_t(j,i)=d_t_old(j,i)+adt
1399 if (itype(i).ne.10) then
1402 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1403 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1404 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1405 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1406 d_t(j,inres)=d_t_old(j,inres)+adt
1413 c---------------------------------------------------------------------
1414 subroutine sddir_verlet2
1415 c Calculating the adjusted velocities for accelerations
1416 implicit real*8 (a-h,o-z)
1417 include 'DIMENSIONS'
1418 include 'COMMON.CONTROL'
1419 include 'COMMON.VAR'
1422 include 'COMMON.LANGEVIN'
1424 include 'COMMON.LANGEVIN.lang0'
1426 include 'COMMON.CHAIN'
1427 include 'COMMON.DERIV'
1428 include 'COMMON.GEO'
1429 include 'COMMON.LOCAL'
1430 include 'COMMON.INTERACT'
1431 include 'COMMON.IOUNITS'
1432 include 'COMMON.NAMES'
1433 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1434 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1435 c Revised 3/31/05 AL: correlation between random contributions to
1436 c position and velocity increments included.
1437 c The correlation coefficients are calculated at low-friction limit.
1438 c Also, friction forces are now not calculated with new velocities.
1440 c call friction_force
1441 call stochastic_force(stochforcvec)
1443 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1444 c forces (d_as_work)
1446 call ginv_mult(stochforcvec, d_as_work1)
1452 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1453 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1458 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1459 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1464 if (itype(i).ne.10) then
1467 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1468 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1469 & +cos60*d_as_work1(ind+j))*d_time
1476 c---------------------------------------------------------------------
1477 subroutine max_accel
1479 c Find the maximum difference in the accelerations of the the sites
1480 c at the beginning and the end of the time step.
1482 implicit real*8 (a-h,o-z)
1483 include 'DIMENSIONS'
1484 include 'COMMON.CONTROL'
1485 include 'COMMON.VAR'
1487 include 'COMMON.CHAIN'
1488 include 'COMMON.DERIV'
1489 include 'COMMON.GEO'
1490 include 'COMMON.LOCAL'
1491 include 'COMMON.INTERACT'
1492 include 'COMMON.IOUNITS'
1493 double precision aux(3),accel(3),accel_old(3),dacc
1495 c aux(j)=d_a(j,0)-d_a_old(j,0)
1496 accel_old(j)=d_a_old(j,0)
1503 c 7/3/08 changed to asymmetric difference
1505 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1506 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1507 accel(j)=accel(j)+0.5d0*d_a(j,i)
1508 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1509 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1510 dacc=dabs(accel(j)-accel_old(j))
1511 if (dacc.gt.amax) amax=dacc
1519 accel_old(j)=d_a_old(j,0)
1524 accel_old(j)=accel_old(j)+d_a_old(j,1)
1525 accel(j)=accel(j)+d_a(j,1)
1529 if (itype(i).ne.10) then
1531 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1532 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1533 accel(j)=accel(j)+d_a(j,i+nres)
1537 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1538 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1539 dacc=dabs(accel(j)-accel_old(j))
1540 if (dacc.gt.amax) amax=dacc
1544 accel_old(j)=accel_old(j)+d_a_old(j,i)
1545 accel(j)=accel(j)+d_a(j,i)
1546 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1551 c---------------------------------------------------------------------
1552 subroutine predict_edrift(epdrift)
1554 c Predict the drift of the potential energy
1556 implicit real*8 (a-h,o-z)
1557 include 'DIMENSIONS'
1558 include 'COMMON.CONTROL'
1559 include 'COMMON.VAR'
1561 include 'COMMON.CHAIN'
1562 include 'COMMON.DERIV'
1563 include 'COMMON.GEO'
1564 include 'COMMON.LOCAL'
1565 include 'COMMON.INTERACT'
1566 include 'COMMON.IOUNITS'
1567 include 'COMMON.MUCA'
1568 double precision epdrift,epdriftij
1569 c Drift of the potential energy
1575 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1576 if (lmuca) epdriftij=epdriftij*factor
1577 c write (iout,*) "back",i,j,epdriftij
1578 if (epdriftij.gt.epdrift) epdrift=epdriftij
1582 if (itype(i).ne.10) then
1585 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1586 if (lmuca) epdriftij=epdriftij*factor
1587 c write (iout,*) "side",i,j,epdriftij
1588 if (epdriftij.gt.epdrift) epdrift=epdriftij
1592 epdrift=0.5d0*epdrift*d_time*d_time
1593 c write (iout,*) "epdrift",epdrift
1596 c-----------------------------------------------------------------------
1597 subroutine verlet_bath
1599 c Coupling to the thermostat by using the Berendsen algorithm
1601 implicit real*8 (a-h,o-z)
1602 include 'DIMENSIONS'
1603 include 'COMMON.CONTROL'
1604 include 'COMMON.VAR'
1606 include 'COMMON.CHAIN'
1607 include 'COMMON.DERIV'
1608 include 'COMMON.GEO'
1609 include 'COMMON.LOCAL'
1610 include 'COMMON.INTERACT'
1611 include 'COMMON.IOUNITS'
1612 include 'COMMON.NAMES'
1613 double precision T_half,fact
1615 T_half=2.0d0/(dimen3*Rb)*EK
1616 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1617 c write(iout,*) "T_half", T_half
1618 c write(iout,*) "EK", EK
1619 c write(iout,*) "fact", fact
1621 d_t(j,0)=fact*d_t(j,0)
1625 d_t(j,i)=fact*d_t(j,i)
1629 if (itype(i).ne.10) then
1632 d_t(j,inres)=fact*d_t(j,inres)
1638 c---------------------------------------------------------
1640 c Set up the initial conditions of a MD simulation
1641 implicit real*8 (a-h,o-z)
1642 include 'DIMENSIONS'
1646 integer IERROR,ERRCODE
1648 include 'COMMON.SETUP'
1649 include 'COMMON.CONTROL'
1650 include 'COMMON.VAR'
1653 include 'COMMON.LANGEVIN'
1655 include 'COMMON.LANGEVIN.lang0'
1657 include 'COMMON.CHAIN'
1658 include 'COMMON.DERIV'
1659 include 'COMMON.GEO'
1660 include 'COMMON.LOCAL'
1661 include 'COMMON.INTERACT'
1662 include 'COMMON.IOUNITS'
1663 include 'COMMON.NAMES'
1664 include 'COMMON.REMD'
1665 real*8 energia_long(0:n_ene),
1666 & energia_short(0:n_ene),vcm(3),incr(3),E_short
1667 double precision cm(3),L(3),xv,sigv,lowb,highb
1668 double precision varia(maxvar)
1676 c write(iout,*) "d_time", d_time
1677 c Compute the standard deviations of stochastic forces for Langevin dynamics
1678 c if the friction coefficients do not depend on surface area
1679 if (lang.gt.0 .and. .not.surfarea) then
1681 stdforcp(i)=stdfp*dsqrt(gamp)
1684 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
1687 c Open the pdb file for snapshotshots
1690 if (ilen(tmpdir).gt.0)
1691 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1692 & liczba(:ilen(liczba))//".pdb")
1694 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1698 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1699 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1700 & liczba(:ilen(liczba))//".x")
1701 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1704 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1705 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1706 & liczba(:ilen(liczba))//".cx")
1707 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1713 if (ilen(tmpdir).gt.0)
1714 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1715 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1717 if (ilen(tmpdir).gt.0)
1718 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1719 cartname=prefix(:ilen(prefix))//"_MD.cx"
1723 write (qstr,'(256(1h ))')
1726 iq = qinfrag(i,iset)*10
1727 iw = wfrag(i,iset)/100
1729 if(me.eq.king.or..not.out1file)
1730 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1731 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1736 iq = qinpair(i,iset)*10
1737 iw = wpair(i,iset)/100
1739 if(me.eq.king.or..not.out1file)
1740 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1741 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1745 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1747 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1749 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1751 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1755 if (restart1file) then
1757 & inquire(file=mremd_rst_name,exist=file_exist)
1758 write (*,*) me," Before broadcast: file_exist",file_exist
1760 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1762 write (*,*) me," After broadcast: file_exist",file_exist
1764 c inquire(file=mremd_rst_name,exist=file_exist)
1765 if(me.eq.king.or..not.out1file)
1766 & write(iout,*) "Initial state read by master and distributed"
1768 if (ilen(tmpdir).gt.0)
1769 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1770 & //liczba(:ilen(liczba))//'.rst')
1771 inquire(file=rest2name,exist=file_exist)
1774 if(.not.restart1file) then
1775 if(me.eq.king.or..not.out1file)
1776 & write(iout,*) "Initial state will be read from file ",
1777 & rest2name(:ilen(rest2name))
1780 call rescale_weights(t_bath)
1782 if(me.eq.king.or..not.out1file)then
1783 if (restart1file) then
1784 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1787 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1790 write(iout,*) "Initial velocities randomly generated"
1796 c Generate initial velocities
1797 if(me.eq.king.or..not.out1file)
1798 & write(iout,*) "Initial velocities randomly generated"
1802 c rest2name = prefix(:ilen(prefix))//'.rst'
1803 if(me.eq.king.or..not.out1file)then
1804 write (iout,*) "Initial velocities"
1806 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1807 & (d_t(j,i+nres),j=1,3)
1810 c Zeroing the total angular momentum of the system
1811 write(iout,*) "Calling the zero-angular
1812 & momentum subroutine"
1815 c Getting the potential energy and forces and velocities and accelerations
1817 c write (iout,*) "velocity of the center of the mass:"
1818 c write (iout,*) (vcm(j),j=1,3)
1820 d_t(j,0)=d_t(j,0)-vcm(j)
1822 c Removing the velocity of the center of mass
1824 if(me.eq.king.or..not.out1file)then
1825 write (iout,*) "vcm right after adjustment:"
1826 write (iout,*) (vcm(j),j=1,3)
1831 write (iout,*) "PREMINIM ",preminim
1832 if(iranconf.ne.0 .or. preminim) then
1834 print *, 'Calling OVERLAP_SC'
1835 call overlap_sc(fail)
1839 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1840 print *,'SC_move',nft_sc,etot
1841 if(me.eq.king.or..not.out1file)
1842 & write(iout,*) 'SC_move',nft_sc,etot
1846 print *, 'Calling MINIM_DC'
1847 call minim_dc(etot,iretcode,nfun)
1849 call geom_to_var(nvar,varia)
1850 print *,'Calling MINIMIZE.'
1851 call minimize(etot,varia,iretcode,nfun)
1852 call var_to_geom(nvar,varia)
1854 if(me.eq.king.or..not.out1file) then
1855 write(iout,*) "Minimized energy is",etot
1856 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1857 call etotal(potEcomp)
1858 call enerprint(potEcomp)
1862 call chainbuild_cart
1867 kinetic_T=2.0d0/(dimen3*Rb)*EK
1868 if(me.eq.king.or..not.out1file)then
1878 call etotal(potEcomp)
1881 t_etotal=t_etotal+MPI_Wtime()-tt0
1883 t_etotal=t_etotal+tcpu()-tt0
1888 if(tnp .or. tnp1) then
1891 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
1893 write(iout,*) 'H0= ',H0
1897 HNose1=Hnose_nh(EK,potE)
1899 write (iout,*) 'H0= ',H0
1905 if(me.eq.king.or..not.out1file)
1906 & write(iout,*) 'HMC',hmc_etot,potE,EK
1917 if (amax*d_time .gt. dvmax) then
1918 d_time=d_time*dvmax/amax
1919 if(me.eq.king.or..not.out1file) write (iout,*)
1920 & "Time step reduced to",d_time,
1921 & " because of too large initial acceleration."
1923 if(me.eq.king.or..not.out1file)then
1924 write(iout,*) "Potential energy and its components"
1925 call enerprint(potEcomp)
1926 c write(iout,*) (potEcomp(i),i=0,n_ene)
1928 potE=potEcomp(0)-potEcomp(20)
1931 if (ntwe.ne.0) call statout(itime)
1932 if(me.eq.king.or..not.out1file)
1933 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1934 & " Kinetic energy",EK," potential energy",potE,
1935 & " total energy",totE," maximum acceleration ",
1938 write (iout,*) "Initial coordinates"
1940 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1941 & (c(j,i+nres),j=1,3)
1943 write (iout,*) "Initial dC"
1945 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1946 & (dc(j,i+nres),j=1,3)
1948 write (iout,*) "Initial velocities"
1949 write (iout,"(13x,' backbone ',23x,' side chain')")
1951 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1952 & (d_t(j,i+nres),j=1,3)
1954 write (iout,*) "Initial accelerations"
1956 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1957 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1958 & (d_a(j,i+nres),j=1,3)
1964 d_t_old(j,i)=d_t(j,i)
1965 d_a_old(j,i)=d_a(j,i)
1967 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1976 call etotal_short(energia_short)
1979 t_eshort=t_eshort+MPI_Wtime()-tt0
1981 t_eshort=t_eshort+tcpu()-tt0
1985 if(tnp .or. tnp1) then
1986 E_short=energia_short(0)
1987 HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen3)
1990 c_new_var_csplit Csplit=H0-E_long
1991 c Csplit = H0-energia_short(0)
1992 write(iout,*) 'Csplit= ',Csplit
1998 if(.not.out1file .and. large) then
1999 write (iout,*) "energia_long",energia_long(0),
2000 & " energia_short",energia_short(0),
2001 & " total",energia_long(0)+energia_short(0)
2002 write (iout,*) "Initial fast-force accelerations"
2004 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2005 & (d_a(j,i+nres),j=1,3)
2008 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2011 d_a_short(j,i)=d_a(j,i)
2020 call etotal_long(energia_long)
2023 t_elong=t_elong+MPI_Wtime()-tt0
2025 t_elong=t_elong+tcpu()-tt0
2030 if(.not.out1file .and. large) then
2031 write (iout,*) "energia_long",energia_long(0)
2032 write (iout,*) "Initial slow-force accelerations"
2034 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2035 & (d_a(j,i+nres),j=1,3)
2039 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2041 t_enegrad=t_enegrad+tcpu()-tt0
2049 c-----------------------------------------------------------
2050 subroutine random_vel
2051 implicit real*8 (a-h,o-z)
2052 include 'DIMENSIONS'
2053 include 'COMMON.CONTROL'
2054 include 'COMMON.VAR'
2057 include 'COMMON.LANGEVIN'
2059 include 'COMMON.LANGEVIN.lang0'
2061 include 'COMMON.CHAIN'
2062 include 'COMMON.DERIV'
2063 include 'COMMON.GEO'
2064 include 'COMMON.LOCAL'
2065 include 'COMMON.INTERACT'
2066 include 'COMMON.IOUNITS'
2067 include 'COMMON.NAMES'
2068 include 'COMMON.TIME1'
2069 double precision xv,sigv,lowb,highb
2070 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2071 c First generate velocities in the eigenspace of the G matrix
2072 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2074 c write (iout,*) "RANDOM_VEL dimen",dimen
2080 sigv=dsqrt((Rb*t_bath)/geigen(i))
2083 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2084 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2085 c & " d_t_work_new",d_t_work_new(ii)
2095 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2098 c write (iout,*) "Ek from eigenvectors",Ek1
2100 c Transform velocities to UNRES coordinate space
2106 d_t_work(ind)=d_t_work(ind)
2107 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2109 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2113 c Transfer to the d_t vector
2115 d_t(j,0)=d_t_work(j)
2121 d_t(j,i)=d_t_work(ind)
2125 c write (iout,*) "d_t",i,(d_t(j,i),j=1,3)
2129 if (itype(i).ne.10) then
2132 d_t(j,i+nres)=d_t_work(ind)
2137 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2138 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2143 c-----------------------------------------------------------
2144 subroutine sd_verlet_p_setup
2145 c Sets up the parameters of stochastic Verlet algorithm
2146 implicit real*8 (a-h,o-z)
2147 include 'DIMENSIONS'
2151 include 'COMMON.CONTROL'
2152 include 'COMMON.VAR'
2155 include 'COMMON.LANGEVIN'
2157 include 'COMMON.LANGEVIN.lang0'
2159 include 'COMMON.CHAIN'
2160 include 'COMMON.DERIV'
2161 include 'COMMON.GEO'
2162 include 'COMMON.LOCAL'
2163 include 'COMMON.INTERACT'
2164 include 'COMMON.IOUNITS'
2165 include 'COMMON.NAMES'
2166 include 'COMMON.TIME1'
2167 double precision emgdt(MAXRES6),
2168 & pterm,vterm,rho,rhoc,vsig,
2169 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2170 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2171 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2172 logical lprn /.false./
2173 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2174 double precision ktm
2181 c AL 8/17/04 Code adapted from tinker
2183 c Get the frictional and random terms for stochastic dynamics in the
2184 c eigenspace of mass-scaled UNRES friction matrix
2187 gdt = fricgam(i) * d_time
2189 c Stochastic dynamics reduces to simple MD for zero friction
2191 if (gdt .le. zero) then
2192 pfric_vec(i) = 1.0d0
2193 vfric_vec(i) = d_time
2194 afric_vec(i) = 0.5d0 * d_time * d_time
2195 prand_vec(i) = 0.0d0
2196 vrand_vec1(i) = 0.0d0
2197 vrand_vec2(i) = 0.0d0
2199 c Analytical expressions when friction coefficient is large
2202 if (gdt .ge. gdt_radius) then
2205 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2206 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2207 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2208 vterm = 1.0d0 - egdt**2
2209 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2211 c Use series expansions when friction coefficient is small
2222 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2223 & - gdt5/120.0d0 + gdt6/720.0d0
2224 & - gdt7/5040.0d0 + gdt8/40320.0d0
2225 & - gdt9/362880.0d0) / fricgam(i)**2
2226 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2227 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2228 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2229 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2230 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2231 & + 127.0d0*gdt9/90720.0d0
2232 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2233 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2234 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2235 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2236 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2237 & - 17.0d0*gdt2/1280.0d0
2238 & + 17.0d0*gdt3/6144.0d0
2239 & + 40967.0d0*gdt4/34406400.0d0
2240 & - 57203.0d0*gdt5/275251200.0d0
2241 & - 1429487.0d0*gdt6/13212057600.0d0)
2244 c Compute the scaling factors of random terms for the nonzero friction case
2246 ktm = 0.5d0*d_time/fricgam(i)
2247 psig = dsqrt(ktm*pterm) / fricgam(i)
2248 vsig = dsqrt(ktm*vterm)
2249 rhoc = dsqrt(1.0d0 - rho*rho)
2251 vrand_vec1(i) = vsig * rho
2252 vrand_vec2(i) = vsig * rhoc
2257 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2260 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2261 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2265 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2268 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2269 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2270 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2271 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2272 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2273 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2276 t_sdsetup=t_sdsetup+MPI_Wtime()
2278 t_sdsetup=t_sdsetup+tcpu()-tt0
2282 c-------------------------------------------------------------
2283 subroutine eigtransf1(n,ndim,ab,d,c)
2286 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2292 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2298 c-------------------------------------------------------------
2299 subroutine eigtransf(n,ndim,a,b,d,c)
2302 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2308 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2314 c-------------------------------------------------------------
2315 subroutine sd_verlet1
2316 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2317 implicit real*8 (a-h,o-z)
2318 include 'DIMENSIONS'
2319 include 'COMMON.CONTROL'
2320 include 'COMMON.VAR'
2323 include 'COMMON.LANGEVIN'
2325 include 'COMMON.LANGEVIN.lang0'
2327 include 'COMMON.CHAIN'
2328 include 'COMMON.DERIV'
2329 include 'COMMON.GEO'
2330 include 'COMMON.LOCAL'
2331 include 'COMMON.INTERACT'
2332 include 'COMMON.IOUNITS'
2333 include 'COMMON.NAMES'
2334 double precision stochforcvec(MAXRES6)
2335 common /stochcalc/ stochforcvec
2336 logical lprn /.false./
2338 c write (iout,*) "dc_old"
2340 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2341 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2344 dc_work(j)=dc_old(j,0)
2345 d_t_work(j)=d_t_old(j,0)
2346 d_a_work(j)=d_a_old(j,0)
2351 dc_work(ind+j)=dc_old(j,i)
2352 d_t_work(ind+j)=d_t_old(j,i)
2353 d_a_work(ind+j)=d_a_old(j,i)
2358 if (itype(i).ne.10) then
2360 dc_work(ind+j)=dc_old(j,i+nres)
2361 d_t_work(ind+j)=d_t_old(j,i+nres)
2362 d_a_work(ind+j)=d_a_old(j,i+nres)
2370 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2374 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2375 & vfric_mat(i,j),afric_mat(i,j),
2376 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2384 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2385 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2386 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2387 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2389 d_t_work_new(i)=ddt1+0.5d0*ddt2
2390 d_t_work(i)=ddt1+ddt2
2395 d_t(j,0)=d_t_work(j)
2400 dc(j,i)=dc_work(ind+j)
2401 d_t(j,i)=d_t_work(ind+j)
2406 if (itype(i).ne.10) then
2409 dc(j,inres)=dc_work(ind+j)
2410 d_t(j,inres)=d_t_work(ind+j)
2417 c--------------------------------------------------------------------------
2418 subroutine sd_verlet2
2419 c Calculating the adjusted velocities for accelerations
2420 implicit real*8 (a-h,o-z)
2421 include 'DIMENSIONS'
2422 include 'COMMON.CONTROL'
2423 include 'COMMON.VAR'
2426 include 'COMMON.LANGEVIN'
2428 include 'COMMON.LANGEVIN.lang0'
2430 include 'COMMON.CHAIN'
2431 include 'COMMON.DERIV'
2432 include 'COMMON.GEO'
2433 include 'COMMON.LOCAL'
2434 include 'COMMON.INTERACT'
2435 include 'COMMON.IOUNITS'
2436 include 'COMMON.NAMES'
2437 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2438 common /stochcalc/ stochforcvec
2440 c Compute the stochastic forces which contribute to velocity change
2442 call stochastic_force(stochforcvecV)
2449 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2450 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2451 & vrand_mat2(i,j)*stochforcvecV(j)
2453 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2457 d_t(j,0)=d_t_work(j)
2462 d_t(j,i)=d_t_work(ind+j)
2467 if (itype(i).ne.10) then
2470 d_t(j,inres)=d_t_work(ind+j)
2477 c-----------------------------------------------------------
2478 subroutine sd_verlet_ciccotti_setup
2479 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2481 implicit real*8 (a-h,o-z)
2482 include 'DIMENSIONS'
2486 include 'COMMON.CONTROL'
2487 include 'COMMON.VAR'
2490 include 'COMMON.LANGEVIN'
2492 include 'COMMON.LANGEVIN.lang0'
2494 include 'COMMON.CHAIN'
2495 include 'COMMON.DERIV'
2496 include 'COMMON.GEO'
2497 include 'COMMON.LOCAL'
2498 include 'COMMON.INTERACT'
2499 include 'COMMON.IOUNITS'
2500 include 'COMMON.NAMES'
2501 include 'COMMON.TIME1'
2502 double precision emgdt(MAXRES6),
2503 & pterm,vterm,rho,rhoc,vsig,
2504 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2505 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2506 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2507 logical lprn /.false./
2508 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2509 double precision ktm
2516 c AL 8/17/04 Code adapted from tinker
2518 c Get the frictional and random terms for stochastic dynamics in the
2519 c eigenspace of mass-scaled UNRES friction matrix
2522 write (iout,*) "i",i," fricgam",fricgam(i)
2523 gdt = fricgam(i) * d_time
2525 c Stochastic dynamics reduces to simple MD for zero friction
2527 if (gdt .le. zero) then
2528 pfric_vec(i) = 1.0d0
2529 vfric_vec(i) = d_time
2530 afric_vec(i) = 0.5d0*d_time*d_time
2531 prand_vec(i) = afric_vec(i)
2532 vrand_vec2(i) = vfric_vec(i)
2534 c Analytical expressions when friction coefficient is large
2539 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2540 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2541 prand_vec(i) = afric_vec(i)
2542 vrand_vec2(i) = vfric_vec(i)
2544 c Compute the scaling factors of random terms for the nonzero friction case
2546 c ktm = 0.5d0*d_time/fricgam(i)
2547 c psig = dsqrt(ktm*pterm) / fricgam(i)
2548 c vsig = dsqrt(ktm*vterm)
2549 c prand_vec(i) = psig*afric_vec(i)
2550 c vrand_vec2(i) = vsig*vfric_vec(i)
2555 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2558 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2559 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2563 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2565 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2566 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2567 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2568 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2569 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2571 t_sdsetup=t_sdsetup+MPI_Wtime()
2573 t_sdsetup=t_sdsetup+tcpu()-tt0
2577 c-------------------------------------------------------------
2578 subroutine sd_verlet1_ciccotti
2579 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2580 implicit real*8 (a-h,o-z)
2581 include 'DIMENSIONS'
2585 include 'COMMON.CONTROL'
2586 include 'COMMON.VAR'
2589 include 'COMMON.LANGEVIN'
2591 include 'COMMON.LANGEVIN.lang0'
2593 include 'COMMON.CHAIN'
2594 include 'COMMON.DERIV'
2595 include 'COMMON.GEO'
2596 include 'COMMON.LOCAL'
2597 include 'COMMON.INTERACT'
2598 include 'COMMON.IOUNITS'
2599 include 'COMMON.NAMES'
2600 double precision stochforcvec(MAXRES6)
2601 common /stochcalc/ stochforcvec
2602 logical lprn /.false./
2604 c write (iout,*) "dc_old"
2606 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2607 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2610 dc_work(j)=dc_old(j,0)
2611 d_t_work(j)=d_t_old(j,0)
2612 d_a_work(j)=d_a_old(j,0)
2617 dc_work(ind+j)=dc_old(j,i)
2618 d_t_work(ind+j)=d_t_old(j,i)
2619 d_a_work(ind+j)=d_a_old(j,i)
2624 if (itype(i).ne.10) then
2626 dc_work(ind+j)=dc_old(j,i+nres)
2627 d_t_work(ind+j)=d_t_old(j,i+nres)
2628 d_a_work(ind+j)=d_a_old(j,i+nres)
2637 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2641 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2642 & vfric_mat(i,j),afric_mat(i,j),
2643 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2651 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2652 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2653 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2654 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2656 d_t_work_new(i)=ddt1+0.5d0*ddt2
2657 d_t_work(i)=ddt1+ddt2
2662 d_t(j,0)=d_t_work(j)
2667 dc(j,i)=dc_work(ind+j)
2668 d_t(j,i)=d_t_work(ind+j)
2673 if (itype(i).ne.10) then
2676 dc(j,inres)=dc_work(ind+j)
2677 d_t(j,inres)=d_t_work(ind+j)
2684 c--------------------------------------------------------------------------
2685 subroutine sd_verlet2_ciccotti
2686 c Calculating the adjusted velocities for accelerations
2687 implicit real*8 (a-h,o-z)
2688 include 'DIMENSIONS'
2689 include 'COMMON.CONTROL'
2690 include 'COMMON.VAR'
2693 include 'COMMON.LANGEVIN'
2695 include 'COMMON.LANGEVIN.lang0'
2697 include 'COMMON.CHAIN'
2698 include 'COMMON.DERIV'
2699 include 'COMMON.GEO'
2700 include 'COMMON.LOCAL'
2701 include 'COMMON.INTERACT'
2702 include 'COMMON.IOUNITS'
2703 include 'COMMON.NAMES'
2704 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2705 common /stochcalc/ stochforcvec
2707 c Compute the stochastic forces which contribute to velocity change
2709 call stochastic_force(stochforcvecV)
2716 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2717 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2718 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2720 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2724 d_t(j,0)=d_t_work(j)
2729 d_t(j,i)=d_t_work(ind+j)
2734 if (itype(i).ne.10) then
2737 d_t(j,inres)=d_t_work(ind+j)
2745 c------------------------------------------------------
2746 double precision function HNose(ek,s,e,pi,Q,t_bath,dimenl)
2748 double precision ek,s,e,pi,Q,t_bath,Rb
2751 HNose=ek+e+pi**2/(2*Q)+dimenl*Rb*t_bath*log(s)
2752 c print '(6f15.5,i5,a2,2f15.5)',ek,s,e,pi,Q,t_bath,dimenl,"--",
2753 c & pi**2/(2*Q),dimenl*Rb*t_bath*log(s)
2756 c-----------------------------------------------------------------
2757 double precision function HNose_nh(eki,e)
2758 implicit real*8 (a-h,o-z)
2759 include 'DIMENSIONS'
2761 HNose_nh=eki+e+dimen3*Rb*t_bath*xlogs(1)+qmass(1)*vlogs(1)**2/2
2763 HNose_nh=HNose_nh+qmass(i)*vlogs(i)**2/2+Rb*t_bath*xlogs(i)
2765 c write(4,'(5e15.5)')
2766 c & vlogs(1),xlogs(1),HNose,eki,e
2769 c-----------------------------------------------------------------
2770 SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8)
2771 implicit real*8 (a-h,o-z)
2772 include 'DIMENSIONS'
2774 double precision akin,gnkt,dt,aa,gkt,scale
2775 double precision wdti(maxyosh),wdti2(maxyosh),
2776 & wdti4(maxyosh),wdti8(maxyosh)
2777 integer i,iresn,iyosh,inos,nnos1
2786 C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE
2787 C INTEGRATION FROM t=0 TO t=DT/2
2788 C GET THE TOTAL KINETIC ENERGY
2790 c CALL GETKINP(MASS,VX,VY,VZ,AKIN)
2792 GLOGS(1) = (AKIN - GNKT)/QMASS(1)
2793 C START THE MULTIPLE TIME STEP PROCEDURE
2796 C UPDATE THE THERMOSTAT VELOCITIES
2797 VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2799 AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) )
2800 VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA
2801 & + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA
2803 C UPDATE THE PARTICLE VELOCITIES
2804 AA = EXP(-WDTI2(IYOSH)*VLOGS(1) )
2807 GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1)
2808 C UPDATE THE THERMOSTAT POSITIONS
2810 XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH)
2812 C UPDATE THE THERMOSTAT VELOCITIES
2814 AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) )
2815 VLOGS(INOS) = VLOGS(INOS)*AA*AA
2816 & + WDTI4(IYOSH)*GLOGS(INOS)*AA
2817 GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS)
2818 & -GKT)/QMASS(INOS+1)
2820 VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2823 C UPDATE THE PARTICLE VELOCITIES
2824 c outside of this subroutine
2826 c VX(I) = VX(I)*SCALE
2827 c VY(I) = VY(I)*SCALE
2828 c VZ(I) = VZ(I)*SCALE
2832 c-----------------------------------------------------------------
2833 subroutine tnp1_respa_i_step1
2834 c Applying Nose-Poincare algorithm - step 1 to coordinates
2835 c JPSJ 70 75 (2001) S. Nose
2837 c d_t is not updated here
2839 implicit real*8 (a-h,o-z)
2840 include 'DIMENSIONS'
2841 include 'COMMON.CONTROL'
2842 include 'COMMON.VAR'
2844 include 'COMMON.CHAIN'
2845 include 'COMMON.DERIV'
2846 include 'COMMON.GEO'
2847 include 'COMMON.LOCAL'
2848 include 'COMMON.INTERACT'
2849 include 'COMMON.IOUNITS'
2850 include 'COMMON.NAMES'
2851 double precision adt,adt2,tmp
2853 tmp=1+pi_np/(2*Q_np)*0.5*d_time
2856 s12_dt=d_time/s12_np
2857 d_time_s12=d_time*0.5*s12_np
2860 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2861 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2865 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2866 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2870 if (itype(i).ne.10) then
2873 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2874 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2880 c---------------------------------------------------------------------
2881 subroutine tnp1_respa_i_step2
2882 c Step 2 of the velocity Verlet algorithm: update velocities
2883 implicit real*8 (a-h,o-z)
2884 include 'DIMENSIONS'
2885 include 'COMMON.CONTROL'
2886 include 'COMMON.VAR'
2888 include 'COMMON.CHAIN'
2889 include 'COMMON.DERIV'
2890 include 'COMMON.GEO'
2891 include 'COMMON.LOCAL'
2892 include 'COMMON.INTERACT'
2893 include 'COMMON.IOUNITS'
2894 include 'COMMON.NAMES'
2896 double precision d_time_s12
2900 d_t(j,i)=d_t_new(j,i)
2907 d_time_s12=0.5d0*s12_np*d_time
2910 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
2914 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
2918 if (itype(i).ne.10) then
2921 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
2926 pistar=pistar+(EK-0.5*(E_old+potE)
2927 & -dimen3*Rb*t_bath*log(s12_np)+Csplit-dimen3*Rb*t_bath)*d_time
2928 tmp=1+pistar/(2*Q_np)*0.5*d_time
2934 c-------------------------------------------------------
2936 subroutine tnp1_step1
2937 c Applying Nose-Poincare algorithm - step 1 to coordinates
2938 c JPSJ 70 75 (2001) S. Nose
2940 c d_t is not updated here
2942 implicit real*8 (a-h,o-z)
2943 include 'DIMENSIONS'
2944 include 'COMMON.CONTROL'
2945 include 'COMMON.VAR'
2947 include 'COMMON.CHAIN'
2948 include 'COMMON.DERIV'
2949 include 'COMMON.GEO'
2950 include 'COMMON.LOCAL'
2951 include 'COMMON.INTERACT'
2952 include 'COMMON.IOUNITS'
2953 include 'COMMON.NAMES'
2954 double precision adt,adt2,tmp
2956 tmp=1+pi_np/(2*Q_np)*0.5*d_time
2959 s12_dt=d_time/s12_np
2960 d_time_s12=d_time*0.5*s12_np
2963 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2964 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2968 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2969 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2973 if (itype(i).ne.10) then
2976 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2977 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2983 c---------------------------------------------------------------------
2984 subroutine tnp1_step2
2985 c Step 2 of the velocity Verlet algorithm: update velocities
2986 implicit real*8 (a-h,o-z)
2987 include 'DIMENSIONS'
2988 include 'COMMON.CONTROL'
2989 include 'COMMON.VAR'
2991 include 'COMMON.CHAIN'
2992 include 'COMMON.DERIV'
2993 include 'COMMON.GEO'
2994 include 'COMMON.LOCAL'
2995 include 'COMMON.INTERACT'
2996 include 'COMMON.IOUNITS'
2997 include 'COMMON.NAMES'
2999 double precision d_time_s12
3003 d_t(j,i)=d_t_new(j,i)
3010 d_time_s12=0.5d0*s12_np*d_time
3013 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
3017 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
3021 if (itype(i).ne.10) then
3024 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
3029 cd write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np
3030 pistar=pistar+(EK-0.5*(E_old+potE)
3031 & -dimen3*Rb*t_bath*log(s12_np)+H0-dimen3*Rb*t_bath)*d_time
3032 tmp=1+pistar/(2*Q_np)*0.5*d_time
3039 c-----------------------------------------------------------------
3040 subroutine tnp_respa_i_step1
3041 c Applying Nose-Poincare algorithm - step 1 to coordinates
3042 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3044 c d_t is not updated here, it is destroyed
3046 implicit real*8 (a-h,o-z)
3047 include 'DIMENSIONS'
3048 include 'COMMON.CONTROL'
3049 include 'COMMON.VAR'
3051 include 'COMMON.CHAIN'
3052 include 'COMMON.DERIV'
3053 include 'COMMON.GEO'
3054 include 'COMMON.LOCAL'
3055 include 'COMMON.INTERACT'
3056 include 'COMMON.IOUNITS'
3057 include 'COMMON.NAMES'
3058 double precision C_np,d_time_s,tmp,d_time_ss
3060 d_time_s=d_time*0.5*s_np
3061 ct2 d_time_s=d_time*0.5*s12_np
3064 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3068 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3072 if (itype(i).ne.10) then
3075 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3082 d_t(j,i)=d_t_new(j,i)
3089 C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-Csplit)
3092 pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3093 tmp=0.5*d_time*pistar/Q_np
3094 s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3096 d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3097 ct2 d_time_ss=d_time/s12_np
3098 c d_time_ss=0.5*d_time*(1.0/sold_np+1.0/s_np)
3101 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3105 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3109 if (itype(i).ne.10) then
3112 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3119 c---------------------------------------------------------------------
3121 subroutine tnp_respa_i_step2
3122 c Step 2 of the velocity Verlet algorithm: update velocities
3123 implicit real*8 (a-h,o-z)
3124 include 'DIMENSIONS'
3125 include 'COMMON.CONTROL'
3126 include 'COMMON.VAR'
3128 include 'COMMON.CHAIN'
3129 include 'COMMON.DERIV'
3130 include 'COMMON.GEO'
3131 include 'COMMON.LOCAL'
3132 include 'COMMON.INTERACT'
3133 include 'COMMON.IOUNITS'
3134 include 'COMMON.NAMES'
3136 double precision d_time_s
3138 EK=EK*(s_np/s12_np)**2
3139 HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3140 pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath
3143 cr print '(a,5f)','i_step2',EK,potE,HNose1,pi_np,E_long
3144 d_time_s=d_time*0.5*s12_np
3145 c d_time_s=d_time*0.5*s_np
3148 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3152 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3156 if (itype(i).ne.10) then
3159 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
3168 c-----------------------------------------------------------------
3169 subroutine tnp_respa_step1
3170 c Applying Nose-Poincare algorithm - step 1 to vel for RESPA
3171 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3173 c d_t is not updated here, it is destroyed
3175 implicit real*8 (a-h,o-z)
3176 include 'DIMENSIONS'
3177 include 'COMMON.CONTROL'
3178 include 'COMMON.VAR'
3180 include 'COMMON.CHAIN'
3181 include 'COMMON.DERIV'
3182 include 'COMMON.GEO'
3183 include 'COMMON.LOCAL'
3184 include 'COMMON.INTERACT'
3185 include 'COMMON.IOUNITS'
3186 include 'COMMON.NAMES'
3187 double precision C_np,d_time_s,tmp,d_time_ss
3188 double precision energia(0:n_ene)
3190 d_time_s=d_time*0.5*s_np
3193 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3197 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3201 if (itype(i).ne.10) then
3204 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3210 c C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3213 c pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3214 c tmp=0.5*d_time*pistar/Q_np
3215 c s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3216 c write(iout,*) 'tnp_respa_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3222 c-------------------------------------
3223 c test of reviewer's comment
3224 pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3225 cr print '(a,3f)','1 pi_np,s_np',pi_np,s_np,E_long
3226 c-------------------------------------
3230 c---------------------------------------------------------------------
3231 subroutine tnp_respa_step2
3232 c Step 2 of the velocity Verlet algorithm: update velocities for RESPA
3233 implicit real*8 (a-h,o-z)
3234 include 'DIMENSIONS'
3235 include 'COMMON.CONTROL'
3236 include 'COMMON.VAR'
3238 include 'COMMON.CHAIN'
3239 include 'COMMON.DERIV'
3240 include 'COMMON.GEO'
3241 include 'COMMON.LOCAL'
3242 include 'COMMON.INTERACT'
3243 include 'COMMON.IOUNITS'
3244 include 'COMMON.NAMES'
3246 double precision d_time_s
3252 ct HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3253 ct pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3254 ct & -0.5*d_time*(HNose1-H0)
3256 c-------------------------------------
3257 c test of reviewer's comment
3258 pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3259 cr print '(a,3f)','2 pi_np,s_np',pi_np,s_np,E_long
3260 c-------------------------------------
3261 d_time_s=d_time*0.5*s_np
3264 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3268 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3272 if (itype(i).ne.10) then
3275 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3284 c---------------------------------------------------------------------
3285 subroutine tnp_step1
3286 c Applying Nose-Poincare algorithm - step 1 to coordinates
3287 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3289 c d_t is not updated here, it is destroyed
3291 implicit real*8 (a-h,o-z)
3292 include 'DIMENSIONS'
3293 include 'COMMON.CONTROL'
3294 include 'COMMON.VAR'
3296 include 'COMMON.CHAIN'
3297 include 'COMMON.DERIV'
3298 include 'COMMON.GEO'
3299 include 'COMMON.LOCAL'
3300 include 'COMMON.INTERACT'
3301 include 'COMMON.IOUNITS'
3302 include 'COMMON.NAMES'
3303 double precision C_np,d_time_s,tmp,d_time_ss
3305 d_time_s=d_time*0.5*s_np
3308 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3312 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3316 if (itype(i).ne.10) then
3319 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3326 d_t(j,i)=d_t_new(j,i)
3333 C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3336 pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3337 tmp=0.5*d_time*pistar/Q_np
3338 s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3339 c write(iout,*) 'tnp_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3341 d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3344 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3348 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3352 if (itype(i).ne.10) then
3355 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3362 c-----------------------------------------------------------------
3363 subroutine tnp_step2
3364 c Step 2 of the velocity Verlet algorithm: update velocities
3365 implicit real*8 (a-h,o-z)
3366 include 'DIMENSIONS'
3367 include 'COMMON.CONTROL'
3368 include 'COMMON.VAR'
3370 include 'COMMON.CHAIN'
3371 include 'COMMON.DERIV'
3372 include 'COMMON.GEO'
3373 include 'COMMON.LOCAL'
3374 include 'COMMON.INTERACT'
3375 include 'COMMON.IOUNITS'
3376 include 'COMMON.NAMES'
3378 double precision d_time_s
3380 EK=EK*(s_np/s12_np)**2
3381 HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3382 pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3383 & -0.5*d_time*(HNose1-H0)
3385 cd write(iout,'(a,4f)') 'mmm',EK,potE,HNose1,pi_np
3386 d_time_s=d_time*0.5*s12_np
3389 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3393 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3397 if (itype(i).ne.10) then
3400 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
3410 subroutine hmc_test(itime)
3411 implicit real*8 (a-h,o-z)
3412 include 'DIMENSIONS'
3413 include 'COMMON.CONTROL'
3415 include 'COMMON.CHAIN'
3418 delta=-(potE+EK-hmc_etot)/(Rb*t_bath)
3419 if (delta .lt. -50.0d0) then
3424 xxx=ran_number(0.0d0,1.0d0)
3426 if (me.eq.king .or. .not. out1file)
3427 & write(iout,'(a8,i5,6f10.4)')
3428 & 'HMC',itime,potE+EK,potE,EK,hmc_etot,delta,xxx
3430 if (delta .le. xxx) then
3439 if (me.eq.king .or. .not. out1file)
3440 & write(iout,*) 'HMC accepting new'
3449 call chainbuild_cart
3453 d_t_old(j,i)=d_t(j,i)
3457 kinetic_T=2.0d0/(dimen3*Rb)*EK
3458 call etotal(potEcomp)
3461 if (me.eq.king .or. .not. out1file)
3462 & write(iout,'(a8,i5,3f10.4)')'HMC new',itime,potE+EK,potE,EK