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,7f)') "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(iabs(itype(i)))
1685 & *dsqrt(gamsc(iabs(itype(i))))
1688 c Open the pdb file for snapshotshots
1691 if (ilen(tmpdir).gt.0)
1692 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1693 & liczba(:ilen(liczba))//".pdb")
1695 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1699 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1700 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1701 & liczba(:ilen(liczba))//".x")
1702 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1705 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1706 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1707 & liczba(:ilen(liczba))//".cx")
1708 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1714 if (ilen(tmpdir).gt.0)
1715 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1716 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1718 if (ilen(tmpdir).gt.0)
1719 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1720 cartname=prefix(:ilen(prefix))//"_MD.cx"
1724 write (qstr,'(256(1h ))')
1727 iq = qinfrag(i,iset)*10
1728 iw = wfrag(i,iset)/100
1730 if(me.eq.king.or..not.out1file)
1731 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1732 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1737 iq = qinpair(i,iset)*10
1738 iw = wpair(i,iset)/100
1740 if(me.eq.king.or..not.out1file)
1741 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1742 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1746 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1748 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1750 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1752 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1756 if (restart1file) then
1758 & inquire(file=mremd_rst_name,exist=file_exist)
1759 write (*,*) me," Before broadcast: file_exist",file_exist
1761 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1763 write (*,*) me," After broadcast: file_exist",file_exist
1765 c inquire(file=mremd_rst_name,exist=file_exist)
1766 if(me.eq.king.or..not.out1file)
1767 & write(iout,*) "Initial state read by master and distributed"
1769 if (ilen(tmpdir).gt.0)
1770 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1771 & //liczba(:ilen(liczba))//'.rst')
1772 inquire(file=rest2name,exist=file_exist)
1775 if(.not.restart1file) then
1776 if(me.eq.king.or..not.out1file)
1777 & write(iout,*) "Initial state will be read from file ",
1778 & rest2name(:ilen(rest2name))
1781 call rescale_weights(t_bath)
1783 if(me.eq.king.or..not.out1file)then
1784 if (restart1file) then
1785 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1788 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1791 write(iout,*) "Initial velocities randomly generated"
1797 c Generate initial velocities
1798 if(me.eq.king.or..not.out1file)
1799 & write(iout,*) "Initial velocities randomly generated"
1803 c rest2name = prefix(:ilen(prefix))//'.rst'
1804 if(me.eq.king.or..not.out1file)then
1805 write (iout,*) "Initial velocities"
1807 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1808 & (d_t(j,i+nres),j=1,3)
1811 c Zeroing the total angular momentum of the system
1812 write(iout,*) "Calling the zero-angular
1813 & momentum subroutine"
1816 c Getting the potential energy and forces and velocities and accelerations
1818 c write (iout,*) "velocity of the center of the mass:"
1819 c write (iout,*) (vcm(j),j=1,3)
1821 d_t(j,0)=d_t(j,0)-vcm(j)
1823 c Removing the velocity of the center of mass
1825 if(me.eq.king.or..not.out1file)then
1826 write (iout,*) "vcm right after adjustment:"
1827 write (iout,*) (vcm(j),j=1,3)
1832 if(iranconf.ne.0) 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)
1855 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1858 call chainbuild_cart
1861 call verlet_bath(EK)
1863 kinetic_T=2.0d0/(dimen3*Rb)*EK
1864 if(me.eq.king.or..not.out1file)then
1874 call etotal(potEcomp)
1877 t_etotal=t_etotal+MPI_Wtime()-tt0
1879 t_etotal=t_etotal+tcpu()-tt0
1884 if(tnp .or. tnp1) then
1887 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
1889 write(iout,*) 'H0= ',H0
1893 HNose1=Hnose_nh(EK,potE)
1895 write (iout,*) 'H0= ',H0
1901 if(me.eq.king.or..not.out1file)
1902 & write(iout,*) 'HMC',hmc_etot,potE,EK
1913 if (amax*d_time .gt. dvmax) then
1914 d_time=d_time*dvmax/amax
1915 if(me.eq.king.or..not.out1file) write (iout,*)
1916 & "Time step reduced to",d_time,
1917 & " because of too large initial acceleration."
1919 if(me.eq.king.or..not.out1file)then
1920 write(iout,*) "Potential energy and its components"
1921 call enerprint(potEcomp)
1922 c write(iout,*) (potEcomp(i),i=0,n_ene)
1924 potE=potEcomp(0)-potEcomp(20)
1927 if (ntwe.ne.0) call statout(itime)
1928 if(me.eq.king.or..not.out1file)
1929 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1930 & " Kinetic energy",EK," potential energy",potE,
1931 & " total energy",totE," maximum acceleration ",
1934 write (iout,*) "Initial coordinates"
1936 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1937 & (c(j,i+nres),j=1,3)
1939 write (iout,*) "Initial dC"
1941 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1942 & (dc(j,i+nres),j=1,3)
1944 write (iout,*) "Initial velocities"
1945 write (iout,"(13x,' backbone ',23x,' side chain')")
1947 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1948 & (d_t(j,i+nres),j=1,3)
1950 write (iout,*) "Initial accelerations"
1952 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1953 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1954 & (d_a(j,i+nres),j=1,3)
1960 d_t_old(j,i)=d_t(j,i)
1961 d_a_old(j,i)=d_a(j,i)
1963 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1972 call etotal_short(energia_short)
1975 t_eshort=t_eshort+MPI_Wtime()-tt0
1977 t_eshort=t_eshort+tcpu()-tt0
1981 if(tnp .or. tnp1) then
1982 E_short=energia_short(0)
1983 HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen3)
1986 c_new_var_csplit Csplit=H0-E_long
1987 c Csplit = H0-energia_short(0)
1988 write(iout,*) 'Csplit= ',Csplit
1994 if(.not.out1file .and. large) then
1995 write (iout,*) "energia_long",energia_long(0),
1996 & " energia_short",energia_short(0),
1997 & " total",energia_long(0)+energia_short(0)
1998 write (iout,*) "Initial fast-force accelerations"
2000 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2001 & (d_a(j,i+nres),j=1,3)
2004 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2007 d_a_short(j,i)=d_a(j,i)
2016 call etotal_long(energia_long)
2019 t_elong=t_elong+MPI_Wtime()-tt0
2021 t_elong=t_elong+tcpu()-tt0
2026 if(.not.out1file .and. large) then
2027 write (iout,*) "energia_long",energia_long(0)
2028 write (iout,*) "Initial slow-force accelerations"
2030 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2031 & (d_a(j,i+nres),j=1,3)
2035 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2037 t_enegrad=t_enegrad+tcpu()-tt0
2045 c-----------------------------------------------------------
2046 subroutine random_vel
2047 implicit real*8 (a-h,o-z)
2048 include 'DIMENSIONS'
2049 include 'COMMON.CONTROL'
2050 include 'COMMON.VAR'
2053 include 'COMMON.LANGEVIN'
2055 include 'COMMON.LANGEVIN.lang0'
2057 include 'COMMON.CHAIN'
2058 include 'COMMON.DERIV'
2059 include 'COMMON.GEO'
2060 include 'COMMON.LOCAL'
2061 include 'COMMON.INTERACT'
2062 include 'COMMON.IOUNITS'
2063 include 'COMMON.NAMES'
2064 include 'COMMON.TIME1'
2065 double precision xv,sigv,lowb,highb
2066 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2067 c First generate velocities in the eigenspace of the G matrix
2068 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2070 c write (iout,*) "RANDOM_VEL dimen",dimen
2076 sigv=dsqrt((Rb*t_bath)/geigen(i))
2079 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2080 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2081 c & " d_t_work_new",d_t_work_new(ii)
2091 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2094 c write (iout,*) "Ek from eigenvectors",Ek1
2096 c Transform velocities to UNRES coordinate space
2102 d_t_work(ind)=d_t_work(ind)
2103 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2105 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2109 c Transfer to the d_t vector
2111 d_t(j,0)=d_t_work(j)
2117 d_t(j,i)=d_t_work(ind)
2121 c write (iout,*) "d_t",i,(d_t(j,i),j=1,3)
2125 if (itype(i).ne.10) then
2128 d_t(j,i+nres)=d_t_work(ind)
2133 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2134 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2139 c-----------------------------------------------------------
2140 subroutine sd_verlet_p_setup
2141 c Sets up the parameters of stochastic Verlet algorithm
2142 implicit real*8 (a-h,o-z)
2143 include 'DIMENSIONS'
2147 include 'COMMON.CONTROL'
2148 include 'COMMON.VAR'
2151 include 'COMMON.LANGEVIN'
2153 include 'COMMON.LANGEVIN.lang0'
2155 include 'COMMON.CHAIN'
2156 include 'COMMON.DERIV'
2157 include 'COMMON.GEO'
2158 include 'COMMON.LOCAL'
2159 include 'COMMON.INTERACT'
2160 include 'COMMON.IOUNITS'
2161 include 'COMMON.NAMES'
2162 include 'COMMON.TIME1'
2163 double precision emgdt(MAXRES6),
2164 & pterm,vterm,rho,rhoc,vsig,
2165 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2166 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2167 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2168 logical lprn /.false./
2169 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2170 double precision ktm
2177 c AL 8/17/04 Code adapted from tinker
2179 c Get the frictional and random terms for stochastic dynamics in the
2180 c eigenspace of mass-scaled UNRES friction matrix
2183 gdt = fricgam(i) * d_time
2185 c Stochastic dynamics reduces to simple MD for zero friction
2187 if (gdt .le. zero) then
2188 pfric_vec(i) = 1.0d0
2189 vfric_vec(i) = d_time
2190 afric_vec(i) = 0.5d0 * d_time * d_time
2191 prand_vec(i) = 0.0d0
2192 vrand_vec1(i) = 0.0d0
2193 vrand_vec2(i) = 0.0d0
2195 c Analytical expressions when friction coefficient is large
2198 if (gdt .ge. gdt_radius) then
2201 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2202 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2203 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2204 vterm = 1.0d0 - egdt**2
2205 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2207 c Use series expansions when friction coefficient is small
2218 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2219 & - gdt5/120.0d0 + gdt6/720.0d0
2220 & - gdt7/5040.0d0 + gdt8/40320.0d0
2221 & - gdt9/362880.0d0) / fricgam(i)**2
2222 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2223 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2224 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2225 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2226 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2227 & + 127.0d0*gdt9/90720.0d0
2228 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2229 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2230 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2231 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2232 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2233 & - 17.0d0*gdt2/1280.0d0
2234 & + 17.0d0*gdt3/6144.0d0
2235 & + 40967.0d0*gdt4/34406400.0d0
2236 & - 57203.0d0*gdt5/275251200.0d0
2237 & - 1429487.0d0*gdt6/13212057600.0d0)
2240 c Compute the scaling factors of random terms for the nonzero friction case
2242 ktm = 0.5d0*d_time/fricgam(i)
2243 psig = dsqrt(ktm*pterm) / fricgam(i)
2244 vsig = dsqrt(ktm*vterm)
2245 rhoc = dsqrt(1.0d0 - rho*rho)
2247 vrand_vec1(i) = vsig * rho
2248 vrand_vec2(i) = vsig * rhoc
2253 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2256 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2257 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2261 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2264 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2265 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2266 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2267 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2268 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2269 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2272 t_sdsetup=t_sdsetup+MPI_Wtime()
2274 t_sdsetup=t_sdsetup+tcpu()-tt0
2278 c-------------------------------------------------------------
2279 subroutine eigtransf1(n,ndim,ab,d,c)
2282 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2288 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2294 c-------------------------------------------------------------
2295 subroutine eigtransf(n,ndim,a,b,d,c)
2298 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2304 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2310 c-------------------------------------------------------------
2311 subroutine sd_verlet1
2312 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2313 implicit real*8 (a-h,o-z)
2314 include 'DIMENSIONS'
2315 include 'COMMON.CONTROL'
2316 include 'COMMON.VAR'
2319 include 'COMMON.LANGEVIN'
2321 include 'COMMON.LANGEVIN.lang0'
2323 include 'COMMON.CHAIN'
2324 include 'COMMON.DERIV'
2325 include 'COMMON.GEO'
2326 include 'COMMON.LOCAL'
2327 include 'COMMON.INTERACT'
2328 include 'COMMON.IOUNITS'
2329 include 'COMMON.NAMES'
2330 double precision stochforcvec(MAXRES6)
2331 common /stochcalc/ stochforcvec
2332 logical lprn /.false./
2334 c write (iout,*) "dc_old"
2336 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2337 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2340 dc_work(j)=dc_old(j,0)
2341 d_t_work(j)=d_t_old(j,0)
2342 d_a_work(j)=d_a_old(j,0)
2347 dc_work(ind+j)=dc_old(j,i)
2348 d_t_work(ind+j)=d_t_old(j,i)
2349 d_a_work(ind+j)=d_a_old(j,i)
2354 if (itype(i).ne.10) then
2356 dc_work(ind+j)=dc_old(j,i+nres)
2357 d_t_work(ind+j)=d_t_old(j,i+nres)
2358 d_a_work(ind+j)=d_a_old(j,i+nres)
2366 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2370 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2371 & vfric_mat(i,j),afric_mat(i,j),
2372 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2380 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2381 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2382 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2383 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2385 d_t_work_new(i)=ddt1+0.5d0*ddt2
2386 d_t_work(i)=ddt1+ddt2
2391 d_t(j,0)=d_t_work(j)
2396 dc(j,i)=dc_work(ind+j)
2397 d_t(j,i)=d_t_work(ind+j)
2402 if (itype(i).ne.10) then
2405 dc(j,inres)=dc_work(ind+j)
2406 d_t(j,inres)=d_t_work(ind+j)
2413 c--------------------------------------------------------------------------
2414 subroutine sd_verlet2
2415 c Calculating the adjusted velocities for accelerations
2416 implicit real*8 (a-h,o-z)
2417 include 'DIMENSIONS'
2418 include 'COMMON.CONTROL'
2419 include 'COMMON.VAR'
2422 include 'COMMON.LANGEVIN'
2424 include 'COMMON.LANGEVIN.lang0'
2426 include 'COMMON.CHAIN'
2427 include 'COMMON.DERIV'
2428 include 'COMMON.GEO'
2429 include 'COMMON.LOCAL'
2430 include 'COMMON.INTERACT'
2431 include 'COMMON.IOUNITS'
2432 include 'COMMON.NAMES'
2433 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2434 common /stochcalc/ stochforcvec
2436 c Compute the stochastic forces which contribute to velocity change
2438 call stochastic_force(stochforcvecV)
2445 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2446 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2447 & vrand_mat2(i,j)*stochforcvecV(j)
2449 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2453 d_t(j,0)=d_t_work(j)
2458 d_t(j,i)=d_t_work(ind+j)
2463 if (itype(i).ne.10) then
2466 d_t(j,inres)=d_t_work(ind+j)
2473 c-----------------------------------------------------------
2474 subroutine sd_verlet_ciccotti_setup
2475 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2477 implicit real*8 (a-h,o-z)
2478 include 'DIMENSIONS'
2482 include 'COMMON.CONTROL'
2483 include 'COMMON.VAR'
2486 include 'COMMON.LANGEVIN'
2488 include 'COMMON.LANGEVIN.lang0'
2490 include 'COMMON.CHAIN'
2491 include 'COMMON.DERIV'
2492 include 'COMMON.GEO'
2493 include 'COMMON.LOCAL'
2494 include 'COMMON.INTERACT'
2495 include 'COMMON.IOUNITS'
2496 include 'COMMON.NAMES'
2497 include 'COMMON.TIME1'
2498 double precision emgdt(MAXRES6),
2499 & pterm,vterm,rho,rhoc,vsig,
2500 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2501 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2502 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2503 logical lprn /.false./
2504 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2505 double precision ktm
2512 c AL 8/17/04 Code adapted from tinker
2514 c Get the frictional and random terms for stochastic dynamics in the
2515 c eigenspace of mass-scaled UNRES friction matrix
2518 write (iout,*) "i",i," fricgam",fricgam(i)
2519 gdt = fricgam(i) * d_time
2521 c Stochastic dynamics reduces to simple MD for zero friction
2523 if (gdt .le. zero) then
2524 pfric_vec(i) = 1.0d0
2525 vfric_vec(i) = d_time
2526 afric_vec(i) = 0.5d0*d_time*d_time
2527 prand_vec(i) = afric_vec(i)
2528 vrand_vec2(i) = vfric_vec(i)
2530 c Analytical expressions when friction coefficient is large
2535 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2536 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2537 prand_vec(i) = afric_vec(i)
2538 vrand_vec2(i) = vfric_vec(i)
2540 c Compute the scaling factors of random terms for the nonzero friction case
2542 c ktm = 0.5d0*d_time/fricgam(i)
2543 c psig = dsqrt(ktm*pterm) / fricgam(i)
2544 c vsig = dsqrt(ktm*vterm)
2545 c prand_vec(i) = psig*afric_vec(i)
2546 c vrand_vec2(i) = vsig*vfric_vec(i)
2551 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2554 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2555 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2559 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2561 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2562 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2563 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2564 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2565 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2567 t_sdsetup=t_sdsetup+MPI_Wtime()
2569 t_sdsetup=t_sdsetup+tcpu()-tt0
2573 c-------------------------------------------------------------
2574 subroutine sd_verlet1_ciccotti
2575 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2576 implicit real*8 (a-h,o-z)
2577 include 'DIMENSIONS'
2581 include 'COMMON.CONTROL'
2582 include 'COMMON.VAR'
2585 include 'COMMON.LANGEVIN'
2587 include 'COMMON.LANGEVIN.lang0'
2589 include 'COMMON.CHAIN'
2590 include 'COMMON.DERIV'
2591 include 'COMMON.GEO'
2592 include 'COMMON.LOCAL'
2593 include 'COMMON.INTERACT'
2594 include 'COMMON.IOUNITS'
2595 include 'COMMON.NAMES'
2596 double precision stochforcvec(MAXRES6)
2597 common /stochcalc/ stochforcvec
2598 logical lprn /.false./
2600 c write (iout,*) "dc_old"
2602 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2603 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2606 dc_work(j)=dc_old(j,0)
2607 d_t_work(j)=d_t_old(j,0)
2608 d_a_work(j)=d_a_old(j,0)
2613 dc_work(ind+j)=dc_old(j,i)
2614 d_t_work(ind+j)=d_t_old(j,i)
2615 d_a_work(ind+j)=d_a_old(j,i)
2620 if (itype(i).ne.10) then
2622 dc_work(ind+j)=dc_old(j,i+nres)
2623 d_t_work(ind+j)=d_t_old(j,i+nres)
2624 d_a_work(ind+j)=d_a_old(j,i+nres)
2633 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2637 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2638 & vfric_mat(i,j),afric_mat(i,j),
2639 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2647 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2648 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2649 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2650 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2652 d_t_work_new(i)=ddt1+0.5d0*ddt2
2653 d_t_work(i)=ddt1+ddt2
2658 d_t(j,0)=d_t_work(j)
2663 dc(j,i)=dc_work(ind+j)
2664 d_t(j,i)=d_t_work(ind+j)
2669 if (itype(i).ne.10) then
2672 dc(j,inres)=dc_work(ind+j)
2673 d_t(j,inres)=d_t_work(ind+j)
2680 c--------------------------------------------------------------------------
2681 subroutine sd_verlet2_ciccotti
2682 c Calculating the adjusted velocities for accelerations
2683 implicit real*8 (a-h,o-z)
2684 include 'DIMENSIONS'
2685 include 'COMMON.CONTROL'
2686 include 'COMMON.VAR'
2689 include 'COMMON.LANGEVIN'
2691 include 'COMMON.LANGEVIN.lang0'
2693 include 'COMMON.CHAIN'
2694 include 'COMMON.DERIV'
2695 include 'COMMON.GEO'
2696 include 'COMMON.LOCAL'
2697 include 'COMMON.INTERACT'
2698 include 'COMMON.IOUNITS'
2699 include 'COMMON.NAMES'
2700 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2701 common /stochcalc/ stochforcvec
2703 c Compute the stochastic forces which contribute to velocity change
2705 call stochastic_force(stochforcvecV)
2712 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2713 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2714 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2716 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2720 d_t(j,0)=d_t_work(j)
2725 d_t(j,i)=d_t_work(ind+j)
2730 if (itype(i).ne.10) then
2733 d_t(j,inres)=d_t_work(ind+j)
2741 c------------------------------------------------------
2742 double precision function HNose(ek,s,e,pi,Q,t_bath,dimenl)
2744 double precision ek,s,e,pi,Q,t_bath,Rb
2747 HNose=ek+e+pi**2/(2*Q)+dimenl*Rb*t_bath*log(s)
2748 c print '(6f15.5,i5,a2,2f15.5)',ek,s,e,pi,Q,t_bath,dimenl,"--",
2749 c & pi**2/(2*Q),dimenl*Rb*t_bath*log(s)
2752 c-----------------------------------------------------------------
2753 double precision function HNose_nh(eki,e)
2754 implicit real*8 (a-h,o-z)
2755 include 'DIMENSIONS'
2757 HNose_nh=eki+e+dimen3*Rb*t_bath*xlogs(1)+qmass(1)*vlogs(1)**2/2
2759 HNose_nh=HNose_nh+qmass(i)*vlogs(i)**2/2+Rb*t_bath*xlogs(i)
2761 c write(4,'(5e15.5)')
2762 c & vlogs(1),xlogs(1),HNose,eki,e
2765 c-----------------------------------------------------------------
2766 SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8)
2767 implicit real*8 (a-h,o-z)
2768 include 'DIMENSIONS'
2770 double precision akin,gnkt,dt,aa,gkt,scale
2771 double precision wdti(maxyosh),wdti2(maxyosh),
2772 & wdti4(maxyosh),wdti8(maxyosh)
2773 integer i,iresn,iyosh,inos,nnos1
2782 C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE
2783 C INTEGRATION FROM t=0 TO t=DT/2
2784 C GET THE TOTAL KINETIC ENERGY
2786 c CALL GETKINP(MASS,VX,VY,VZ,AKIN)
2788 GLOGS(1) = (AKIN - GNKT)/QMASS(1)
2789 C START THE MULTIPLE TIME STEP PROCEDURE
2792 C UPDATE THE THERMOSTAT VELOCITIES
2793 VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2795 AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) )
2796 VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA
2797 & + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA
2799 C UPDATE THE PARTICLE VELOCITIES
2800 AA = EXP(-WDTI2(IYOSH)*VLOGS(1) )
2803 GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1)
2804 C UPDATE THE THERMOSTAT POSITIONS
2806 XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH)
2808 C UPDATE THE THERMOSTAT VELOCITIES
2810 AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) )
2811 VLOGS(INOS) = VLOGS(INOS)*AA*AA
2812 & + WDTI4(IYOSH)*GLOGS(INOS)*AA
2813 GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS)
2814 & -GKT)/QMASS(INOS+1)
2816 VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2819 C UPDATE THE PARTICLE VELOCITIES
2820 c outside of this subroutine
2822 c VX(I) = VX(I)*SCALE
2823 c VY(I) = VY(I)*SCALE
2824 c VZ(I) = VZ(I)*SCALE
2828 c-----------------------------------------------------------------
2829 subroutine tnp1_respa_i_step1
2830 c Applying Nose-Poincare algorithm - step 1 to coordinates
2831 c JPSJ 70 75 (2001) S. Nose
2833 c d_t is not updated here
2835 implicit real*8 (a-h,o-z)
2836 include 'DIMENSIONS'
2837 include 'COMMON.CONTROL'
2838 include 'COMMON.VAR'
2840 include 'COMMON.CHAIN'
2841 include 'COMMON.DERIV'
2842 include 'COMMON.GEO'
2843 include 'COMMON.LOCAL'
2844 include 'COMMON.INTERACT'
2845 include 'COMMON.IOUNITS'
2846 include 'COMMON.NAMES'
2847 double precision adt,adt2,tmp
2849 tmp=1+pi_np/(2*Q_np)*0.5*d_time
2852 s12_dt=d_time/s12_np
2853 d_time_s12=d_time*0.5*s12_np
2856 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2857 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2861 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2862 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2866 if (itype(i).ne.10) then
2869 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2870 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2876 c---------------------------------------------------------------------
2877 subroutine tnp1_respa_i_step2
2878 c Step 2 of the velocity Verlet algorithm: update velocities
2879 implicit real*8 (a-h,o-z)
2880 include 'DIMENSIONS'
2881 include 'COMMON.CONTROL'
2882 include 'COMMON.VAR'
2884 include 'COMMON.CHAIN'
2885 include 'COMMON.DERIV'
2886 include 'COMMON.GEO'
2887 include 'COMMON.LOCAL'
2888 include 'COMMON.INTERACT'
2889 include 'COMMON.IOUNITS'
2890 include 'COMMON.NAMES'
2892 double precision d_time_s12
2896 d_t(j,i)=d_t_new(j,i)
2903 d_time_s12=0.5d0*s12_np*d_time
2906 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
2910 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
2914 if (itype(i).ne.10) then
2917 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
2922 pistar=pistar+(EK-0.5*(E_old+potE)
2923 & -dimen3*Rb*t_bath*log(s12_np)+Csplit-dimen3*Rb*t_bath)*d_time
2924 tmp=1+pistar/(2*Q_np)*0.5*d_time
2930 c-------------------------------------------------------
2932 subroutine tnp1_step1
2933 c Applying Nose-Poincare algorithm - step 1 to coordinates
2934 c JPSJ 70 75 (2001) S. Nose
2936 c d_t is not updated here
2938 implicit real*8 (a-h,o-z)
2939 include 'DIMENSIONS'
2940 include 'COMMON.CONTROL'
2941 include 'COMMON.VAR'
2943 include 'COMMON.CHAIN'
2944 include 'COMMON.DERIV'
2945 include 'COMMON.GEO'
2946 include 'COMMON.LOCAL'
2947 include 'COMMON.INTERACT'
2948 include 'COMMON.IOUNITS'
2949 include 'COMMON.NAMES'
2950 double precision adt,adt2,tmp
2952 tmp=1+pi_np/(2*Q_np)*0.5*d_time
2955 s12_dt=d_time/s12_np
2956 d_time_s12=d_time*0.5*s12_np
2959 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2960 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2964 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2965 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2969 if (itype(i).ne.10) then
2972 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2973 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2979 c---------------------------------------------------------------------
2980 subroutine tnp1_step2
2981 c Step 2 of the velocity Verlet algorithm: update velocities
2982 implicit real*8 (a-h,o-z)
2983 include 'DIMENSIONS'
2984 include 'COMMON.CONTROL'
2985 include 'COMMON.VAR'
2987 include 'COMMON.CHAIN'
2988 include 'COMMON.DERIV'
2989 include 'COMMON.GEO'
2990 include 'COMMON.LOCAL'
2991 include 'COMMON.INTERACT'
2992 include 'COMMON.IOUNITS'
2993 include 'COMMON.NAMES'
2995 double precision d_time_s12
2999 d_t(j,i)=d_t_new(j,i)
3006 d_time_s12=0.5d0*s12_np*d_time
3009 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
3013 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
3017 if (itype(i).ne.10) then
3020 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
3025 cd write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np
3026 pistar=pistar+(EK-0.5*(E_old+potE)
3027 & -dimen3*Rb*t_bath*log(s12_np)+H0-dimen3*Rb*t_bath)*d_time
3028 tmp=1+pistar/(2*Q_np)*0.5*d_time
3035 c-----------------------------------------------------------------
3036 subroutine tnp_respa_i_step1
3037 c Applying Nose-Poincare algorithm - step 1 to coordinates
3038 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3040 c d_t is not updated here, it is destroyed
3042 implicit real*8 (a-h,o-z)
3043 include 'DIMENSIONS'
3044 include 'COMMON.CONTROL'
3045 include 'COMMON.VAR'
3047 include 'COMMON.CHAIN'
3048 include 'COMMON.DERIV'
3049 include 'COMMON.GEO'
3050 include 'COMMON.LOCAL'
3051 include 'COMMON.INTERACT'
3052 include 'COMMON.IOUNITS'
3053 include 'COMMON.NAMES'
3054 double precision C_np,d_time_s,tmp,d_time_ss
3056 d_time_s=d_time*0.5*s_np
3057 ct2 d_time_s=d_time*0.5*s12_np
3060 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3064 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3068 if (itype(i).ne.10) then
3071 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3078 d_t(j,i)=d_t_new(j,i)
3085 C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-Csplit)
3088 pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3089 tmp=0.5*d_time*pistar/Q_np
3090 s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3092 d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3093 ct2 d_time_ss=d_time/s12_np
3094 c d_time_ss=0.5*d_time*(1.0/sold_np+1.0/s_np)
3097 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3101 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3105 if (itype(i).ne.10) then
3108 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3115 c---------------------------------------------------------------------
3117 subroutine tnp_respa_i_step2
3118 c Step 2 of the velocity Verlet algorithm: update velocities
3119 implicit real*8 (a-h,o-z)
3120 include 'DIMENSIONS'
3121 include 'COMMON.CONTROL'
3122 include 'COMMON.VAR'
3124 include 'COMMON.CHAIN'
3125 include 'COMMON.DERIV'
3126 include 'COMMON.GEO'
3127 include 'COMMON.LOCAL'
3128 include 'COMMON.INTERACT'
3129 include 'COMMON.IOUNITS'
3130 include 'COMMON.NAMES'
3132 double precision d_time_s
3134 EK=EK*(s_np/s12_np)**2
3135 HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3136 pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath
3139 cr print '(a,5f)','i_step2',EK,potE,HNose1,pi_np,E_long
3140 d_time_s=d_time*0.5*s12_np
3141 c d_time_s=d_time*0.5*s_np
3144 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3148 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3152 if (itype(i).ne.10) then
3155 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
3164 c-----------------------------------------------------------------
3165 subroutine tnp_respa_step1
3166 c Applying Nose-Poincare algorithm - step 1 to vel for RESPA
3167 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3169 c d_t is not updated here, it is destroyed
3171 implicit real*8 (a-h,o-z)
3172 include 'DIMENSIONS'
3173 include 'COMMON.CONTROL'
3174 include 'COMMON.VAR'
3176 include 'COMMON.CHAIN'
3177 include 'COMMON.DERIV'
3178 include 'COMMON.GEO'
3179 include 'COMMON.LOCAL'
3180 include 'COMMON.INTERACT'
3181 include 'COMMON.IOUNITS'
3182 include 'COMMON.NAMES'
3183 double precision C_np,d_time_s,tmp,d_time_ss
3184 double precision energia(0:n_ene)
3186 d_time_s=d_time*0.5*s_np
3189 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3193 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3197 if (itype(i).ne.10) then
3200 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3206 c C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3209 c pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3210 c tmp=0.5*d_time*pistar/Q_np
3211 c s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3212 c write(iout,*) 'tnp_respa_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3218 c-------------------------------------
3219 c test of reviewer's comment
3220 pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3221 cr print '(a,3f)','1 pi_np,s_np',pi_np,s_np,E_long
3222 c-------------------------------------
3226 c---------------------------------------------------------------------
3227 subroutine tnp_respa_step2
3228 c Step 2 of the velocity Verlet algorithm: update velocities for RESPA
3229 implicit real*8 (a-h,o-z)
3230 include 'DIMENSIONS'
3231 include 'COMMON.CONTROL'
3232 include 'COMMON.VAR'
3234 include 'COMMON.CHAIN'
3235 include 'COMMON.DERIV'
3236 include 'COMMON.GEO'
3237 include 'COMMON.LOCAL'
3238 include 'COMMON.INTERACT'
3239 include 'COMMON.IOUNITS'
3240 include 'COMMON.NAMES'
3242 double precision d_time_s
3248 ct HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3249 ct pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3250 ct & -0.5*d_time*(HNose1-H0)
3252 c-------------------------------------
3253 c test of reviewer's comment
3254 pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3255 cr print '(a,3f)','2 pi_np,s_np',pi_np,s_np,E_long
3256 c-------------------------------------
3257 d_time_s=d_time*0.5*s_np
3260 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3264 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3268 if (itype(i).ne.10) then
3271 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3280 c---------------------------------------------------------------------
3281 subroutine tnp_step1
3282 c Applying Nose-Poincare algorithm - step 1 to coordinates
3283 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3285 c d_t is not updated here, it is destroyed
3287 implicit real*8 (a-h,o-z)
3288 include 'DIMENSIONS'
3289 include 'COMMON.CONTROL'
3290 include 'COMMON.VAR'
3292 include 'COMMON.CHAIN'
3293 include 'COMMON.DERIV'
3294 include 'COMMON.GEO'
3295 include 'COMMON.LOCAL'
3296 include 'COMMON.INTERACT'
3297 include 'COMMON.IOUNITS'
3298 include 'COMMON.NAMES'
3299 double precision C_np,d_time_s,tmp,d_time_ss
3301 d_time_s=d_time*0.5*s_np
3304 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3308 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3312 if (itype(i).ne.10) then
3315 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3322 d_t(j,i)=d_t_new(j,i)
3329 C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3332 pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3333 tmp=0.5*d_time*pistar/Q_np
3334 s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3335 c write(iout,*) 'tnp_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3337 d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3340 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3344 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3348 if (itype(i).ne.10) then
3351 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3358 c-----------------------------------------------------------------
3359 subroutine tnp_step2
3360 c Step 2 of the velocity Verlet algorithm: update velocities
3361 implicit real*8 (a-h,o-z)
3362 include 'DIMENSIONS'
3363 include 'COMMON.CONTROL'
3364 include 'COMMON.VAR'
3366 include 'COMMON.CHAIN'
3367 include 'COMMON.DERIV'
3368 include 'COMMON.GEO'
3369 include 'COMMON.LOCAL'
3370 include 'COMMON.INTERACT'
3371 include 'COMMON.IOUNITS'
3372 include 'COMMON.NAMES'
3374 double precision d_time_s
3376 EK=EK*(s_np/s12_np)**2
3377 HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3378 pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3379 & -0.5*d_time*(HNose1-H0)
3381 cd write(iout,'(a,4f)') 'mmm',EK,potE,HNose1,pi_np
3382 d_time_s=d_time*0.5*s12_np
3385 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3389 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3393 if (itype(i).ne.10) then
3396 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
3406 subroutine hmc_test(itime)
3407 implicit real*8 (a-h,o-z)
3408 include 'DIMENSIONS'
3409 include 'COMMON.CONTROL'
3411 include 'COMMON.CHAIN'
3414 delta=-(potE+EK-hmc_etot)/(Rb*t_bath)
3415 if (delta .lt. -50.0d0) then
3420 xxx=ran_number(0.0d0,1.0d0)
3422 if (me.eq.king .or. .not. out1file)
3423 & write(iout,'(a8,i5,6f10.4)')
3424 & 'HMC',itime,potE+EK,potE,EK,hmc_etot,delta,xxx
3426 if (delta .le. xxx) then
3435 if (me.eq.king .or. .not. out1file)
3436 & write(iout,*) 'HMC accepting new'
3445 call chainbuild_cart
3449 d_t_old(j,i)=d_t(j,i)
3453 kinetic_T=2.0d0/(dimen3*Rb)*EK
3454 call etotal(potEcomp)
3457 if (me.eq.king .or. .not. out1file)
3458 & write(iout,'(a8,i5,3f10.4)')'HMC new',itime,potE+EK,potE,EK