2 c------------------------------------------------
3 c The driver for molecular dynamics subroutines
4 c------------------------------------------------
5 implicit real*8 (a-h,o-z)
7 include 'COMMON.CONTROL'
10 include 'COMMON.LANGEVIN'
11 include 'COMMON.CHAIN'
12 include 'COMMON.DERIV'
14 include 'COMMON.LOCAL'
15 include 'COMMON.INTERACT'
16 include 'COMMON.IOUNITS'
17 include 'COMMON.NAMES'
18 include 'COMMON.TIME1'
19 double precision cm(3),L(3),vcm(3)
20 double precision energia(0:n_ene)
24 common /gucio/ cm,energia
32 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
34 c Determine the inverse of the inertia matrix.
35 call setup_MD_matrices
38 t_MDsetup = tcpu()-tt0
40 c Entering the MD loop
42 if (lang.eq.2 .or. lang.eq.3) then
45 call sd_verlet_p_setup
47 call sd_verlet_ciccotti_setup
51 pfric0_mat(i,j,0)=pfric_mat(i,j)
52 afric0_mat(i,j,0)=afric_mat(i,j)
53 vfric0_mat(i,j,0)=vfric_mat(i,j)
54 prand0_mat(i,j,0)=prand_mat(i,j)
55 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
56 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
63 else if (lang.eq.1 .or. lang.eq.4) then
66 t_langsetup=tcpu()-tt0
70 if (lang.gt.0 .and. surfarea .and.
71 & mod(itime,reset_fricmat).eq.0) then
72 if (lang.eq.2 .or. lang.eq.3) then
75 call sd_verlet_p_setup
77 call sd_verlet_ciccotti_setup
81 pfric0_mat(i,j,0)=pfric_mat(i,j)
82 afric0_mat(i,j,0)=afric_mat(i,j)
83 vfric0_mat(i,j,0)=vfric_mat(i,j)
84 prand0_mat(i,j,0)=prand_mat(i,j)
85 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
86 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
93 else if (lang.eq.1 .or. lang.eq.4) then
96 write (iout,'(a,i10)')
97 & "Friction matrix reset based on surface area, itime",itime
99 if (reset_vel .and. tbf .and. lang.eq.0
100 & .and. mod(itime,count_reset_vel).eq.0) then
102 write(iout,'(a,f20.2)')
103 & "Velocities reset to random values, time",totT
106 d_t_old(j,i)=d_t(j,i)
110 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
114 d_t(j,0)=d_t(j,0)-vcm(j)
117 kinetic_T=2.0d0/(dimen*Rb)*EK
118 scalfac=dsqrt(T_bath/kinetic_T)
119 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
122 d_t_old(j,i)=scalfac*d_t(j,i)
128 c Time-reversible RESPA algorithm
129 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
130 call RESPA_step(itime)
132 c Variable time step algorithm.
133 call velverlet_step(itime)
136 call brown_step(itime)
138 if (mod(itime,ntwe).eq.0) call statout(itime)
139 if (mod(itime,ntwx).eq.0) then
140 write (tytul,'("time",f8.2)') totT
142 call pdbout(potE,tytul,ipdb)
147 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
148 open(irest2,file=rest2name,status='unknown')
149 write(irest2,*) totT,EK,potE,totE
151 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
154 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
161 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
163 & 'MD calculations setup:',t_MDsetup,
164 & 'Energy & gradient evaluation:',t_enegrad,
165 & 'Stochastic MD setup:',t_langsetup,
166 & 'Stochastic MD step setup:',t_sdsetup,
168 write (iout,'(/28(1h=),a25,27(1h=))')
169 & ' End of MD calculation '
172 c-------------------------------------------------------------------------------
173 subroutine brown_step(itime)
174 c------------------------------------------------
175 c Perform a single Euler integration step of Brownian dynamics
176 c------------------------------------------------
177 implicit real*8 (a-h,o-z)
179 include 'COMMON.CONTROL'
182 include 'COMMON.LANGEVIN'
183 include 'COMMON.CHAIN'
184 include 'COMMON.DERIV'
186 include 'COMMON.LOCAL'
187 include 'COMMON.INTERACT'
188 include 'COMMON.IOUNITS'
189 include 'COMMON.NAMES'
190 include 'COMMON.TIME1'
191 double precision energia(0:n_ene),zapas(MAXRES6)
192 integer ilen,rstcount
194 double precision stochforcvec(MAXRES6)
195 double precision Bmat(MAXRES6,MAXRES2),Cmat(maxres2,maxres2),
196 & Cinv(maxres2,maxres2),GBmat(MAXRES6,MAXRES2),
197 & Tmat(MAXRES6,MAXRES2),Pmat(maxres6,maxres6),Td(maxres6),
199 common /stochcalc/ stochforcvec
200 common /gucio/ cm, energia
202 logical lprn /.false./,lprn1 /.false./
204 double precision difftol /1.0d-5/
207 if (itype(i).ne.10) nbond=nbond+1
211 write (iout,*) "Generalized inverse of fricmat"
212 call matout(dimen,dimen,MAXRES6,MAXRES6,fricmat)
224 Bmat(ind+j,ind1)=dC_norm(j,i)
229 if (itype(i).ne.10) then
232 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
238 write (iout,*) "Matrix Bmat"
239 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Bmat)
245 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
250 write (iout,*) "Matrix GBmat"
251 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Gbmat)
257 Cmat(i,j)=Cmat(i,j)+Bmat(k,i)*GBmat(k,j)
262 write (iout,*) "Matrix Cmat"
263 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
265 call matinvert(nbond,MAXRES2,Cmat,Cinv)
267 write (iout,*) "Matrix Cinv"
268 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cinv)
274 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
279 write (iout,*) "Matrix Tmat"
280 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Tmat)
290 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
295 write (iout,*) "Matrix Pmat"
296 call MATOUT(dimen,dimen,MAXRES6,MAXRES6,Pmat)
303 Td(i)=Td(i)+vbl*Tmat(i,ind)
306 if (itype(k).ne.10) then
308 Td(i)=Td(i)+vbldsc0(itype(k))*Tmat(i,ind)
313 write (iout,*) "Vector Td"
315 write (iout,'(i5,f10.5)') i,Td(i)
318 call stochastic_force(stochforcvec)
320 write (iout,*) "stochforcvec"
322 write (iout,*) i,stochforcvec(i)
326 zapas(j)=-gcart(j,0)+stochforcvec(j)
328 dC_work(j)=dC_old(j,0)
334 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
335 dC_work(ind)=dC_old(j,i)
339 if (itype(i).ne.10) then
342 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
343 dC_work(ind)=dC_old(j,i+nres)
349 write (iout,*) "Initial d_t_work"
351 write (iout,*) i,d_t_work(i)
358 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
365 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
369 write (iout,*) "Final d_t_work and zapas"
371 write (iout,*) i,d_t_work(i),zapas(i)
385 dc_work(ind+j)=dc(j,i)
391 d_t(j,i+nres)=d_t_work(ind+j)
392 dc(j,i+nres)=zapas(ind+j)
393 dc_work(ind+j)=dc(j,i+nres)
399 write (iout,*) "Before correction for rotational lengthening"
400 write (iout,*) "New coordinates",
401 & " and differences between actual and standard bond lengths"
406 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
407 & i,(dC(j,i),j=1,3),xx
410 if (itype(i).ne.10) then
412 xx=vbld(i+nres)-vbldsc0(itype(i))
413 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
414 & i,(dC(j,i+nres),j=1,3),xx
418 c Second correction (rotational lengthening)
424 blen2 = scalar(dc(1,i),dc(1,i))
425 ppvec(ind)=2*vbl**2-blen2
426 diffbond=dabs(vbl-dsqrt(blen2))
427 if (diffbond.gt.diffmax) diffmax=diffbond
428 if (ppvec(ind).gt.0.0d0) then
429 ppvec(ind)=dsqrt(ppvec(ind))
434 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
438 if (itype(i).ne.10) then
440 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
441 ppvec(ind)=2*vbldsc0(itype(i))**2-blen2
442 diffbond=dabs(vbldsc0(itype(i))-dsqrt(blen2))
443 if (diffbond.gt.diffmax) diffmax=diffbond
444 if (ppvec(ind).gt.0.0d0) then
445 ppvec(ind)=dsqrt(ppvec(ind))
450 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
454 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
455 if (diffmax.lt.difftol) goto 10
459 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
465 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
476 dc_work(ind+j)=zapas(ind+j)
481 if (itype(i).ne.10) then
483 dc(j,i+nres)=zapas(ind+j)
484 dc_work(ind+j)=zapas(ind+j)
489 c Building the chain from the newly calculated coordinates
491 if (large.and. mod(itime,ntwe).eq.0) then
492 write (iout,*) "Cartesian and internal coordinates: step 1"
495 write (iout,'(a)') "Potential forces"
497 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),
498 & (-gxcart(j,i),j=1,3)
500 write (iout,'(a)') "Stochastic forces"
502 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),
503 & (stochforc(j,i+nres),j=1,3)
505 write (iout,'(a)') "Velocities"
507 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
508 & (d_t(j,i+nres),j=1,3)
512 write (iout,*) "After correction for rotational lengthening"
513 write (iout,*) "New coordinates",
514 & " and differences between actual and standard bond lengths"
519 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
520 & i,(dC(j,i),j=1,3),xx
523 if (itype(i).ne.10) then
525 xx=vbld(i+nres)-vbldsc0(itype(i))
526 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
527 & i,(dC(j,i+nres),j=1,3),xx
532 c write (iout,*) "Too many attempts at correcting the bonds"
536 c Calculate energy and forces
539 potE=energia(0)-energia(20)
542 c Calculate the kinetic and total energy and the kinetic temperature
544 t_enegrad=t_enegrad+tcpu()-tt0
546 kinetic_T=2.0d0/(dimen*Rb)*EK
549 c-------------------------------------------------------------------------------
550 subroutine velverlet_step(itime)
551 c-------------------------------------------------------------------------------
552 c Perform a single velocity Verlet step; the time step can be rescaled if
553 c increments in accelerations exceed the threshold
554 c-------------------------------------------------------------------------------
555 implicit real*8 (a-h,o-z)
557 include 'COMMON.CONTROL'
560 include 'COMMON.LANGEVIN'
561 include 'COMMON.CHAIN'
562 include 'COMMON.DERIV'
564 include 'COMMON.LOCAL'
565 include 'COMMON.INTERACT'
566 include 'COMMON.IOUNITS'
567 include 'COMMON.NAMES'
568 include 'COMMON.TIME1'
569 include 'COMMON.MUCA'
570 double precision energia(0:n_ene),vcm(3),incr(3)
571 double precision cm(3),L(3)
572 integer ilen,count,rstcount
575 integer maxcount_scale /20/
576 common /gucio/ cm, energia
577 double precision stochforcvec(MAXRES6)
578 common /stochcalc/ stochforcvec
581 double precision HNose1,HNose,HNose_nh,H,vtnp(maxres6)
582 double precision vtnp_(maxres6),vtnp_a(maxres6)
588 else if (lang.eq.2 .or. lang.eq.3) then
589 call stochastic_force(stochforcvec)
593 icount_scale=icount_scale+1
594 if (icount_scale.gt.maxcount_scale) then
596 & "ERROR: too many attempts at scaling down the time step. ",
597 & "amax=",amax,"epdrift=",epdrift,
598 & "damax=",damax,"edriftmax=",edriftmax,
602 c First step of the velocity Verlet algorithm
605 else if (lang.eq.3) then
606 call sd_verlet1_ciccotti
607 else if (lang.eq.1) then
616 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
620 d_t_old(j,i)=d_t_old(j,i)*scale_nh
626 c Build the chain from the newly calculated coordinates
628 if (rattle) call rattle1
629 if (large.and. mod(itime,ntwe).eq.0) then
630 write (iout,*) "Cartesian and internal coordinates: step 1"
633 write (iout,*) "Accelerations"
635 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
636 & (d_a(j,i+nres),j=1,3)
638 write (iout,*) "Velocities, step 1"
640 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
641 & (d_t(j,i+nres),j=1,3)
645 c Calculate energy and forces
649 potE=energia(0)-energia(20)
651 c Get the new accelerations
653 t_enegrad=t_enegrad+tcpu()-tt0
654 c Determine maximum acceleration and scale down the timestep if needed
656 call predict_edrift(epdrift)
657 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
658 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
660 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
662 itime_scal=itime_scal+ifac_time
663 c fac_time=dmin1(damax/amax,0.5d0)
664 fac_time=0.5d0**ifac_time
665 d_time=d_time*fac_time
666 if (lang.eq.2 .or. lang.eq.3) then
667 c write (iout,*) "Calling sd_verlet_setup: 1"
668 c Rescale the stochastic forces and recalculate or restore
669 c the matrices of tinker integrator
670 if (itime_scal.gt.maxflag_stoch) then
671 if (large) write (iout,'(a,i5,a)')
672 & "Calculate matrices for stochastic step;",
673 & " itime_scal ",itime_scal
675 call sd_verlet_p_setup
677 call sd_verlet_ciccotti_setup
679 write (iout,'(2a,i3,a,i3,1h.)')
680 & "Warning: cannot store matrices for stochastic",
681 & " integration because the index",itime_scal,
682 & " is greater than",maxflag_stoch
683 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
684 & " integration Langevin algorithm for better efficiency."
685 else if (flag_stoch(itime_scal)) then
686 if (large) write (iout,'(a,i5,a,l1)')
687 & "Restore matrices for stochastic step; itime_scal ",
688 & itime_scal," flag ",flag_stoch(itime_scal)
691 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
692 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
693 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
694 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
695 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
696 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
700 if (large) write (iout,'(2a,i5,a,l1)')
701 & "Calculate & store matrices for stochastic step;",
702 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
704 call sd_verlet_p_setup
706 call sd_verlet_ciccotti_setup
708 flag_stoch(ifac_time)=.true.
711 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
712 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
713 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
714 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
715 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
716 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
720 fac_time=1.0d0/dsqrt(fac_time)
722 stochforcvec(i)=fac_time*stochforcvec(i)
724 else if (lang.eq.1) then
725 c Rescale the accelerations due to stochastic forces
726 fac_time=1.0d0/dsqrt(fac_time)
728 d_as_work(i)=d_as_work(i)*fac_time
731 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
732 & "itime",itime," Timestep scaled down to ",
733 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
735 c Second step of the velocity Verlet algorithm
738 else if (lang.eq.3) then
739 call sd_verlet2_ciccotti
740 else if (lang.eq.1) then
750 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
753 d_t(j,i)=d_t(j,i)*scale_nh
758 if (rattle) call rattle2
760 if (d_time.ne.d_time0) then
762 if (lang.eq.2 .or. lang.eq.3) then
763 if (large) write (iout,'(a)')
764 & "Restore original matrices for stochastic step"
765 c write (iout,*) "Calling sd_verlet_setup: 2"
766 c Restore the matrices of tinker integrator if the time step has been restored
769 pfric_mat(i,j)=pfric0_mat(i,j,0)
770 afric_mat(i,j)=afric0_mat(i,j,0)
771 vfric_mat(i,j)=vfric0_mat(i,j,0)
772 prand_mat(i,j)=prand0_mat(i,j,0)
773 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
774 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
782 c Calculate the kinetic and the total energy and the kinetic temperature
783 if (tnp .or. tnp1) then
786 d_t_old(j,i)=d_t(j,i)
787 d_t(j,i)=d_t(j,i)/s_np
795 c write (iout,*) "step",itime," EK",EK," EK1",EK1
797 c Couple the system to Berendsen bath if needed
798 if (tbf .and. lang.eq.0) then
801 kinetic_T=2.0d0/(dimen*Rb)*EK
802 c Backup the coordinates, velocities, and accelerations
806 if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
807 d_a_old(j,i)=d_a(j,i)
810 c if (mod(itime,ntwe).eq.0 .and. large) then
811 if (mod(itime,ntwe).eq.0) then
813 if(tnp .or. tnp1) then
814 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
816 cd write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
817 cd & ,EK+potE+pi_np**2/(2*Q_np)+dimen*0.001986d0*t_bath*log(s_np)
818 write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
822 HNose1=Hnose_nh(EK,potE)
824 write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
840 if (itype(i).ne.10) then
844 vtnp(itnp)=d_t(j,inres)
849 c Transform velocities from UNRES coordinate space to cartesian and Gvec
856 vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
857 vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
859 vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
863 write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
870 c-------------------------------------------------------------------------------
871 subroutine RESPA_step(itime)
872 c-------------------------------------------------------------------------------
873 c Perform a single RESPA step.
874 c-------------------------------------------------------------------------------
875 implicit real*8 (a-h,o-z)
877 include 'COMMON.CONTROL'
880 include 'COMMON.LANGEVIN'
881 include 'COMMON.CHAIN'
882 include 'COMMON.DERIV'
884 include 'COMMON.LOCAL'
885 include 'COMMON.INTERACT'
886 include 'COMMON.IOUNITS'
887 include 'COMMON.NAMES'
888 include 'COMMON.TIME1'
889 double precision energia(0:n_ene),energia_short(0:n_ene),
890 & energia_long(0:n_ene)
891 double precision cm(3),L(3),vcm(3),incr(3)
892 integer ilen,count,rstcount
895 integer maxcount_scale /10/
897 double precision stochforcvec(MAXRES6)
898 double precision grad_tmp(3,0:maxres2)
899 common /stochcalc/ stochforcvec
902 double precision vtnp(maxres6), vtnp_(maxres6), vtnp_a(maxres6)
903 if (large.and. mod(itime,ntwe).eq.0) then
904 write (iout,*) "***************** RESPA itime",itime
905 write (iout,*) "Cartesian and internal coordinates: step 0"
908 write (iout,*) "Accelerations"
910 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
911 & (d_a(j,i+nres),j=1,3)
913 write (iout,*) "Velocities, step 0"
915 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
916 & (d_t(j,i+nres),j=1,3)
920 c Perform the initial RESPA step (increment velocities)
921 c write (iout,*) "*********************** RESPA ini"
924 creview call tnp1_respa_step1
929 if (tnh.and..not.xiresp) then
930 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
933 d_t(j,i)=d_t(j,i)*scale_nh
940 cd if(tnp .or. tnp1) then
941 cd write (iout,'(a,3f)') "EE1 NP S, pi",totT, s_np, pi_np
944 if (mod(itime,ntwe).eq.0 .and. large) then
945 write (iout,*) "Velocities, end"
947 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
948 & (d_t(j,i+nres),j=1,3)
955 cr grad_tmp(j,i)=gcart(j,i)
956 cr grad_tmp(j,i+nres)=gxcart(j,i)
961 c Compute the short-range forces
963 call etotal_short(energia_short)
964 if (tnp.or.tnp1) potE=energia_short(0)
970 if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
971 d_a_old(j,i)=d_a(j,i)
975 c Split the time step
976 d_time=d_time/ntime_split
978 c Perform the short-range RESPSA steps (velocity Verlet increments of
979 c positions and velocities using short-range forces)
980 c write (iout,*) "*********************** RESPA split"
982 do itsplit=1,ntime_split
985 else if (lang.eq.2 .or. lang.eq.3) then
986 call stochastic_force(stochforcvec)
988 c First step of the velocity Verlet algorithm
991 else if (lang.eq.3) then
992 call sd_verlet1_ciccotti
993 else if (lang.eq.1) then
996 call tnp1_respa_i_step1
997 cr if(itsplit.eq.1)then
998 cr d_time_s12=d_time0*0.5*s12_np
1000 cr d_t_half(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
1004 cr d_t_half(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
1008 cr if (itype(i).ne.10) then
1011 cr d_t_half(j,inres)=d_t_old(j,inres)
1012 cr & +d_a_old(j,inres)*d_time_s12
1018 call tnp_respa_i_step1
1020 if (tnh.and.xiresp) then
1022 call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
1025 d_t_old(j,i)=d_t_old(j,i)*scale_nh
1028 cd write(iout,*) "SSS1",itsplit,EK,scale_nh
1032 c Build the chain from the newly calculated coordinates
1033 call chainbuild_cart
1034 if (rattle) call rattle1
1035 if (large.and. mod(itime,ntwe).eq.0) then
1036 write (iout,*) "Cartesian and internal coordinates: step 1"
1039 write (iout,*) "Accelerations"
1041 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1042 & (d_a(j,i+nres),j=1,3)
1044 write (iout,*) "Velocities, step 1"
1046 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1047 & (d_t(j,i+nres),j=1,3)
1051 c Calculate energy and forces
1053 c E_long aproximation
1054 cr if (tnp .or. tnp1) then
1055 cr dtmp=0.5*d_time*(1.0/s12_np+1.0/s_np)
1058 cr E_long=E_long+d_t_new(j,i)*dtmp*grad_tmp(j,i)
1062 c-------------------------------------
1063 c test of reviewer's comment
1065 c-------------------------------------
1069 ctest call etotal_long(energia_long)
1070 ctest E_long=energia_long(0)
1074 call etotal_short(energia_short)
1076 potE=energia_short(0)
1078 c if(tnp .or. tnp1) then
1079 c write (iout,*) "kkk",E_long2,E_long
1084 c Get the new accelerations
1086 t_enegrad=t_enegrad+tcpu()-tt0
1087 c Second step of the velocity Verlet algorithm
1090 else if (lang.eq.3) then
1091 call sd_verlet2_ciccotti
1092 else if (lang.eq.1) then
1095 call tnp1_respa_i_step2
1097 call tnp_respa_i_step2
1100 if (tnh.and.xiresp) then
1102 call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
1105 d_t(j,i)=d_t(j,i)*scale_nh
1108 cd write(iout,*) "SSS2",itsplit,EK,scale_nh
1112 if (rattle) call rattle2
1113 c Backup the coordinates, velocities, and accelerations
1114 if (tnp .or. tnp1) then
1117 d_t_old(j,i)=d_t(j,i)
1118 if (tnp) d_t(j,i)=d_t(j,i)/s_np
1119 if (tnp1) d_t(j,i)=d_t(j,i)/s_np
1127 if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
1128 d_a_old(j,i)=d_a(j,i)
1132 cd if(tnp .or. tnp1) then
1134 cd HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
1135 cd H=(HNose1-H0)*s_np
1136 cd write (iout,*) "jjj",EK,potE
1137 cd write (iout,*) "iii H=",H,abs(HNose1-H0)/H0
1138 cd write (iout,'(a,3f)')
1139 cd & "III NP S, pi",totT+itsplit*d_time, s_np, pi_np
1144 c Restore the time step
1146 c Compute long-range forces
1148 call etotal_long(energia_long)
1149 E_long=energia_long(0)
1150 potE=energia_short(0)+energia_long(0)
1152 c Compute accelerations from long-range forces
1154 if (large.and. mod(itime,ntwe).eq.0) then
1155 write (iout,*) "Cartesian and internal coordinates: step 2"
1158 write (iout,*) "Accelerations"
1160 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1161 & (d_a(j,i+nres),j=1,3)
1163 write (iout,*) "Velocities, step 2"
1165 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1166 & (d_t(j,i+nres),j=1,3)
1169 c Compute the final RESPA step (increment velocities)
1170 c write (iout,*) "*********************** RESPA fin"
1172 creview call tnp1_respa_step2
1173 call tnp_respa_step2
1175 call tnp_respa_step2
1178 if (tnh.and..not.xiresp) then
1180 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
1183 d_t(j,i)=d_t(j,i)*scale_nh
1189 if (tnp .or. tnp1) then
1192 d_t(j,i)=d_t_old(j,i)/s_np
1197 c Compute the complete potential energy
1198 cc potE=energia_short(0)+energia_long(0)
1200 c Calculate the kinetic and the total energy and the kinetic temperature
1204 c Couple the system to Berendsen bath if needed
1205 if (tbf .and. lang.eq.0) then
1208 kinetic_T=2.0d0/(dimen*Rb)*EK
1209 c Backup the coordinates, velocities, and accelerations
1210 if (mod(itime,ntwe).eq.0 .and. large) then
1211 write (iout,*) "Velocities, end"
1213 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1214 & (d_t(j,i+nres),j=1,3)
1219 c if(tnp .or. tnp1) then
1220 c HNose1=Hnose(EK,s_np,energia_short(0),pi_np,Q_np,t_bath,dimen)
1221 c_new_var_csplit Csplit=H0-E_long
1222 c Csplit=H0-energia_short(0)
1227 if (mod(itime,ntwe).eq.0) then
1229 if(tnp .or. tnp1) then
1230 write (iout,'(a3,7f)') "TTT",EK,s_np,potE,pi_np,Csplit,
1231 & E_long,energia_short(0)
1232 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
1234 cd write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
1235 cd & ,EK+potE+pi_np**2/(2*Q_np)+dimen*0.001986d0*t_bath*log(s_np)
1236 write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
1237 cd write (iout,'(a,3f)') "EE2 NP S, pi",totT, s_np, pi_np
1241 HNose1=Hnose_nh(EK,potE)
1243 write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
1260 if (itype(i).ne.10) then
1264 vtnp(itnp)=d_t(j,inres)
1269 c Transform velocities from UNRES coordinate space to cartesian and Gvec
1276 vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
1277 vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
1279 vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
1283 write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
1290 c---------------------------------------------------------------------
1291 subroutine RESPA_vel
1292 c First and last RESPA step (incrementing velocities using long-range
1294 implicit real*8 (a-h,o-z)
1295 include 'DIMENSIONS'
1296 include 'COMMON.CONTROL'
1297 include 'COMMON.VAR'
1299 include 'COMMON.CHAIN'
1300 include 'COMMON.DERIV'
1301 include 'COMMON.GEO'
1302 include 'COMMON.LOCAL'
1303 include 'COMMON.INTERACT'
1304 include 'COMMON.IOUNITS'
1305 include 'COMMON.NAMES'
1307 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1311 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1315 if (itype(i).ne.10) then
1318 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1324 c-----------------------------------------------------------------
1326 c Applying velocity Verlet algorithm - step 1 to coordinates
1327 implicit real*8 (a-h,o-z)
1328 include 'DIMENSIONS'
1329 include 'COMMON.CONTROL'
1330 include 'COMMON.VAR'
1332 include 'COMMON.CHAIN'
1333 include 'COMMON.DERIV'
1334 include 'COMMON.GEO'
1335 include 'COMMON.LOCAL'
1336 include 'COMMON.INTERACT'
1337 include 'COMMON.IOUNITS'
1338 include 'COMMON.NAMES'
1339 double precision adt,adt2
1342 adt=d_a_old(j,0)*d_time
1344 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1345 d_t_new(j,0)=d_t_old(j,0)+adt2
1346 d_t(j,0)=d_t_old(j,0)+adt
1350 adt=d_a_old(j,i)*d_time
1352 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1353 d_t_new(j,i)=d_t_old(j,i)+adt2
1354 d_t(j,i)=d_t_old(j,i)+adt
1358 if (itype(i).ne.10) then
1361 adt=d_a_old(j,inres)*d_time
1363 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1364 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1365 d_t(j,inres)=d_t_old(j,inres)+adt
1371 c---------------------------------------------------------------------
1373 c Step 2 of the velocity Verlet algorithm: update velocities
1374 implicit real*8 (a-h,o-z)
1375 include 'DIMENSIONS'
1376 include 'COMMON.CONTROL'
1377 include 'COMMON.VAR'
1379 include 'COMMON.CHAIN'
1380 include 'COMMON.DERIV'
1381 include 'COMMON.GEO'
1382 include 'COMMON.LOCAL'
1383 include 'COMMON.INTERACT'
1384 include 'COMMON.IOUNITS'
1385 include 'COMMON.NAMES'
1387 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1391 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1395 if (itype(i).ne.10) then
1398 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1405 subroutine sddir_precalc
1406 c Applying velocity Verlet algorithm - step 1 to coordinates
1407 implicit real*8 (a-h,o-z)
1408 include 'DIMENSIONS'
1409 include 'COMMON.CONTROL'
1410 include 'COMMON.VAR'
1412 include 'COMMON.LANGEVIN'
1413 include 'COMMON.CHAIN'
1414 include 'COMMON.DERIV'
1415 include 'COMMON.GEO'
1416 include 'COMMON.LOCAL'
1417 include 'COMMON.INTERACT'
1418 include 'COMMON.IOUNITS'
1419 include 'COMMON.NAMES'
1420 double precision stochforcvec(MAXRES6)
1421 common /stochcalc/ stochforcvec
1423 c Compute friction and stochastic forces
1426 call stochastic_force(stochforcvec)
1428 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1429 c forces (d_as_work)
1435 d_af_work(i)=d_af_work(i)+Ginv(i,j)*fric_work(j)
1436 d_as_work(i)=d_as_work(i)+Ginv(i,j)*stochforcvec(j)
1441 c---------------------------------------------------------------------
1442 subroutine sddir_verlet1
1443 c Applying velocity Verlet algorithm - step 1 to velocities
1444 implicit real*8 (a-h,o-z)
1445 include 'DIMENSIONS'
1446 include 'COMMON.CONTROL'
1447 include 'COMMON.VAR'
1449 include 'COMMON.LANGEVIN'
1450 include 'COMMON.CHAIN'
1451 include 'COMMON.DERIV'
1452 include 'COMMON.GEO'
1453 include 'COMMON.LOCAL'
1454 include 'COMMON.INTERACT'
1455 include 'COMMON.IOUNITS'
1456 include 'COMMON.NAMES'
1457 c Revised 3/31/05 AL: correlation between random contributions to
1458 c position and velocity increments included.
1459 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1460 double precision adt,adt2
1462 c Add the contribution from BOTH friction and stochastic force to the
1463 c coordinates, but ONLY the contribution from the friction forces to velocities
1466 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1467 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1468 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1469 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1470 d_t(j,0)=d_t_old(j,0)+adt
1475 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1476 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1477 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1478 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1479 d_t(j,i)=d_t_old(j,i)+adt
1484 if (itype(i).ne.10) then
1487 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1488 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1489 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1490 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1491 d_t(j,inres)=d_t_old(j,inres)+adt
1498 c---------------------------------------------------------------------
1499 subroutine sddir_verlet2
1500 c Calculating the adjusted velocities for accelerations
1501 implicit real*8 (a-h,o-z)
1502 include 'DIMENSIONS'
1503 include 'COMMON.CONTROL'
1504 include 'COMMON.VAR'
1506 include 'COMMON.LANGEVIN'
1507 include 'COMMON.CHAIN'
1508 include 'COMMON.DERIV'
1509 include 'COMMON.GEO'
1510 include 'COMMON.LOCAL'
1511 include 'COMMON.INTERACT'
1512 include 'COMMON.IOUNITS'
1513 include 'COMMON.NAMES'
1514 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1515 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1516 c Revised 3/31/05 AL: correlation between random contributions to
1517 c position and velocity increments included.
1518 c The correlation coefficients are calculated at low-friction limit.
1519 c Also, friction forces are now not calculated with new velocities.
1521 c call friction_force
1522 call stochastic_force(stochforcvec)
1524 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1525 c forces (d_as_work)
1528 c d_af_work(i)=0.0d0
1531 c d_af_work(i)=d_af_work(i)+Ginv(i,j)*fric_work(j)
1532 d_as_work1(i)=d_as_work1(i)+Ginv(i,j)*stochforcvec(j)
1539 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1540 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1545 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1546 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1551 if (itype(i).ne.10) then
1554 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1555 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1556 & +cos60*d_as_work1(ind+j))*d_time
1563 c---------------------------------------------------------------------
1564 subroutine max_accel
1566 c Find the maximum difference in the accelerations of the the sites
1567 c at the beginning and the end of the time step.
1569 implicit real*8 (a-h,o-z)
1570 include 'DIMENSIONS'
1571 include 'COMMON.CONTROL'
1572 include 'COMMON.VAR'
1574 include 'COMMON.CHAIN'
1575 include 'COMMON.DERIV'
1576 include 'COMMON.GEO'
1577 include 'COMMON.LOCAL'
1578 include 'COMMON.INTERACT'
1579 include 'COMMON.IOUNITS'
1580 double precision aux(3),accel(3)
1582 aux(j)=d_a(j,0)-d_a_old(j,0)
1589 accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1590 if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1597 if (itype(i).ne.10) then
1599 accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1603 if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1606 aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1611 c---------------------------------------------------------------------
1612 subroutine predict_edrift(epdrift)
1614 c Predict the drift of the potential energy
1616 implicit real*8 (a-h,o-z)
1617 include 'DIMENSIONS'
1618 include 'COMMON.CONTROL'
1619 include 'COMMON.VAR'
1621 include 'COMMON.CHAIN'
1622 include 'COMMON.DERIV'
1623 include 'COMMON.GEO'
1624 include 'COMMON.LOCAL'
1625 include 'COMMON.INTERACT'
1626 include 'COMMON.IOUNITS'
1627 include 'COMMON.MUCA'
1628 double precision epdrift,epdriftij
1629 c Drift of the potential energy
1635 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1636 if (lmuca) epdriftij=epdriftij*factor
1637 c write (iout,*) "back",i,j,epdriftij
1638 if (epdriftij.gt.epdrift) epdrift=epdriftij
1642 if (itype(i).ne.10) then
1645 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1646 if (lmuca) epdriftij=epdriftij*factor
1647 c write (iout,*) "side",i,j,epdriftij
1648 if (epdriftij.gt.epdrift) epdrift=epdriftij
1652 epdrift=0.5d0*epdrift*d_time*d_time
1653 c write (iout,*) "epdrift",epdrift
1656 c-----------------------------------------------------------------------
1657 subroutine verlet_bath
1659 c Coupling to the thermostat by using the Berendsen algorithm
1661 implicit real*8 (a-h,o-z)
1662 include 'DIMENSIONS'
1663 include 'COMMON.CONTROL'
1664 include 'COMMON.VAR'
1666 include 'COMMON.CHAIN'
1667 include 'COMMON.DERIV'
1668 include 'COMMON.GEO'
1669 include 'COMMON.LOCAL'
1670 include 'COMMON.INTERACT'
1671 include 'COMMON.IOUNITS'
1672 include 'COMMON.NAMES'
1673 double precision T_half,fact
1675 T_half=2.0d0/(dimen*Rb)*EK
1676 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1677 c write(iout,*) "T_half", T_half
1678 c write(iout,*) "EK", EK
1679 c write(iout,*) "fact", fact
1681 d_t(j,0)=fact*d_t(j,0)
1685 d_t(j,i)=fact*d_t(j,i)
1689 if (itype(i).ne.10) then
1692 d_t(j,inres)=fact*d_t(j,inres)
1698 c---------------------------------------------------------
1700 c Set up the initial conditions of a MD simulation
1701 implicit real*8 (a-h,o-z)
1702 include 'DIMENSIONS'
1705 include 'COMMON.INFO'
1706 include 'COMMON.SETUP'
1710 include 'COMMON.CONTROL'
1711 include 'COMMON.VAR'
1713 include 'COMMON.LANGEVIN'
1714 include 'COMMON.CHAIN'
1715 include 'COMMON.DERIV'
1716 include 'COMMON.GEO'
1717 include 'COMMON.LOCAL'
1718 include 'COMMON.INTERACT'
1719 include 'COMMON.IOUNITS'
1720 include 'COMMON.NAMES'
1721 real*8 energia(0:n_ene),energia_long(0:n_ene),
1722 & energia_short(0:n_ene),vcm(3),incr(3),E_short
1723 double precision cm(3),L(3),xv,sigv,lowb,highb
1730 write(iout,*) "d_time", d_time
1731 c Compute the standard deviations of stochastic forces for Langevin dynamics
1732 c if the friction coefficients do not depend on surface area
1733 if (lang.gt.0 .and. .not.surfarea) then
1735 stdforcp(i)=stdfp*dsqrt(gamp)
1738 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
1741 c Open the pdb file for snapshotshots
1743 if (nprocs.eq.1) then
1746 npos = dlog10(dfloat(nprocs-1))+1
1748 if (npos.lt.3) npos=3
1749 write (liczba,'(i1)') npos
1750 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1752 write (liczba,form) myrank
1755 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1758 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1763 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1765 cartname=prefix(:ilen(prefix))//"_MD.cx"
1769 write (qstr,'(256(1h ))')
1775 write (iout,*) "Frag",qinfrag(i),wfrag(i),iq,iw
1776 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1784 write (iout,*) "Pair",i,qinpair(i),wpair(i),iq,iw
1785 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1789 pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1790 cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1791 statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1795 write(iout,*) "Initial state will be read from file ",
1796 & rest2name(:ilen(rest2name))
1799 c Generate initial velocities
1800 write(iout,*) "Initial velocities randomly generated"
1804 c rest2name = prefix(:ilen(prefix))//'.rst'
1805 write(iout,*) "Initial backbone velocities"
1807 write(iout,*) (d_t(j,i),j=1,3)
1809 write(iout,*) "Initial side-chain velocities"
1811 write(iout,*) (d_t(j,i+nres),j=1,3)
1813 c Zeroing the total angular momentum of the system
1814 write(iout,*) "Calling the zero-angular
1815 & momentum subroutine"
1817 c Getting the potential energy and forces and velocities and accelerations
1819 c write (iout,*) "velocity of the center of the mass:"
1820 c write (iout,*) (vcm(j),j=1,3)
1822 d_t(j,0)=d_t(j,0)-vcm(j)
1824 c Removing the velocity of the center of mass
1826 write (iout,*) "vcm right after adjustment:"
1827 write (iout,*) (vcm(j),j=1,3)
1831 call chainbuild_cart
1834 call verlet_bath(EK)
1836 kinetic_T=2.0d0/(dimen*Rb)*EK
1840 call etotal(energia)
1843 if(tnp .or. tnp1) then
1846 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
1848 write(iout,*) 'H0= ',H0
1852 HNose1=Hnose_nh(EK,potE)
1854 write (iout,*) 'H0= ',H0
1861 if (amax*d_time .gt. dvmax) d_time=d_time*dvmax/amax
1862 write(iout,*) "Potential energy"
1863 write(iout,*) (energia(i),i=0,n_ene)
1864 potE=energia(0)-energia(20)
1868 write (iout,*) "Initial:",
1869 & " Kinetic energy",EK," potential energy",potE,
1870 & " total energy",totE," maximum acceleration ",
1873 write (iout,*) "Initial velocities"
1875 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1876 & (d_t(j,i+nres),j=1,3)
1878 write (iout,*) "Initial accelerations"
1880 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1881 & (d_a(j,i+nres),j=1,3)
1887 d_t_old(j,i)=d_t(j,i)
1888 d_a_old(j,i)=d_a(j,i)
1890 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1894 call etotal_long(energia_long)
1895 E_long=energia_long(0)
1896 if(tnp .or. tnp1) then
1897 call etotal_short(energia_short)
1898 E_short=energia_short(0)
1899 HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen)
1902 c_new_var_csplit Csplit=H0-E_long
1903 c Csplit = H0-energia_short(0)
1904 write(iout,*) 'Csplit= ',Csplit
1908 c call etotal_short(energia_short)
1909 c write (iout,*) "energia_long",energia_long(0),
1910 c & " energia_short",energia_short(0),
1911 c & " total",energia_long(0)+energia_short(0)
1915 c-----------------------------------------------------------
1916 subroutine random_vel
1917 implicit real*8 (a-h,o-z)
1918 include 'DIMENSIONS'
1919 include 'COMMON.CONTROL'
1920 include 'COMMON.VAR'
1922 include 'COMMON.LANGEVIN'
1923 include 'COMMON.CHAIN'
1924 include 'COMMON.DERIV'
1925 include 'COMMON.GEO'
1926 include 'COMMON.LOCAL'
1927 include 'COMMON.INTERACT'
1928 include 'COMMON.IOUNITS'
1929 include 'COMMON.NAMES'
1930 include 'COMMON.TIME1'
1931 double precision xv,sigv,lowb,highb
1932 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1933 c First generate velocities in the eigenspace of the G matrix
1934 c write (iout,*) "Calling random_vel"
1937 sigv=dsqrt((Rb*t_bath)/geigen(i))
1940 d_t_work_new(i)=anorm_distr(xv,sigv,lowb,highb)
1944 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(i)**2
1946 c Transform velocities to UNRES coordinate space
1950 d_t_work(i)=d_t_work(i)+Gvec(i,j)*d_t_work_new(j)
1953 c Transfer to the d_t vector
1955 d_t(j,0)=d_t_work(j)
1961 d_t(j,i)=d_t_work(ind)
1965 if (itype(i).ne.10) then
1968 d_t(j,i+nres)=d_t_work(ind)
1973 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1974 c & 2.0d0/(dimen*Rb)*EK,2.0d0/(dimen*Rb)*EK1
1977 c-----------------------------------------------------------
1978 subroutine sd_verlet_p_setup
1979 c Sets up the parameters of stochastic Verlet algorithm
1980 implicit real*8 (a-h,o-z)
1981 include 'DIMENSIONS'
1982 include 'COMMON.CONTROL'
1983 include 'COMMON.VAR'
1985 include 'COMMON.LANGEVIN'
1986 include 'COMMON.CHAIN'
1987 include 'COMMON.DERIV'
1988 include 'COMMON.GEO'
1989 include 'COMMON.LOCAL'
1990 include 'COMMON.INTERACT'
1991 include 'COMMON.IOUNITS'
1992 include 'COMMON.NAMES'
1993 include 'COMMON.TIME1'
1994 double precision emgdt(MAXRES6),
1995 & pterm,vterm,rho,rhoc,vsig,
1996 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1997 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1998 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1999 logical lprn /.false./
2000 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2001 double precision ktm
2004 c AL 8/17/04 Code adapted from tinker
2006 c Get the frictional and random terms for stochastic dynamics in the
2007 c eigenspace of mass-scaled UNRES friction matrix
2010 gdt = fricgam(i) * d_time
2012 c Stochastic dynamics reduces to simple MD for zero friction
2014 if (gdt .le. zero) then
2015 pfric_vec(i) = 1.0d0
2016 vfric_vec(i) = d_time
2017 afric_vec(i) = 0.5d0 * d_time * d_time
2018 prand_vec(i) = 0.0d0
2019 vrand_vec1(i) = 0.0d0
2020 vrand_vec2(i) = 0.0d0
2022 c Analytical expressions when friction coefficient is large
2025 if (gdt .ge. gdt_radius) then
2028 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2029 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2030 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2031 vterm = 1.0d0 - egdt**2
2032 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2034 c Use series expansions when friction coefficient is small
2045 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2046 & - gdt5/120.0d0 + gdt6/720.0d0
2047 & - gdt7/5040.0d0 + gdt8/40320.0d0
2048 & - gdt9/362880.0d0) / fricgam(i)**2
2049 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2050 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2051 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2052 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2053 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2054 & + 127.0d0*gdt9/90720.0d0
2055 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2056 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2057 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2058 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2059 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2060 & - 17.0d0*gdt2/1280.0d0
2061 & + 17.0d0*gdt3/6144.0d0
2062 & + 40967.0d0*gdt4/34406400.0d0
2063 & - 57203.0d0*gdt5/275251200.0d0
2064 & - 1429487.0d0*gdt6/13212057600.0d0)
2067 c Compute the scaling factors of random terms for the nonzero friction case
2069 ktm = 0.5d0*d_time/fricgam(i)
2070 psig = dsqrt(ktm*pterm) / fricgam(i)
2071 vsig = dsqrt(ktm*vterm)
2072 rhoc = dsqrt(1.0d0 - rho*rho)
2074 vrand_vec1(i) = vsig * rho
2075 vrand_vec2(i) = vsig * rhoc
2080 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2083 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2084 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2088 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2090 call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
2091 call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
2092 call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
2093 call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
2094 call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec1,vrand_mat1)
2095 call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
2096 c call eigtransf1(dimen,maxres6,mt3mt2,pfric_vec,pfric_mat)
2097 c call eigtransf1(dimen,maxres6,mt3mt2,vfric_vec,vfric_mat)
2098 c call eigtransf1(dimen,maxres6,mt3mt2,afric_vec,afric_mat)
2099 c call eigtransf1(dimen,maxres6,mt3mt1,prand_vec,prand_mat)
2100 c call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec1,vrand_mat1)
2101 c call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec2,vrand_mat2)
2102 t_sdsetup=t_sdsetup+tcpu()-tt0
2105 c-------------------------------------------------------------
2106 subroutine eigtransf1(n,ndim,ab,d,c)
2109 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2115 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2121 c-------------------------------------------------------------
2122 subroutine eigtransf(n,ndim,a,b,d,c)
2125 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2131 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2137 c-------------------------------------------------------------
2138 subroutine sd_verlet1
2139 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2140 implicit real*8 (a-h,o-z)
2141 include 'DIMENSIONS'
2142 include 'COMMON.CONTROL'
2143 include 'COMMON.VAR'
2145 include 'COMMON.LANGEVIN'
2146 include 'COMMON.CHAIN'
2147 include 'COMMON.DERIV'
2148 include 'COMMON.GEO'
2149 include 'COMMON.LOCAL'
2150 include 'COMMON.INTERACT'
2151 include 'COMMON.IOUNITS'
2152 include 'COMMON.NAMES'
2153 double precision stochforcvec(MAXRES6)
2154 common /stochcalc/ stochforcvec
2155 logical lprn /.false./
2157 c write (iout,*) "dc_old"
2159 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2160 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2163 dc_work(j)=dc_old(j,0)
2164 d_t_work(j)=d_t_old(j,0)
2165 d_a_work(j)=d_a_old(j,0)
2170 dc_work(ind+j)=dc_old(j,i)
2171 d_t_work(ind+j)=d_t_old(j,i)
2172 d_a_work(ind+j)=d_a_old(j,i)
2177 if (itype(i).ne.10) then
2179 dc_work(ind+j)=dc_old(j,i+nres)
2180 d_t_work(ind+j)=d_t_old(j,i+nres)
2181 d_a_work(ind+j)=d_a_old(j,i+nres)
2189 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2193 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2194 & vfric_mat(i,j),afric_mat(i,j),
2195 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2203 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2204 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2205 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2206 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2208 d_t_work_new(i)=ddt1+0.5d0*ddt2
2209 d_t_work(i)=ddt1+ddt2
2213 d_t(j,0)=d_t_work(j)
2218 dc(j,i)=dc_work(ind+j)
2219 d_t(j,i)=d_t_work(ind+j)
2224 if (itype(i).ne.10) then
2227 dc(j,inres)=dc_work(ind+j)
2228 d_t(j,inres)=d_t_work(ind+j)
2235 c--------------------------------------------------------------------------
2236 subroutine sd_verlet2
2237 c Calculating the adjusted velocities for accelerations
2238 implicit real*8 (a-h,o-z)
2239 include 'DIMENSIONS'
2240 include 'COMMON.CONTROL'
2241 include 'COMMON.VAR'
2243 include 'COMMON.LANGEVIN'
2244 include 'COMMON.CHAIN'
2245 include 'COMMON.DERIV'
2246 include 'COMMON.GEO'
2247 include 'COMMON.LOCAL'
2248 include 'COMMON.INTERACT'
2249 include 'COMMON.IOUNITS'
2250 include 'COMMON.NAMES'
2251 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2252 common /stochcalc/ stochforcvec
2254 c Compute the stochastic forces which contribute to velocity change
2256 call stochastic_force(stochforcvecV)
2262 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2263 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2264 & vrand_mat2(i,j)*stochforcvecV(j)
2266 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2269 d_t(j,0)=d_t_work(j)
2274 d_t(j,i)=d_t_work(ind+j)
2279 if (itype(i).ne.10) then
2282 d_t(j,inres)=d_t_work(ind+j)
2289 c-----------------------------------------------------------
2290 subroutine sd_verlet_ciccotti_setup
2291 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2293 implicit real*8 (a-h,o-z)
2294 include 'DIMENSIONS'
2295 include 'COMMON.CONTROL'
2296 include 'COMMON.VAR'
2298 include 'COMMON.LANGEVIN'
2299 include 'COMMON.CHAIN'
2300 include 'COMMON.DERIV'
2301 include 'COMMON.GEO'
2302 include 'COMMON.LOCAL'
2303 include 'COMMON.INTERACT'
2304 include 'COMMON.IOUNITS'
2305 include 'COMMON.NAMES'
2306 include 'COMMON.TIME1'
2307 double precision emgdt(MAXRES6),
2308 & pterm,vterm,rho,rhoc,vsig,
2309 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2310 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2311 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2312 logical lprn /.false./
2313 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2314 double precision ktm
2317 c AL 8/17/04 Code adapted from tinker
2319 c Get the frictional and random terms for stochastic dynamics in the
2320 c eigenspace of mass-scaled UNRES friction matrix
2323 write (iout,*) "i",i," fricgam",fricgam(i)
2324 gdt = fricgam(i) * d_time
2326 c Stochastic dynamics reduces to simple MD for zero friction
2328 if (gdt .le. zero) then
2329 pfric_vec(i) = 1.0d0
2330 vfric_vec(i) = d_time
2331 afric_vec(i) = 0.5d0*d_time*d_time
2332 prand_vec(i) = afric_vec(i)
2333 vrand_vec2(i) = vfric_vec(i)
2335 c Analytical expressions when friction coefficient is large
2340 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2341 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2342 prand_vec(i) = afric_vec(i)
2343 vrand_vec2(i) = vfric_vec(i)
2345 c Compute the scaling factors of random terms for the nonzero friction case
2347 c ktm = 0.5d0*d_time/fricgam(i)
2348 c psig = dsqrt(ktm*pterm) / fricgam(i)
2349 c vsig = dsqrt(ktm*vterm)
2350 c prand_vec(i) = psig*afric_vec(i)
2351 c vrand_vec2(i) = vsig*vfric_vec(i)
2356 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2359 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2360 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2364 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2366 call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
2367 call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
2368 call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
2369 call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
2370 call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
2371 t_sdsetup=t_sdsetup+tcpu()-tt0
2374 c-------------------------------------------------------------
2375 subroutine sd_verlet1_ciccotti
2376 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2377 implicit real*8 (a-h,o-z)
2378 include 'DIMENSIONS'
2379 include 'COMMON.CONTROL'
2380 include 'COMMON.VAR'
2382 include 'COMMON.LANGEVIN'
2383 include 'COMMON.CHAIN'
2384 include 'COMMON.DERIV'
2385 include 'COMMON.GEO'
2386 include 'COMMON.LOCAL'
2387 include 'COMMON.INTERACT'
2388 include 'COMMON.IOUNITS'
2389 include 'COMMON.NAMES'
2390 double precision stochforcvec(MAXRES6)
2391 common /stochcalc/ stochforcvec
2392 logical lprn /.false./
2394 c write (iout,*) "dc_old"
2396 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2397 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2400 dc_work(j)=dc_old(j,0)
2401 d_t_work(j)=d_t_old(j,0)
2402 d_a_work(j)=d_a_old(j,0)
2407 dc_work(ind+j)=dc_old(j,i)
2408 d_t_work(ind+j)=d_t_old(j,i)
2409 d_a_work(ind+j)=d_a_old(j,i)
2414 if (itype(i).ne.10) then
2416 dc_work(ind+j)=dc_old(j,i+nres)
2417 d_t_work(ind+j)=d_t_old(j,i+nres)
2418 d_a_work(ind+j)=d_a_old(j,i+nres)
2426 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2430 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2431 & vfric_mat(i,j),afric_mat(i,j),
2432 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2440 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2441 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2442 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2443 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2445 d_t_work_new(i)=ddt1+0.5d0*ddt2
2446 d_t_work(i)=ddt1+ddt2
2450 d_t(j,0)=d_t_work(j)
2455 dc(j,i)=dc_work(ind+j)
2456 d_t(j,i)=d_t_work(ind+j)
2461 if (itype(i).ne.10) then
2464 dc(j,inres)=dc_work(ind+j)
2465 d_t(j,inres)=d_t_work(ind+j)
2472 c--------------------------------------------------------------------------
2473 subroutine sd_verlet2_ciccotti
2474 c Calculating the adjusted velocities for accelerations
2475 implicit real*8 (a-h,o-z)
2476 include 'DIMENSIONS'
2477 include 'COMMON.CONTROL'
2478 include 'COMMON.VAR'
2480 include 'COMMON.LANGEVIN'
2481 include 'COMMON.CHAIN'
2482 include 'COMMON.DERIV'
2483 include 'COMMON.GEO'
2484 include 'COMMON.LOCAL'
2485 include 'COMMON.INTERACT'
2486 include 'COMMON.IOUNITS'
2487 include 'COMMON.NAMES'
2488 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2489 common /stochcalc/ stochforcvec
2491 c Compute the stochastic forces which contribute to velocity change
2493 call stochastic_force(stochforcvecV)
2499 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2500 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2501 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2503 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2506 d_t(j,0)=d_t_work(j)
2511 d_t(j,i)=d_t_work(ind+j)
2516 if (itype(i).ne.10) then
2519 d_t(j,inres)=d_t_work(ind+j)
2526 c------------------------------------------------------
2527 double precision function HNose(ek,s,e,pi,Q,t_bath,dimen)
2529 double precision ek,s,e,pi,Q,t_bath,Rb
2532 HNose=ek+e+pi**2/(2*Q)+dimen*Rb*t_bath*log(s)
2533 c print '(6f15.5,i5,a2,2f15.5)',ek,s,e,pi,Q,t_bath,dimen,"--",
2534 c & pi**2/(2*Q),dimen*Rb*t_bath*log(s)
2537 c-----------------------------------------------------------------
2538 double precision function HNose_nh(eki,e)
2539 implicit real*8 (a-h,o-z)
2540 include 'DIMENSIONS'
2542 HNose_nh=eki+e+dimen*Rb*t_bath*xlogs(1)+qmass(1)*vlogs(1)**2/2
2544 HNose_nh=HNose_nh+qmass(i)*vlogs(i)**2/2+Rb*t_bath*xlogs(i)
2546 c write(4,'(5e15.5)')
2547 c & vlogs(1),xlogs(1),HNose,eki,e
2550 c-------------------------------------------------------
2552 subroutine tnp1_step1
2553 c Applying Nose-Poincare algorithm - step 1 to coordinates
2554 c JPSJ 70 75 (2001) S. Nose
2556 c d_t is not updated here
2558 implicit real*8 (a-h,o-z)
2559 include 'DIMENSIONS'
2560 include 'COMMON.CONTROL'
2561 include 'COMMON.VAR'
2563 include 'COMMON.CHAIN'
2564 include 'COMMON.DERIV'
2565 include 'COMMON.GEO'
2566 include 'COMMON.LOCAL'
2567 include 'COMMON.INTERACT'
2568 include 'COMMON.IOUNITS'
2569 include 'COMMON.NAMES'
2570 double precision adt,adt2,tmp
2572 tmp=1+pi_np/(2*Q_np)*0.5*d_time
2575 s12_dt=d_time/s12_np
2576 d_time_s12=d_time*0.5*s12_np
2579 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2580 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2584 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2585 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2589 if (itype(i).ne.10) then
2592 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2593 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2599 c---------------------------------------------------------------------
2600 subroutine tnp1_step2
2601 c Step 2 of the velocity Verlet algorithm: update velocities
2602 implicit real*8 (a-h,o-z)
2603 include 'DIMENSIONS'
2604 include 'COMMON.CONTROL'
2605 include 'COMMON.VAR'
2607 include 'COMMON.CHAIN'
2608 include 'COMMON.DERIV'
2609 include 'COMMON.GEO'
2610 include 'COMMON.LOCAL'
2611 include 'COMMON.INTERACT'
2612 include 'COMMON.IOUNITS'
2613 include 'COMMON.NAMES'
2615 double precision d_time_s12
2619 d_t(j,i)=d_t_new(j,i)
2626 d_time_s12=0.5d0*s12_np*d_time
2629 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
2633 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
2637 if (itype(i).ne.10) then
2640 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
2645 cd write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np
2646 pistar=pistar+(EK-0.5*(E_old+potE)
2647 & -dimen*Rb*t_bath*log(s12_np)+H0-dimen*Rb*t_bath)*d_time
2648 tmp=1+pistar/(2*Q_np)*0.5*d_time
2654 c---------------------------------------------------------------------
2655 subroutine tnp_step1
2656 c Applying Nose-Poincare algorithm - step 1 to coordinates
2657 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
2659 c d_t is not updated here, it is destroyed
2661 implicit real*8 (a-h,o-z)
2662 include 'DIMENSIONS'
2663 include 'COMMON.CONTROL'
2664 include 'COMMON.VAR'
2666 include 'COMMON.CHAIN'
2667 include 'COMMON.DERIV'
2668 include 'COMMON.GEO'
2669 include 'COMMON.LOCAL'
2670 include 'COMMON.INTERACT'
2671 include 'COMMON.IOUNITS'
2672 include 'COMMON.NAMES'
2673 double precision C_np,d_time_s,tmp,d_time_ss
2675 d_time_s=d_time*0.5*s_np
2678 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
2682 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
2686 if (itype(i).ne.10) then
2689 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
2696 d_t(j,i)=d_t_new(j,i)
2703 C_np=0.5*d_time*(dimen*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
2706 pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
2707 tmp=0.5*d_time*pistar/Q_np
2708 s12_np=s_np*(1.0+tmp)/(1.0-tmp)
2709 c write(iout,*) 'tnp_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
2711 d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
2714 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
2718 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
2722 if (itype(i).ne.10) then
2725 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
2732 c-----------------------------------------------------------------
2733 subroutine tnp_step2
2734 c Step 2 of the velocity Verlet algorithm: update velocities
2735 implicit real*8 (a-h,o-z)
2736 include 'DIMENSIONS'
2737 include 'COMMON.CONTROL'
2738 include 'COMMON.VAR'
2740 include 'COMMON.CHAIN'
2741 include 'COMMON.DERIV'
2742 include 'COMMON.GEO'
2743 include 'COMMON.LOCAL'
2744 include 'COMMON.INTERACT'
2745 include 'COMMON.IOUNITS'
2746 include 'COMMON.NAMES'
2748 double precision d_time_s
2750 EK=EK*(s_np/s12_np)**2
2751 HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen)
2752 pi_np=pistar+0.5*d_time*(2*EK-dimen*Rb*t_bath)
2753 & -0.5*d_time*(HNose1-H0)
2755 cd write(iout,'(a,4f)') 'mmm',EK,potE,HNose1,pi_np
2756 d_time_s=d_time*0.5*s12_np
2759 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
2763 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
2767 if (itype(i).ne.10) then
2770 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
2779 c-----------------------------------------------------------------
2780 subroutine tnp1_respa_i_step1
2781 c Applying Nose-Poincare algorithm - step 1 to coordinates
2782 c JPSJ 70 75 (2001) S. Nose
2784 c d_t is not updated here
2786 implicit real*8 (a-h,o-z)
2787 include 'DIMENSIONS'
2788 include 'COMMON.CONTROL'
2789 include 'COMMON.VAR'
2791 include 'COMMON.CHAIN'
2792 include 'COMMON.DERIV'
2793 include 'COMMON.GEO'
2794 include 'COMMON.LOCAL'
2795 include 'COMMON.INTERACT'
2796 include 'COMMON.IOUNITS'
2797 include 'COMMON.NAMES'
2798 double precision adt,adt2,tmp
2800 tmp=1+pi_np/(2*Q_np)*0.5*d_time
2803 s12_dt=d_time/s12_np
2804 d_time_s12=d_time*0.5*s12_np
2807 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2808 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2812 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2813 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2817 if (itype(i).ne.10) then
2820 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2821 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2827 c---------------------------------------------------------------------
2828 c-----------------------------------------------------------------
2829 subroutine tnp1_respa_step1_
2830 c Applying Nose-Poincare algorithm - step 1 to vel for RESPA
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_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s12
2860 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s12
2864 if (itype(i).ne.10) then
2867 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s12
2873 c---------------------------------------------------------------------
2874 subroutine tnp1_respa_i_step2
2875 c Step 2 of the velocity Verlet algorithm: update velocities
2876 implicit real*8 (a-h,o-z)
2877 include 'DIMENSIONS'
2878 include 'COMMON.CONTROL'
2879 include 'COMMON.VAR'
2881 include 'COMMON.CHAIN'
2882 include 'COMMON.DERIV'
2883 include 'COMMON.GEO'
2884 include 'COMMON.LOCAL'
2885 include 'COMMON.INTERACT'
2886 include 'COMMON.IOUNITS'
2887 include 'COMMON.NAMES'
2889 double precision d_time_s12
2893 d_t(j,i)=d_t_new(j,i)
2900 d_time_s12=0.5d0*s12_np*d_time
2903 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
2907 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
2911 if (itype(i).ne.10) then
2914 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
2919 pistar=pistar+(EK-0.5*(E_old+potE)
2920 & -dimen*Rb*t_bath*log(s12_np)+Csplit-dimen*Rb*t_bath)*d_time
2921 tmp=1+pistar/(2*Q_np)*0.5*d_time
2927 c---------------------------------------------------------------------
2928 subroutine tnp1_respa_step2_
2929 c Step 2 of the velocity Verlet algorithm: update velocities for RESPA
2930 implicit real*8 (a-h,o-z)
2931 include 'DIMENSIONS'
2932 include 'COMMON.CONTROL'
2933 include 'COMMON.VAR'
2935 include 'COMMON.CHAIN'
2936 include 'COMMON.DERIV'
2937 include 'COMMON.GEO'
2938 include 'COMMON.LOCAL'
2939 include 'COMMON.INTERACT'
2940 include 'COMMON.IOUNITS'
2941 include 'COMMON.NAMES'
2943 double precision d_time_s12
2947 d_t(j,i)=d_t_half(j,i)
2954 d_time_s12=0.5d0*s12_np*d_time
2957 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s12
2961 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s12
2965 if (itype(i).ne.10) then
2968 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s12
2973 cd write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np
2974 pistar=pistar+(EK-0.5*(E_old+potE)
2975 & -dimen*Rb*t_bath*log(s12_np)+H0-dimen*Rb*t_bath)*d_time
2976 tmp=1+pistar/(2*Q_np)*0.5*d_time
2983 c-----------------------------------------------------------------
2984 subroutine tnp_respa_i_step1
2985 c Applying Nose-Poincare algorithm - step 1 to coordinates
2986 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
2988 c d_t is not updated here, it is destroyed
2990 implicit real*8 (a-h,o-z)
2991 include 'DIMENSIONS'
2992 include 'COMMON.CONTROL'
2993 include 'COMMON.VAR'
2995 include 'COMMON.CHAIN'
2996 include 'COMMON.DERIV'
2997 include 'COMMON.GEO'
2998 include 'COMMON.LOCAL'
2999 include 'COMMON.INTERACT'
3000 include 'COMMON.IOUNITS'
3001 include 'COMMON.NAMES'
3002 double precision C_np,d_time_s,tmp,d_time_ss
3004 d_time_s=d_time*0.5*s_np
3005 ct2 d_time_s=d_time*0.5*s12_np
3008 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3012 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3016 if (itype(i).ne.10) then
3019 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3026 d_t(j,i)=d_t_new(j,i)
3033 C_np=0.5*d_time*(dimen*Rb*t_bath*(1.0+log(s_np))-EK+potE-Csplit)
3036 pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3037 tmp=0.5*d_time*pistar/Q_np
3038 s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3040 d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3041 ct2 d_time_ss=d_time/s12_np
3042 c d_time_ss=0.5*d_time*(1.0/sold_np+1.0/s_np)
3045 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3049 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3053 if (itype(i).ne.10) then
3056 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3063 c-----------------------------------------------------------------
3064 subroutine tnp_respa_step1
3065 c Applying Nose-Poincare algorithm - step 1 to vel for RESPA
3066 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3068 c d_t is not updated here, it is destroyed
3070 implicit real*8 (a-h,o-z)
3071 include 'DIMENSIONS'
3072 include 'COMMON.CONTROL'
3073 include 'COMMON.VAR'
3075 include 'COMMON.CHAIN'
3076 include 'COMMON.DERIV'
3077 include 'COMMON.GEO'
3078 include 'COMMON.LOCAL'
3079 include 'COMMON.INTERACT'
3080 include 'COMMON.IOUNITS'
3081 include 'COMMON.NAMES'
3082 double precision C_np,d_time_s,tmp,d_time_ss
3083 double precision energia(0:n_ene)
3085 d_time_s=d_time*0.5*s_np
3088 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3092 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3096 if (itype(i).ne.10) then
3099 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3105 c C_np=0.5*d_time*(dimen*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3108 c pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3109 c tmp=0.5*d_time*pistar/Q_np
3110 c s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3111 c write(iout,*) 'tnp_respa_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3117 c-------------------------------------
3118 c test of reviewer's comment
3119 pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3120 cr print '(a,3f)','1 pi_np,s_np',pi_np,s_np,E_long
3121 c-------------------------------------
3125 c---------------------------------------------------------------------
3127 subroutine tnp_respa_i_step2
3128 c Step 2 of the velocity Verlet algorithm: update velocities
3129 implicit real*8 (a-h,o-z)
3130 include 'DIMENSIONS'
3131 include 'COMMON.CONTROL'
3132 include 'COMMON.VAR'
3134 include 'COMMON.CHAIN'
3135 include 'COMMON.DERIV'
3136 include 'COMMON.GEO'
3137 include 'COMMON.LOCAL'
3138 include 'COMMON.INTERACT'
3139 include 'COMMON.IOUNITS'
3140 include 'COMMON.NAMES'
3142 double precision d_time_s
3144 EK=EK*(s_np/s12_np)**2
3145 HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen)
3146 pi_np=pistar+0.5*d_time*(2*EK-dimen*Rb*t_bath
3149 cr print '(a,5f)','i_step2',EK,potE,HNose1,pi_np,E_long
3150 d_time_s=d_time*0.5*s12_np
3151 c d_time_s=d_time*0.5*s_np
3154 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3158 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3162 if (itype(i).ne.10) then
3165 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
3174 c---------------------------------------------------------------------
3175 subroutine tnp_respa_step2
3176 c Step 2 of the velocity Verlet algorithm: update velocities for RESPA
3177 implicit real*8 (a-h,o-z)
3178 include 'DIMENSIONS'
3179 include 'COMMON.CONTROL'
3180 include 'COMMON.VAR'
3182 include 'COMMON.CHAIN'
3183 include 'COMMON.DERIV'
3184 include 'COMMON.GEO'
3185 include 'COMMON.LOCAL'
3186 include 'COMMON.INTERACT'
3187 include 'COMMON.IOUNITS'
3188 include 'COMMON.NAMES'
3190 double precision d_time_s
3196 ct HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen)
3197 ct pi_np=pistar+0.5*d_time*(2*EK-dimen*Rb*t_bath)
3198 ct & -0.5*d_time*(HNose1-H0)
3200 c-------------------------------------
3201 c test of reviewer's comment
3202 pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3203 cr print '(a,3f)','2 pi_np,s_np',pi_np,s_np,E_long
3204 c-------------------------------------
3205 d_time_s=d_time*0.5*s_np
3208 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3212 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3216 if (itype(i).ne.10) then
3219 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3228 c-----------------------------------------------------------------
3229 SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8)
3230 implicit real*8 (a-h,o-z)
3231 include 'DIMENSIONS'
3233 double precision akin,gnkt,dt,aa,gkt,scale
3234 double precision wdti(maxyosh),wdti2(maxyosh),
3235 & wdti4(maxyosh),wdti8(maxyosh)
3236 integer i,iresn,iyosh,inos,nnos1
3245 C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE
3246 C INTEGRATION FROM t=0 TO t=DT/2
3247 C GET THE TOTAL KINETIC ENERGY
3249 c CALL GETKINP(MASS,VX,VY,VZ,AKIN)
3251 GLOGS(1) = (AKIN - GNKT)/QMASS(1)
3252 C START THE MULTIPLE TIME STEP PROCEDURE
3255 C UPDATE THE THERMOSTAT VELOCITIES
3256 VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
3258 AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) )
3259 VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA
3260 & + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA
3262 C UPDATE THE PARTICLE VELOCITIES
3263 AA = EXP(-WDTI2(IYOSH)*VLOGS(1) )
3266 GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1)
3267 C UPDATE THE THERMOSTAT POSITIONS
3269 XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH)
3271 C UPDATE THE THERMOSTAT VELOCITIES
3273 AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) )
3274 VLOGS(INOS) = VLOGS(INOS)*AA*AA
3275 & + WDTI4(IYOSH)*GLOGS(INOS)*AA
3276 GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS)
3277 & -GKT)/QMASS(INOS+1)
3279 VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
3282 C UPDATE THE PARTICLE VELOCITIES
3283 c outside of this subroutine
3285 c VX(I) = VX(I)*SCALE
3286 c VY(I) = VY(I)*SCALE
3287 c VZ(I) = VZ(I)*SCALE