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 include 'COMMON.HAIRPIN'
29 double precision cm(3),L(3),vcm(3)
31 double precision v_work(maxres6),v_transf(maxres6)
41 if (ilen(tmpdir).gt.0)
42 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
43 & //liczba(:ilen(liczba))//'.rst')
45 if (ilen(tmpdir).gt.0)
46 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
53 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
59 c Determine the inverse of the inertia matrix.
60 call setup_MD_matrices
64 t_MDsetup = MPI_Wtime()-tt0
66 t_MDsetup = tcpu()-tt0
69 c Entering the MD loop
75 if (lang.eq.2 .or. lang.eq.3) then
79 call sd_verlet_p_setup
81 call sd_verlet_ciccotti_setup
85 pfric0_mat(i,j,0)=pfric_mat(i,j)
86 afric0_mat(i,j,0)=afric_mat(i,j)
87 vfric0_mat(i,j,0)=vfric_mat(i,j)
88 prand0_mat(i,j,0)=prand_mat(i,j)
89 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
90 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
99 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
101 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
105 else if (lang.eq.1 .or. lang.eq.4) then
109 t_langsetup=MPI_Wtime()-tt0
112 t_langsetup=tcpu()-tt0
115 do itime=1,n_timestep
117 if (large.and. mod(itime,ntwe).eq.0)
118 & write (iout,*) "itime",itime
120 if (lang.gt.0 .and. surfarea .and.
121 & mod(itime,reset_fricmat).eq.0) then
122 if (lang.eq.2 .or. lang.eq.3) then
126 call sd_verlet_p_setup
128 call sd_verlet_ciccotti_setup
132 pfric0_mat(i,j,0)=pfric_mat(i,j)
133 afric0_mat(i,j,0)=afric_mat(i,j)
134 vfric0_mat(i,j,0)=vfric_mat(i,j)
135 prand0_mat(i,j,0)=prand_mat(i,j)
136 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
137 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
142 flag_stoch(i)=.false.
145 else if (lang.eq.1 .or. lang.eq.4) then
148 write (iout,'(a,i10)')
149 & "Friction matrix reset based on surface area, itime",itime
151 if (reset_vel .and. tbf .and. lang.eq.0
152 & .and. mod(itime,count_reset_vel).eq.0) then
154 write(iout,'(a,f20.2)')
155 & "Velocities reset to random values, time",totT
158 d_t_old(j,i)=d_t(j,i)
162 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
166 d_t(j,0)=d_t(j,0)-vcm(j)
169 kinetic_T=2.0d0/(dimen3*Rb)*EK
170 scalfac=dsqrt(T_bath/kinetic_T)
171 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
174 d_t_old(j,i)=scalfac*d_t(j,i)
180 c Time-reversible RESPA algorithm
181 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
182 call RESPA_step(itime)
184 c Variable time step algorithm.
185 call velverlet_step(itime)
189 call brown_step(itime)
191 print *,"Brown dynamics not here!"
193 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
199 if (mod(itime,ntwe).eq.0) then
201 C call enerprint(potEcomp)
202 C print *,itime,'AFM',Eafmforc,etot
216 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
219 v_work(ind)=d_t(j,i+nres)
224 write (66,'(80f10.5)')
225 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
229 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
231 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
233 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
236 if (mod(itime,ntwx).eq.0) then
237 write(iout,*) 'time=',itime
238 C call check_ecartint
240 write (tytul,'("time",f8.2)') totT
242 call hairpin(.true.,nharp,iharp)
243 call secondary2(.true.)
244 call pdbout(potE,tytul,ipdb)
249 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
250 open(irest2,file=rest2name,status='unknown')
251 write(irest2,*) totT,EK,potE,totE,t_bath
253 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
256 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
267 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
269 & 'MD calculations setup:',t_MDsetup,
270 & 'Energy & gradient evaluation:',t_enegrad,
271 & 'Stochastic MD setup:',t_langsetup,
272 & 'Stochastic MD step setup:',t_sdsetup,
274 write (iout,'(/28(1h=),a25,27(1h=))')
275 & ' End of MD calculation '
277 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
279 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
280 & " time_fricmatmult",time_fricmatmult," time_fsample ",
285 c-------------------------------------------------------------------------------
286 subroutine velverlet_step(itime)
287 c-------------------------------------------------------------------------------
288 c Perform a single velocity Verlet step; the time step can be rescaled if
289 c increments in accelerations exceed the threshold
290 c-------------------------------------------------------------------------------
291 implicit real*8 (a-h,o-z)
295 integer ierror,ierrcode
297 include 'COMMON.SETUP'
298 include 'COMMON.CONTROL'
302 include 'COMMON.LANGEVIN'
304 include 'COMMON.LANGEVIN.lang0'
306 include 'COMMON.CHAIN'
307 include 'COMMON.DERIV'
309 include 'COMMON.LOCAL'
310 include 'COMMON.INTERACT'
311 include 'COMMON.IOUNITS'
312 include 'COMMON.NAMES'
313 include 'COMMON.TIME1'
314 include 'COMMON.MUCA'
315 double precision vcm(3),incr(3)
316 double precision cm(3),L(3)
317 integer ilen,count,rstcount
320 integer maxcount_scale /20/
322 double precision stochforcvec(MAXRES6)
323 common /stochcalc/ stochforcvec
326 double precision HNose1,HNose,HNose_nh,H,vtnp(maxres6)
327 double precision vtnp_(maxres6),vtnp_a(maxres6)
333 else if (lang.eq.2 .or. lang.eq.3) then
335 call stochastic_force(stochforcvec)
338 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
340 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
347 icount_scale=icount_scale+1
348 if (icount_scale.gt.maxcount_scale) then
350 & "ERROR: too many attempts at scaling down the time step. ",
351 & "amax=",amax,"epdrift=",epdrift,
352 & "damax=",damax,"edriftmax=",edriftmax,
356 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
360 c First step of the velocity Verlet algorithm
365 else if (lang.eq.3) then
367 call sd_verlet1_ciccotti
369 else if (lang.eq.1) then
371 C else if (tnp1) then
377 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
380 d_t_old(j,i)=d_t_old(j,i)*scale_nh
386 c Build the chain from the newly calculated coordinates
388 if (rattle) call rattle1
390 if (large.and. mod(itime,ntwe).eq.0) then
391 write (iout,*) "Cartesian and internal coordinates: step 1"
396 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
397 & (dc(j,i+nres),j=1,3)
399 write (iout,*) "Accelerations"
401 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
402 & (d_a(j,i+nres),j=1,3)
404 write (iout,*) "Velocities, step 1"
406 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
407 & (d_t(j,i+nres),j=1,3)
416 c Calculate energy and forces
418 call etotal(potEcomp)
419 if (large.and. mod(itime,ntwe).eq.0)
420 & call enerprint(potEcomp)
423 t_etotal=t_etotal+MPI_Wtime()-tt0
425 t_etotal=t_etotal+tcpu()-tt0
429 potE=potEcomp(0)-potEcomp(20)
431 c Get the new accelerations
434 t_enegrad=t_enegrad+MPI_Wtime()-tt0
436 t_enegrad=t_enegrad+tcpu()-tt0
438 c Determine maximum acceleration and scale down the timestep if needed
440 amax=amax/(itime_scal+1)**2
441 call predict_edrift(epdrift)
442 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
443 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
445 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
447 itime_scal=itime_scal+ifac_time
448 c fac_time=dmin1(damax/amax,0.5d0)
449 fac_time=0.5d0**ifac_time
450 d_time=d_time*fac_time
451 if (lang.eq.2 .or. lang.eq.3) then
453 c write (iout,*) "Calling sd_verlet_setup: 1"
454 c Rescale the stochastic forces and recalculate or restore
455 c the matrices of tinker integrator
456 if (itime_scal.gt.maxflag_stoch) then
457 if (large) write (iout,'(a,i5,a)')
458 & "Calculate matrices for stochastic step;",
459 & " itime_scal ",itime_scal
461 call sd_verlet_p_setup
463 call sd_verlet_ciccotti_setup
465 write (iout,'(2a,i3,a,i3,1h.)')
466 & "Warning: cannot store matrices for stochastic",
467 & " integration because the index",itime_scal,
468 & " is greater than",maxflag_stoch
469 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
470 & " integration Langevin algorithm for better efficiency."
471 else if (flag_stoch(itime_scal)) then
472 if (large) write (iout,'(a,i5,a,l1)')
473 & "Restore matrices for stochastic step; itime_scal ",
474 & itime_scal," flag ",flag_stoch(itime_scal)
477 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
478 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
479 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
480 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
481 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
482 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
486 if (large) write (iout,'(2a,i5,a,l1)')
487 & "Calculate & store matrices for stochastic step;",
488 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
490 call sd_verlet_p_setup
492 call sd_verlet_ciccotti_setup
494 flag_stoch(ifac_time)=.true.
497 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
498 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
499 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
500 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
501 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
502 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
506 fac_time=1.0d0/dsqrt(fac_time)
508 stochforcvec(i)=fac_time*stochforcvec(i)
511 else if (lang.eq.1) then
512 c Rescale the accelerations due to stochastic forces
513 fac_time=1.0d0/dsqrt(fac_time)
515 d_as_work(i)=d_as_work(i)*fac_time
518 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
519 & "itime",itime," Timestep scaled down to ",
520 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
522 c Second step of the velocity Verlet algorithm
527 else if (lang.eq.3) then
529 call sd_verlet2_ciccotti
531 else if (lang.eq.1) then
533 C> else if (tnp1) then
535 C> else if (tnp) then
541 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
544 d_t(j,i)=d_t(j,i)*scale_nh
549 if (rattle) call rattle2
552 C print *,totTafm,"TU?"
553 if (d_time.ne.d_time0) then
556 if (lang.eq.2 .or. lang.eq.3) then
557 if (large) write (iout,'(a)')
558 & "Restore original matrices for stochastic step"
559 c write (iout,*) "Calling sd_verlet_setup: 2"
560 c Restore the matrices of tinker integrator if the time step has been restored
563 pfric_mat(i,j)=pfric0_mat(i,j,0)
564 afric_mat(i,j)=afric0_mat(i,j,0)
565 vfric_mat(i,j)=vfric0_mat(i,j,0)
566 prand_mat(i,j)=prand0_mat(i,j,0)
567 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
568 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
577 if (tnp .or. tnp1) then
580 d_t_old(j,i)=d_t(j,i)
581 d_t(j,i)=d_t(j,i)/s_np
586 c Calculate the kinetic and the total energy and the kinetic temperature
591 c write (iout,*) "step",itime," EK",EK," EK1",EK1
593 c Couple the system to Berendsen bath if needed
594 if (tbf .and. lang.eq.0) then
597 kinetic_T=2.0d0/(dimen3*Rb)*EK
598 c Backup the coordinates, velocities, and accelerations
602 if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
603 d_a_old(j,i)=d_a(j,i)
607 if (mod(itime,ntwe).eq.0 .and. large) then
608 if(tnp .or. tnp1) then
609 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
611 cd write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
612 cd & ,EK+potE+pi_np**2/(2*Q_np)+dimen3*0.001986d0*t_bath*log(s_np)
613 cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
618 HNose1=Hnose_nh(EK,potE)
621 cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
637 if (itype(i).ne.10) then
641 vtnp(itnp)=d_t(j,inres)
646 c Transform velocities from UNRES coordinate space to cartesian and Gvec
653 vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
654 vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
656 vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
660 write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
663 write (iout,*) "Velocities, step 2"
665 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
666 & (d_t(j,i+nres),j=1,3)
673 c-------------------------------------------------------------------------------
674 subroutine RESPA_step(itime)
675 c-------------------------------------------------------------------------------
676 c Perform a single RESPA step.
677 c-------------------------------------------------------------------------------
678 implicit real*8 (a-h,o-z)
682 integer IERROR,ERRCODE
684 include 'COMMON.SETUP'
685 include 'COMMON.CONTROL'
689 include 'COMMON.LANGEVIN'
691 include 'COMMON.LANGEVIN.lang0'
693 include 'COMMON.CHAIN'
694 include 'COMMON.DERIV'
696 include 'COMMON.LOCAL'
697 include 'COMMON.INTERACT'
698 include 'COMMON.IOUNITS'
699 include 'COMMON.NAMES'
700 include 'COMMON.TIME1'
701 double precision energia_short(0:n_ene),
702 & energia_long(0:n_ene)
703 double precision cm(3),L(3),vcm(3),incr(3)
704 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
705 & d_a_old0(3,0:maxres2)
706 logical PRINT_AMTS_MSG /.false./
707 integer ilen,count,rstcount
710 integer maxcount_scale /10/
712 double precision stochforcvec(MAXRES6)
713 common /stochcalc/ stochforcvec
715 double precision vtnp(maxres6), vtnp_(maxres6), vtnp_a(maxres6)
717 common /cipiszcze/ itt
720 if (large.and. mod(itime,ntwe).eq.0) then
721 write (iout,*) "***************** RESPA itime",itime
722 write (iout,*) "Cartesian and internal coordinates: step 0"
724 call pdbout(0.0d0,"cipiszcze",iout)
726 write (iout,*) "Accelerations from long-range forces"
728 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
729 & (d_a(j,i+nres),j=1,3)
731 write (iout,*) "Velocities, step 0"
733 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
734 & (d_t(j,i+nres),j=1,3)
739 c Perform the initial RESPA step (increment velocities)
740 c write (iout,*) "*********************** RESPA ini"
746 if (tnh.and..not.xiresp) then
747 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
750 d_t(j,i)=d_t(j,i)*scale_nh
757 cd if(tnp .or. tnp1) then
758 cd write (iout,'(a,3f)') "EE1 NP S, pi",totT, s_np, pi_np
762 if (mod(itime,ntwe).eq.0 .and. large) then
763 write (iout,*) "Velocities, end"
765 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
766 & (d_t(j,i+nres),j=1,3)
770 c Compute the short-range forces
776 C 7/2/2009 commented out
777 if (tnp.or.tnp1) potE=energia_short(0)
779 c call etotal_short(energia_short)
782 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
785 d_a(j,i)=d_a_short(j,i)
789 if (large.and. mod(itime,ntwe).eq.0) then
790 write (iout,*) "energia_short",energia_short(0)
791 write (iout,*) "Accelerations from short-range forces"
793 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
794 & (d_a(j,i+nres),j=1,3)
799 t_enegrad=t_enegrad+MPI_Wtime()-tt0
801 t_enegrad=t_enegrad+tcpu()-tt0
806 C d_t_old(j,i)=d_t(j,i)
807 if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
808 d_a_old(j,i)=d_a(j,i)
811 c 6/30/08 A-MTS: attempt at increasing the split number
814 dc_old0(j,i)=dc_old(j,i)
815 d_t_old0(j,i)=d_t_old(j,i)
816 d_a_old0(j,i)=d_a_old(j,i)
819 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
820 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
827 c write (iout,*) "itime",itime," ntime_split",ntime_split
828 c Split the time step
829 d_time=d_time0/ntime_split
830 c Perform the short-range RESPA steps (velocity Verlet increments of
831 c positions and velocities using short-range forces)
832 c write (iout,*) "*********************** RESPA split"
833 do itsplit=1,ntime_split
836 else if (lang.eq.2 .or. lang.eq.3) then
838 call stochastic_force(stochforcvec)
841 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
843 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
848 c First step of the velocity Verlet algorithm
853 else if (lang.eq.3) then
855 call sd_verlet1_ciccotti
857 else if (lang.eq.1) then
860 call tnp1_respa_i_step1
862 call tnp_respa_i_step1
864 if (tnh.and.xiresp) then
866 call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
869 d_t_old(j,i)=d_t_old(j,i)*scale_nh
872 cd write(iout,*) "SSS1",itsplit,EK,scale_nh
876 c Build the chain from the newly calculated coordinates
878 if (rattle) call rattle1
880 if (large.and. mod(itime,ntwe).eq.0) then
881 write (iout,*) "***** ITSPLIT",itsplit
882 write (iout,*) "Cartesian and internal coordinates: step 1"
883 call pdbout(0.0d0,"cipiszcze",iout)
886 write (iout,*) "Velocities, step 1"
888 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
889 & (d_t(j,i+nres),j=1,3)
898 c Calculate energy and forces
900 call etotal_short(energia_short)
901 if (large.and. mod(itime,ntwe).eq.0)
902 & call enerprint(energia_short)
904 potE=energia_short(0)
908 t_eshort=t_eshort+MPI_Wtime()-tt0
910 t_eshort=t_eshort+tcpu()-tt0
914 c Get the new accelerations
916 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
919 d_a_short(j,i)=d_a(j,i)
923 if (large.and. mod(itime,ntwe).eq.0) then
924 write (iout,*)"energia_short",energia_short(0)
925 write (iout,*) "Accelerations from short-range forces"
927 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
928 & (d_a(j,i+nres),j=1,3)
933 c Determine maximum acceleration and scale down the timestep if needed
935 amax=amax/ntime_split**2
936 call predict_edrift(epdrift)
937 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
938 & write (iout,*) "amax",amax," damax",damax,
939 & " epdrift",epdrift," epdriftmax",epdriftmax
940 c Exit loop and try with increased split number if the change of
941 c acceleration is too big
942 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
943 if (ntime_split.lt.maxtime_split) then
945 ntime_split=ntime_split*2
948 dc_old(j,i)=dc_old0(j,i)
949 d_t_old(j,i)=d_t_old0(j,i)
950 d_a_old(j,i)=d_a_old0(j,i)
953 if (PRINT_AMTS_MSG) then
954 write (iout,*) "acceleration/energy drift too large",amax,
955 & epdrift," split increased to ",ntime_split," itime",itime,
961 & "Uh-hu. Bumpy landscape. Maximum splitting number",
963 & " already reached!!! Trying to carry on!"
967 t_enegrad=t_enegrad+MPI_Wtime()-tt0
969 t_enegrad=t_enegrad+tcpu()-tt0
971 c Second step of the velocity Verlet algorithm
976 else if (lang.eq.3) then
978 call sd_verlet2_ciccotti
980 else if (lang.eq.1) then
983 call tnp1_respa_i_step2
985 call tnp_respa_i_step2
988 if (tnh.and.xiresp) then
990 call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
993 d_t(j,i)=d_t(j,i)*scale_nh
996 cd write(iout,*) "SSS2",itsplit,EK,scale_nh
999 if (rattle) call rattle2
1000 if (tnp .or. tnp1) then
1003 d_t_old(j,i)=d_t(j,i)
1004 if (tnp) d_t(j,i)=d_t(j,i)/s_np
1005 if (tnp1) d_t(j,i)=d_t(j,i)/s_np
1011 c Backup the coordinates, velocities, and accelerations
1015 if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
1016 d_a_old(j,i)=d_a(j,i)
1023 c Restore the time step
1025 c Compute long-range forces
1032 call etotal_long(energia_long)
1033 if (large.and. mod(itime,ntwe).eq.0)
1034 & call enerprint(energia_long)
1035 E_long=energia_long(0)
1036 potE=energia_short(0)+energia_long(0)
1039 t_elong=t_elong+MPI_Wtime()-tt0
1041 t_elong=t_elong+tcpu()-tt0
1047 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1049 t_enegrad=t_enegrad+tcpu()-tt0
1051 c Compute accelerations from long-range forces
1053 if (large.and. mod(itime,ntwe).eq.0) then
1054 write (iout,*) "energia_long",energia_long(0)
1055 write (iout,*) "Cartesian and internal coordinates: step 2"
1057 call pdbout(0.0d0,"cipiszcze",iout)
1059 write (iout,*) "Accelerations from long-range forces"
1061 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1062 & (d_a(j,i+nres),j=1,3)
1064 write (iout,*) "Velocities, step 2"
1066 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1067 & (d_t(j,i+nres),j=1,3)
1071 c Compute the final RESPA step (increment velocities)
1072 c write (iout,*) "*********************** RESPA fin"
1075 call tnp_respa_step2
1077 call tnp_respa_step2
1080 if (tnh.and..not.xiresp) then
1082 call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
1085 d_t(j,i)=d_t(j,i)*scale_nh
1091 if (tnp .or. tnp1) then
1094 d_t(j,i)=d_t_old(j,i)/s_np
1099 c Compute the complete potential energy
1101 potEcomp(i)=energia_short(i)+energia_long(i)
1103 potE=potEcomp(0)-potEcomp(20)
1104 c potE=energia_short(0)+energia_long(0)
1107 c Calculate the kinetic and the total energy and the kinetic temperature
1110 c Couple the system to Berendsen bath if needed
1111 if (tbf .and. lang.eq.0) then
1114 kinetic_T=2.0d0/(dimen3*Rb)*EK
1115 c Backup the coordinates, velocities, and accelerations
1117 if (mod(itime,ntwe).eq.0 .and. large) then
1118 write (iout,*) "Velocities, end"
1120 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1121 & (d_t(j,i+nres),j=1,3)
1124 if (mod(itime,ntwe).eq.0) then
1126 if(tnp .or. tnp1) then
1128 write (iout,'(a3,7f20.10)') "TTT",EK,s_np,potE,pi_np,Csplit,
1129 & E_long,energia_short(0)
1131 write (iout,'(a3,7f20.10)') "TTT",EK,s_np,potE,pi_np,Csplit,
1132 & E_long,energia_short(0)
1134 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
1136 cd write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
1137 cd & ,EK+potE+pi_np**2/(2*Q_np)+dimen3*0.001986d0*t_bath*log(s_np)
1138 cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
1140 cd write (iout,'(a,3f)') "EE2 NP S, pi",totT, s_np, pi_np
1144 HNose1=Hnose_nh(EK,potE)
1146 cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
1164 if (itype(i).ne.10) then
1168 vtnp(itnp)=d_t(j,inres)
1172 c Transform velocities from UNRES coordinate space to cartesian and Gvec
1179 vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
1180 vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
1182 vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
1186 write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
1194 c---------------------------------------------------------------------
1195 subroutine RESPA_vel
1196 c First and last RESPA step (incrementing velocities using long-range
1198 implicit real*8 (a-h,o-z)
1199 include 'DIMENSIONS'
1200 include 'COMMON.CONTROL'
1201 include 'COMMON.VAR'
1203 include 'COMMON.CHAIN'
1204 include 'COMMON.DERIV'
1205 include 'COMMON.GEO'
1206 include 'COMMON.LOCAL'
1207 include 'COMMON.INTERACT'
1208 include 'COMMON.IOUNITS'
1209 include 'COMMON.NAMES'
1211 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1215 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1219 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1222 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1228 c-----------------------------------------------------------------
1230 c Applying velocity Verlet algorithm - step 1 to coordinates
1231 implicit real*8 (a-h,o-z)
1232 include 'DIMENSIONS'
1233 include 'COMMON.CONTROL'
1234 include 'COMMON.VAR'
1236 include 'COMMON.CHAIN'
1237 include 'COMMON.DERIV'
1238 include 'COMMON.GEO'
1239 include 'COMMON.LOCAL'
1240 include 'COMMON.INTERACT'
1241 include 'COMMON.IOUNITS'
1242 include 'COMMON.NAMES'
1243 double precision adt,adt2
1246 write (iout,*) "VELVERLET1 START: DC"
1248 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1249 & (dc(j,i+nres),j=1,3)
1253 adt=d_a_old(j,0)*d_time
1255 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1256 d_t_new(j,0)=d_t_old(j,0)+adt2
1257 d_t(j,0)=d_t_old(j,0)+adt
1263 adt=d_a_old(j,i)*d_time
1265 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1266 d_t_new(j,i)=d_t_old(j,i)+adt2
1267 d_t(j,i)=d_t_old(j,i)+adt
1272 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1275 adt=d_a_old(j,inres)*d_time
1277 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1278 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1279 d_t(j,inres)=d_t_old(j,inres)+adt
1284 write (iout,*) "VELVERLET1 END: DC"
1286 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1287 & (dc(j,i+nres),j=1,3)
1292 c---------------------------------------------------------------------
1294 c Step 2 of the velocity Verlet algorithm: update velocities
1295 implicit real*8 (a-h,o-z)
1296 include 'DIMENSIONS'
1297 include 'COMMON.CONTROL'
1298 include 'COMMON.VAR'
1300 include 'COMMON.CHAIN'
1301 include 'COMMON.DERIV'
1302 include 'COMMON.GEO'
1303 include 'COMMON.LOCAL'
1304 include 'COMMON.INTERACT'
1305 include 'COMMON.IOUNITS'
1306 include 'COMMON.NAMES'
1308 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1312 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1316 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1319 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1325 c-----------------------------------------------------------------
1326 subroutine sddir_precalc
1327 c Applying velocity Verlet algorithm - step 1 to coordinates
1328 implicit real*8 (a-h,o-z)
1329 include 'DIMENSIONS'
1333 include 'COMMON.CONTROL'
1334 include 'COMMON.VAR'
1337 include 'COMMON.LANGEVIN'
1339 include 'COMMON.LANGEVIN.lang0'
1341 include 'COMMON.CHAIN'
1342 include 'COMMON.DERIV'
1343 include 'COMMON.GEO'
1344 include 'COMMON.LOCAL'
1345 include 'COMMON.INTERACT'
1346 include 'COMMON.IOUNITS'
1347 include 'COMMON.NAMES'
1348 include 'COMMON.TIME1'
1349 double precision stochforcvec(MAXRES6)
1350 common /stochcalc/ stochforcvec
1352 c Compute friction and stochastic forces
1361 time_fric=time_fric+MPI_Wtime()-time00
1364 time_fric=time_fric+tcpu()-time00
1367 call stochastic_force(stochforcvec)
1369 time_stoch=time_stoch+MPI_Wtime()-time00
1371 time_stoch=time_stoch+tcpu()-time00
1374 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1375 c forces (d_as_work)
1377 call ginv_mult(fric_work, d_af_work)
1378 call ginv_mult(stochforcvec, d_as_work)
1381 c---------------------------------------------------------------------
1382 subroutine sddir_verlet1
1383 c Applying velocity Verlet algorithm - step 1 to velocities
1384 implicit real*8 (a-h,o-z)
1385 include 'DIMENSIONS'
1386 include 'COMMON.CONTROL'
1387 include 'COMMON.VAR'
1390 include 'COMMON.LANGEVIN'
1392 include 'COMMON.LANGEVIN.lang0'
1394 include 'COMMON.CHAIN'
1395 include 'COMMON.DERIV'
1396 include 'COMMON.GEO'
1397 include 'COMMON.LOCAL'
1398 include 'COMMON.INTERACT'
1399 include 'COMMON.IOUNITS'
1400 include 'COMMON.NAMES'
1401 c Revised 3/31/05 AL: correlation between random contributions to
1402 c position and velocity increments included.
1403 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1404 double precision adt,adt2
1406 c Add the contribution from BOTH friction and stochastic force to the
1407 c coordinates, but ONLY the contribution from the friction forces to velocities
1410 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1411 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1412 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1413 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1414 d_t(j,0)=d_t_old(j,0)+adt
1419 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1420 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1421 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1422 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1423 d_t(j,i)=d_t_old(j,i)+adt
1428 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1431 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1432 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1433 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1434 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1435 d_t(j,inres)=d_t_old(j,inres)+adt
1442 c---------------------------------------------------------------------
1443 subroutine sddir_verlet2
1444 c Calculating the adjusted velocities for accelerations
1445 implicit real*8 (a-h,o-z)
1446 include 'DIMENSIONS'
1447 include 'COMMON.CONTROL'
1448 include 'COMMON.VAR'
1451 include 'COMMON.LANGEVIN'
1453 include 'COMMON.LANGEVIN.lang0'
1455 include 'COMMON.CHAIN'
1456 include 'COMMON.DERIV'
1457 include 'COMMON.GEO'
1458 include 'COMMON.LOCAL'
1459 include 'COMMON.INTERACT'
1460 include 'COMMON.IOUNITS'
1461 include 'COMMON.NAMES'
1462 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1463 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1464 c Revised 3/31/05 AL: correlation between random contributions to
1465 c position and velocity increments included.
1466 c The correlation coefficients are calculated at low-friction limit.
1467 c Also, friction forces are now not calculated with new velocities.
1469 c call friction_force
1470 call stochastic_force(stochforcvec)
1472 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1473 c forces (d_as_work)
1475 call ginv_mult(stochforcvec, d_as_work1)
1481 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1482 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1487 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1488 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1493 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1496 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1497 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1498 & +cos60*d_as_work1(ind+j))*d_time
1505 c---------------------------------------------------------------------
1506 subroutine max_accel
1508 c Find the maximum difference in the accelerations of the the sites
1509 c at the beginning and the end of the time step.
1511 implicit real*8 (a-h,o-z)
1512 include 'DIMENSIONS'
1513 include 'COMMON.CONTROL'
1514 include 'COMMON.VAR'
1516 include 'COMMON.CHAIN'
1517 include 'COMMON.DERIV'
1518 include 'COMMON.GEO'
1519 include 'COMMON.LOCAL'
1520 include 'COMMON.INTERACT'
1521 include 'COMMON.IOUNITS'
1522 double precision aux(3),accel(3),accel_old(3),dacc
1524 c aux(j)=d_a(j,0)-d_a_old(j,0)
1525 accel_old(j)=d_a_old(j,0)
1532 c 7/3/08 changed to asymmetric difference
1534 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1535 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1536 accel(j)=accel(j)+0.5d0*d_a(j,i)
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 c write (iout,*) i,dacc
1541 if (dacc.gt.amax) amax=dacc
1549 accel_old(j)=d_a_old(j,0)
1554 accel_old(j)=accel_old(j)+d_a_old(j,1)
1555 accel(j)=accel(j)+d_a(j,1)
1559 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1561 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1562 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1563 accel(j)=accel(j)+d_a(j,i+nres)
1567 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1568 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1569 dacc=dabs(accel(j)-accel_old(j))
1570 c write (iout,*) "side-chain",i,dacc
1571 if (dacc.gt.amax) amax=dacc
1575 accel_old(j)=accel_old(j)+d_a_old(j,i)
1576 accel(j)=accel(j)+d_a(j,i)
1577 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1582 c---------------------------------------------------------------------
1583 subroutine predict_edrift(epdrift)
1585 c Predict the drift of the potential energy
1587 implicit real*8 (a-h,o-z)
1588 include 'DIMENSIONS'
1589 include 'COMMON.CONTROL'
1590 include 'COMMON.VAR'
1592 include 'COMMON.CHAIN'
1593 include 'COMMON.DERIV'
1594 include 'COMMON.GEO'
1595 include 'COMMON.LOCAL'
1596 include 'COMMON.INTERACT'
1597 include 'COMMON.IOUNITS'
1598 include 'COMMON.MUCA'
1599 double precision epdrift,epdriftij
1600 c Drift of the potential energy
1606 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1607 if (lmuca) epdriftij=epdriftij*factor
1608 c write (iout,*) "back",i,j,epdriftij
1609 if (epdriftij.gt.epdrift) epdrift=epdriftij
1613 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1616 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1617 if (lmuca) epdriftij=epdriftij*factor
1618 c write (iout,*) "side",i,j,epdriftij
1619 if (epdriftij.gt.epdrift) epdrift=epdriftij
1623 epdrift=0.5d0*epdrift*d_time*d_time
1624 c write (iout,*) "epdrift",epdrift
1627 c-----------------------------------------------------------------------
1628 subroutine verlet_bath
1630 c Coupling to the thermostat by using the Berendsen algorithm
1632 implicit real*8 (a-h,o-z)
1633 include 'DIMENSIONS'
1634 include 'COMMON.CONTROL'
1635 include 'COMMON.VAR'
1637 include 'COMMON.CHAIN'
1638 include 'COMMON.DERIV'
1639 include 'COMMON.GEO'
1640 include 'COMMON.LOCAL'
1641 include 'COMMON.INTERACT'
1642 include 'COMMON.IOUNITS'
1643 include 'COMMON.NAMES'
1644 double precision T_half,fact
1646 T_half=2.0d0/(dimen3*Rb)*EK
1647 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1648 c write(iout,*) "T_half", T_half
1649 c write(iout,*) "EK", EK
1650 c write(iout,*) "fact", fact
1652 d_t(j,0)=fact*d_t(j,0)
1656 d_t(j,i)=fact*d_t(j,i)
1660 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1663 d_t(j,inres)=fact*d_t(j,inres)
1669 c---------------------------------------------------------
1671 c Set up the initial conditions of a MD simulation
1672 implicit real*8 (a-h,o-z)
1673 include 'DIMENSIONS'
1677 integer IERROR,ERRCODE
1679 include 'COMMON.SETUP'
1680 include 'COMMON.CONTROL'
1681 include 'COMMON.VAR'
1684 include 'COMMON.LANGEVIN'
1686 include 'COMMON.LANGEVIN.lang0'
1688 include 'COMMON.CHAIN'
1689 include 'COMMON.DERIV'
1690 include 'COMMON.GEO'
1691 include 'COMMON.LOCAL'
1692 include 'COMMON.INTERACT'
1693 include 'COMMON.IOUNITS'
1694 include 'COMMON.NAMES'
1695 include 'COMMON.REMD'
1696 real*8 energia_long(0:n_ene),
1697 & energia_short(0:n_ene),vcm(3),incr(3)
1698 double precision cm(3),L(3),xv,sigv,lowb,highb
1699 double precision varia(maxvar)
1707 c write(iout,*) "d_time", d_time
1708 c Compute the standard deviations of stochastic forces for Langevin dynamics
1709 c if the friction coefficients do not depend on surface area
1710 if (lang.gt.0 .and. .not.surfarea) then
1712 stdforcp(i)=stdfp*dsqrt(gamp)
1715 stdforcsc(i)=stdfsc(iabs(itype(i)))
1716 & *dsqrt(gamsc(iabs(itype(i))))
1719 c Open the pdb file for snapshotshots
1722 if (ilen(tmpdir).gt.0)
1723 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1724 & liczba(:ilen(liczba))//".pdb")
1726 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1730 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1731 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1732 & liczba(:ilen(liczba))//".x")
1733 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1736 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1737 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1738 & liczba(:ilen(liczba))//".cx")
1739 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1745 if (ilen(tmpdir).gt.0)
1746 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1747 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1749 if (ilen(tmpdir).gt.0)
1750 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1751 cartname=prefix(:ilen(prefix))//"_MD.cx"
1755 write (qstr,'(256(1h ))')
1758 iq = qinfrag(i,iset)*10
1759 iw = wfrag(i,iset)/100
1761 if(me.eq.king.or..not.out1file)
1762 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1763 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1768 iq = qinpair(i,iset)*10
1769 iw = wpair(i,iset)/100
1771 if(me.eq.king.or..not.out1file)
1772 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1773 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1777 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1779 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1781 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1783 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1787 if (restart1file) then
1789 & inquire(file=mremd_rst_name,exist=file_exist)
1791 write (*,*) me," Before broadcast: file_exist",file_exist
1792 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1794 write (*,*) me," After broadcast: file_exist",file_exist
1795 c inquire(file=mremd_rst_name,exist=file_exist)
1797 if(me.eq.king.or..not.out1file)
1798 & write(iout,*) "Initial state read by master and distributed"
1800 if (ilen(tmpdir).gt.0)
1801 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1802 & //liczba(:ilen(liczba))//'.rst')
1803 inquire(file=rest2name,exist=file_exist)
1806 if(.not.restart1file) then
1807 if(me.eq.king.or..not.out1file)
1808 & write(iout,*) "Initial state will be read from file ",
1809 & rest2name(:ilen(rest2name))
1812 call rescale_weights(t_bath)
1814 if(me.eq.king.or..not.out1file)then
1815 if (restart1file) then
1816 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1819 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1822 write(iout,*) "Initial velocities randomly generated"
1829 c Generate initial velocities
1830 if(me.eq.king.or..not.out1file)
1831 & write(iout,*) "Initial velocities randomly generated"
1834 CtotTafm is the variable for AFM time which eclipsed during
1837 c rest2name = prefix(:ilen(prefix))//'.rst'
1838 if(me.eq.king.or..not.out1file)then
1839 write (iout,*) "Initial velocities"
1841 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1842 & (d_t(j,i+nres),j=1,3)
1844 c Zeroing the total angular momentum of the system
1845 write(iout,*) "Calling the zero-angular
1846 & momentum subroutine"
1849 c Getting the potential energy and forces and velocities and accelerations
1851 c write (iout,*) "velocity of the center of the mass:"
1852 c write (iout,*) (vcm(j),j=1,3)
1854 d_t(j,0)=d_t(j,0)-vcm(j)
1856 c Removing the velocity of the center of mass
1858 if(me.eq.king.or..not.out1file)then
1859 write (iout,*) "vcm right after adjustment:"
1860 write (iout,*) (vcm(j),j=1,3)
1864 if(iranconf.ne.0) then
1866 print *, 'Calling OVERLAP_SC'
1867 call overlap_sc(fail)
1871 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1872 print *,'SC_move',nft_sc,etot
1873 if(me.eq.king.or..not.out1file)
1874 & write(iout,*) 'SC_move',nft_sc,etot
1878 print *, 'Calling MINIM_DC'
1879 call minim_dc(etot,iretcode,nfun)
1881 call geom_to_var(nvar,varia)
1882 print *,'Calling MINIMIZE.'
1883 call minimize(etot,varia,iretcode,nfun)
1884 call var_to_geom(nvar,varia)
1886 if(me.eq.king.or..not.out1file)
1887 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1890 call chainbuild_cart
1895 kinetic_T=2.0d0/(dimen3*Rb)*EK
1896 if(me.eq.king.or..not.out1file)then
1906 call etotal(potEcomp)
1907 call enerprint(potEcomp)
1908 if (large) call enerprint(potEcomp)
1911 t_etotal=t_etotal+MPI_Wtime()-tt0
1913 t_etotal=t_etotal+tcpu()-tt0
1917 if(tnp .or. tnp1) then
1920 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
1922 write(iout,*) 'H0= ',H0
1926 HNose1=Hnose_nh(EK,potE)
1928 write (iout,*) 'H0= ',H0
1934 if (amax*d_time .gt. dvmax) then
1935 d_time=d_time*dvmax/amax
1936 if(me.eq.king.or..not.out1file) write (iout,*)
1937 & "Time step reduced to",d_time,
1938 & " because of too large initial acceleration."
1940 C if(me.eq.king.or..not.out1file)then
1941 C write(iout,*) "Potential energy and its components"
1942 C call enerprint(potEcomp)
1943 c write(iout,*) (potEcomp(i),i=0,n_ene)
1945 potE=potEcomp(0)-potEcomp(20)
1948 if (ntwe.ne.0) call statout(itime)
1949 if(me.eq.king.or..not.out1file)
1950 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1951 & " Kinetic energy",EK," Potential energy",potE,
1952 & " Total energy",totE," Maximum acceleration ",
1955 write (iout,*) "Initial coordinates"
1957 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1958 & (c(j,i+nres),j=1,3)
1960 write (iout,*) "Initial dC"
1962 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1963 & (dc(j,i+nres),j=1,3)
1965 write (iout,*) "Initial velocities"
1966 write (iout,"(13x,' backbone ',23x,' side chain')")
1968 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1969 & (d_t(j,i+nres),j=1,3)
1971 write (iout,*) "Initial accelerations"
1973 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1974 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1975 & (d_a(j,i+nres),j=1,3)
1981 d_t_old(j,i)=d_t(j,i)
1982 d_a_old(j,i)=d_a(j,i)
1984 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1993 call etotal_short(energia_short)
1994 if (large) call enerprint(potEcomp)
1997 t_eshort=t_eshort+MPI_Wtime()-tt0
1999 t_eshort=t_eshort+tcpu()-tt0
2002 if(tnp .or. tnp1) then
2003 E_short=energia_short(0)
2004 HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen3)
2007 c_new_var_csplit Csplit=H0-E_long
2008 c Csplit = H0-energia_short(0)
2009 write(iout,*) 'Csplit= ',Csplit
2014 if(.not.out1file .and. large) then
2015 write (iout,*) "energia_long",energia_long(0),
2016 & " energia_short",energia_short(0),
2017 & " total",energia_long(0)+energia_short(0)
2018 write (iout,*) "Initial fast-force accelerations"
2020 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2021 & (d_a(j,i+nres),j=1,3)
2024 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2027 d_a_short(j,i)=d_a(j,i)
2036 call etotal_long(energia_long)
2037 if (large) call enerprint(potEcomp)
2040 t_elong=t_elong+MPI_Wtime()-tt0
2042 t_elong=t_elong+tcpu()-tt0
2047 if(.not.out1file .and. large) then
2048 write (iout,*) "energia_long",energia_long(0)
2049 write (iout,*) "Initial slow-force accelerations"
2051 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2052 & (d_a(j,i+nres),j=1,3)
2056 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2058 t_enegrad=t_enegrad+tcpu()-tt0
2063 c-----------------------------------------------------------
2064 subroutine random_vel
2065 implicit real*8 (a-h,o-z)
2066 include 'DIMENSIONS'
2067 include 'COMMON.CONTROL'
2068 include 'COMMON.VAR'
2071 include 'COMMON.LANGEVIN'
2073 include 'COMMON.LANGEVIN.lang0'
2075 include 'COMMON.CHAIN'
2076 include 'COMMON.DERIV'
2077 include 'COMMON.GEO'
2078 include 'COMMON.LOCAL'
2079 include 'COMMON.INTERACT'
2080 include 'COMMON.IOUNITS'
2081 include 'COMMON.NAMES'
2082 include 'COMMON.TIME1'
2083 double precision xv,sigv,lowb,highb,vec_afm(3)
2084 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2085 c First generate velocities in the eigenspace of the G matrix
2086 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2093 sigv=dsqrt((Rb*t_bath)/geigen(i))
2096 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2098 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2099 c & " d_t_work_new",d_t_work_new(ii)
2102 C if (SELFGUIDE.gt.0) then
2105 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
2106 C distance=distance+vec_afm(j)**2
2108 C distance=dsqrt(distance)
2110 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
2111 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
2112 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
2113 C & d_t_work_new(j+(afmend-1)*3)
2124 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2127 c write (iout,*) "Ek from eigenvectors",Ek1
2129 c Transform velocities to UNRES coordinate space
2135 d_t_work(ind)=d_t_work(ind)
2136 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2138 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2142 c Transfer to the d_t vector
2144 d_t(j,0)=d_t_work(j)
2150 d_t(j,i)=d_t_work(ind)
2154 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2157 d_t(j,i+nres)=d_t_work(ind)
2162 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2163 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2168 c-----------------------------------------------------------
2169 subroutine sd_verlet_p_setup
2170 c Sets up the parameters of stochastic Verlet algorithm
2171 implicit real*8 (a-h,o-z)
2172 include 'DIMENSIONS'
2176 include 'COMMON.CONTROL'
2177 include 'COMMON.VAR'
2180 include 'COMMON.LANGEVIN'
2182 include 'COMMON.LANGEVIN.lang0'
2184 include 'COMMON.CHAIN'
2185 include 'COMMON.DERIV'
2186 include 'COMMON.GEO'
2187 include 'COMMON.LOCAL'
2188 include 'COMMON.INTERACT'
2189 include 'COMMON.IOUNITS'
2190 include 'COMMON.NAMES'
2191 include 'COMMON.TIME1'
2192 double precision emgdt(MAXRES6),
2193 & pterm,vterm,rho,rhoc,vsig,
2194 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2195 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2196 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2197 logical lprn /.false./
2198 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2199 double precision ktm
2206 c AL 8/17/04 Code adapted from tinker
2208 c Get the frictional and random terms for stochastic dynamics in the
2209 c eigenspace of mass-scaled UNRES friction matrix
2212 gdt = fricgam(i) * d_time
2214 c Stochastic dynamics reduces to simple MD for zero friction
2216 if (gdt .le. zero) then
2217 pfric_vec(i) = 1.0d0
2218 vfric_vec(i) = d_time
2219 afric_vec(i) = 0.5d0 * d_time * d_time
2220 prand_vec(i) = 0.0d0
2221 vrand_vec1(i) = 0.0d0
2222 vrand_vec2(i) = 0.0d0
2224 c Analytical expressions when friction coefficient is large
2227 if (gdt .ge. gdt_radius) then
2230 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2231 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2232 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2233 vterm = 1.0d0 - egdt**2
2234 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2236 c Use series expansions when friction coefficient is small
2247 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2248 & - gdt5/120.0d0 + gdt6/720.0d0
2249 & - gdt7/5040.0d0 + gdt8/40320.0d0
2250 & - gdt9/362880.0d0) / fricgam(i)**2
2251 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2252 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2253 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2254 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2255 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2256 & + 127.0d0*gdt9/90720.0d0
2257 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2258 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2259 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2260 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2261 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2262 & - 17.0d0*gdt2/1280.0d0
2263 & + 17.0d0*gdt3/6144.0d0
2264 & + 40967.0d0*gdt4/34406400.0d0
2265 & - 57203.0d0*gdt5/275251200.0d0
2266 & - 1429487.0d0*gdt6/13212057600.0d0)
2269 c Compute the scaling factors of random terms for the nonzero friction case
2271 ktm = 0.5d0*d_time/fricgam(i)
2272 psig = dsqrt(ktm*pterm) / fricgam(i)
2273 vsig = dsqrt(ktm*vterm)
2274 rhoc = dsqrt(1.0d0 - rho*rho)
2276 vrand_vec1(i) = vsig * rho
2277 vrand_vec2(i) = vsig * rhoc
2282 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2285 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2286 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2290 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2293 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2294 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2295 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2296 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2297 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2298 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2301 t_sdsetup=t_sdsetup+MPI_Wtime()
2303 t_sdsetup=t_sdsetup+tcpu()-tt0
2307 c-------------------------------------------------------------
2308 subroutine eigtransf1(n,ndim,ab,d,c)
2311 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2317 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2323 c-------------------------------------------------------------
2324 subroutine eigtransf(n,ndim,a,b,d,c)
2327 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2333 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2339 c-------------------------------------------------------------
2340 subroutine sd_verlet1
2341 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2342 implicit real*8 (a-h,o-z)
2343 include 'DIMENSIONS'
2344 include 'COMMON.CONTROL'
2345 include 'COMMON.VAR'
2348 include 'COMMON.LANGEVIN'
2350 include 'COMMON.LANGEVIN.lang0'
2352 include 'COMMON.CHAIN'
2353 include 'COMMON.DERIV'
2354 include 'COMMON.GEO'
2355 include 'COMMON.LOCAL'
2356 include 'COMMON.INTERACT'
2357 include 'COMMON.IOUNITS'
2358 include 'COMMON.NAMES'
2359 double precision stochforcvec(MAXRES6)
2360 common /stochcalc/ stochforcvec
2361 logical lprn /.false./
2363 c write (iout,*) "dc_old"
2365 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2366 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2369 dc_work(j)=dc_old(j,0)
2370 d_t_work(j)=d_t_old(j,0)
2371 d_a_work(j)=d_a_old(j,0)
2376 dc_work(ind+j)=dc_old(j,i)
2377 d_t_work(ind+j)=d_t_old(j,i)
2378 d_a_work(ind+j)=d_a_old(j,i)
2383 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2385 dc_work(ind+j)=dc_old(j,i+nres)
2386 d_t_work(ind+j)=d_t_old(j,i+nres)
2387 d_a_work(ind+j)=d_a_old(j,i+nres)
2395 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2399 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2400 & vfric_mat(i,j),afric_mat(i,j),
2401 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2409 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2410 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2411 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2412 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2414 d_t_work_new(i)=ddt1+0.5d0*ddt2
2415 d_t_work(i)=ddt1+ddt2
2420 d_t(j,0)=d_t_work(j)
2425 dc(j,i)=dc_work(ind+j)
2426 d_t(j,i)=d_t_work(ind+j)
2431 if (itype(i).ne.10) then
2434 dc(j,inres)=dc_work(ind+j)
2435 d_t(j,inres)=d_t_work(ind+j)
2442 c--------------------------------------------------------------------------
2443 subroutine sd_verlet2
2444 c Calculating the adjusted velocities for accelerations
2445 implicit real*8 (a-h,o-z)
2446 include 'DIMENSIONS'
2447 include 'COMMON.CONTROL'
2448 include 'COMMON.VAR'
2451 include 'COMMON.LANGEVIN'
2453 include 'COMMON.LANGEVIN.lang0'
2455 include 'COMMON.CHAIN'
2456 include 'COMMON.DERIV'
2457 include 'COMMON.GEO'
2458 include 'COMMON.LOCAL'
2459 include 'COMMON.INTERACT'
2460 include 'COMMON.IOUNITS'
2461 include 'COMMON.NAMES'
2462 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2463 common /stochcalc/ stochforcvec
2465 c Compute the stochastic forces which contribute to velocity change
2467 call stochastic_force(stochforcvecV)
2474 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2475 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2476 & vrand_mat2(i,j)*stochforcvecV(j)
2478 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2482 d_t(j,0)=d_t_work(j)
2487 d_t(j,i)=d_t_work(ind+j)
2492 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2495 d_t(j,inres)=d_t_work(ind+j)
2502 c-----------------------------------------------------------
2503 subroutine sd_verlet_ciccotti_setup
2504 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2506 implicit real*8 (a-h,o-z)
2507 include 'DIMENSIONS'
2511 include 'COMMON.CONTROL'
2512 include 'COMMON.VAR'
2515 include 'COMMON.LANGEVIN'
2517 include 'COMMON.LANGEVIN.lang0'
2519 include 'COMMON.CHAIN'
2520 include 'COMMON.DERIV'
2521 include 'COMMON.GEO'
2522 include 'COMMON.LOCAL'
2523 include 'COMMON.INTERACT'
2524 include 'COMMON.IOUNITS'
2525 include 'COMMON.NAMES'
2526 include 'COMMON.TIME1'
2527 double precision emgdt(MAXRES6),
2528 & pterm,vterm,rho,rhoc,vsig,
2529 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2530 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2531 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2532 logical lprn /.false./
2533 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2534 double precision ktm
2541 c AL 8/17/04 Code adapted from tinker
2543 c Get the frictional and random terms for stochastic dynamics in the
2544 c eigenspace of mass-scaled UNRES friction matrix
2547 write (iout,*) "i",i," fricgam",fricgam(i)
2548 gdt = fricgam(i) * d_time
2550 c Stochastic dynamics reduces to simple MD for zero friction
2552 if (gdt .le. zero) then
2553 pfric_vec(i) = 1.0d0
2554 vfric_vec(i) = d_time
2555 afric_vec(i) = 0.5d0*d_time*d_time
2556 prand_vec(i) = afric_vec(i)
2557 vrand_vec2(i) = vfric_vec(i)
2559 c Analytical expressions when friction coefficient is large
2564 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2565 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2566 prand_vec(i) = afric_vec(i)
2567 vrand_vec2(i) = vfric_vec(i)
2569 c Compute the scaling factors of random terms for the nonzero friction case
2571 c ktm = 0.5d0*d_time/fricgam(i)
2572 c psig = dsqrt(ktm*pterm) / fricgam(i)
2573 c vsig = dsqrt(ktm*vterm)
2574 c prand_vec(i) = psig*afric_vec(i)
2575 c vrand_vec2(i) = vsig*vfric_vec(i)
2580 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2583 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2584 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2588 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2590 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2591 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2592 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2593 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2594 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2596 t_sdsetup=t_sdsetup+MPI_Wtime()
2598 t_sdsetup=t_sdsetup+tcpu()-tt0
2602 c-------------------------------------------------------------
2603 subroutine sd_verlet1_ciccotti
2604 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2605 implicit real*8 (a-h,o-z)
2606 include 'DIMENSIONS'
2610 include 'COMMON.CONTROL'
2611 include 'COMMON.VAR'
2614 include 'COMMON.LANGEVIN'
2616 include 'COMMON.LANGEVIN.lang0'
2618 include 'COMMON.CHAIN'
2619 include 'COMMON.DERIV'
2620 include 'COMMON.GEO'
2621 include 'COMMON.LOCAL'
2622 include 'COMMON.INTERACT'
2623 include 'COMMON.IOUNITS'
2624 include 'COMMON.NAMES'
2625 double precision stochforcvec(MAXRES6)
2626 common /stochcalc/ stochforcvec
2627 logical lprn /.false./
2629 c write (iout,*) "dc_old"
2631 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2632 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2635 dc_work(j)=dc_old(j,0)
2636 d_t_work(j)=d_t_old(j,0)
2637 d_a_work(j)=d_a_old(j,0)
2642 dc_work(ind+j)=dc_old(j,i)
2643 d_t_work(ind+j)=d_t_old(j,i)
2644 d_a_work(ind+j)=d_a_old(j,i)
2649 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2651 dc_work(ind+j)=dc_old(j,i+nres)
2652 d_t_work(ind+j)=d_t_old(j,i+nres)
2653 d_a_work(ind+j)=d_a_old(j,i+nres)
2662 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2666 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2667 & vfric_mat(i,j),afric_mat(i,j),
2668 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2676 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2677 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2678 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2679 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2681 d_t_work_new(i)=ddt1+0.5d0*ddt2
2682 d_t_work(i)=ddt1+ddt2
2687 d_t(j,0)=d_t_work(j)
2692 dc(j,i)=dc_work(ind+j)
2693 d_t(j,i)=d_t_work(ind+j)
2698 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2701 dc(j,inres)=dc_work(ind+j)
2702 d_t(j,inres)=d_t_work(ind+j)
2709 c--------------------------------------------------------------------------
2710 subroutine sd_verlet2_ciccotti
2711 c Calculating the adjusted velocities for accelerations
2712 implicit real*8 (a-h,o-z)
2713 include 'DIMENSIONS'
2714 include 'COMMON.CONTROL'
2715 include 'COMMON.VAR'
2718 include 'COMMON.LANGEVIN'
2720 include 'COMMON.LANGEVIN.lang0'
2722 include 'COMMON.CHAIN'
2723 include 'COMMON.DERIV'
2724 include 'COMMON.GEO'
2725 include 'COMMON.LOCAL'
2726 include 'COMMON.INTERACT'
2727 include 'COMMON.IOUNITS'
2728 include 'COMMON.NAMES'
2729 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2730 common /stochcalc/ stochforcvec
2732 c Compute the stochastic forces which contribute to velocity change
2734 call stochastic_force(stochforcvecV)
2741 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2742 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2743 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2745 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2749 d_t(j,0)=d_t_work(j)
2754 d_t(j,i)=d_t_work(ind+j)
2759 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2762 d_t(j,inres)=d_t_work(ind+j)
2770 double precision function HNose(ek,s,e,pi,Q,t_bath,dimenl)
2772 double precision ek,s,e,pi,Q,t_bath,Rb
2775 HNose=ek+e+pi**2/(2*Q)+dimenl*Rb*t_bath*log(s)
2776 c print '(6f15.5,i5,a2,2f15.5)',ek,s,e,pi,Q,t_bath,dimenl,"--",
2777 c & pi**2/(2*Q),dimenl*Rb*t_bath*log(s)
2780 c-----------------------------------------------------------------
2781 double precision function HNose_nh(eki,e)
2782 implicit real*8 (a-h,o-z)
2783 include 'DIMENSIONS'
2785 HNose_nh=eki+e+dimen3*Rb*t_bath*xlogs(1)+qmass(1)*vlogs(1)**2/2
2787 HNose_nh=HNose_nh+qmass(i)*vlogs(i)**2/2+Rb*t_bath*xlogs(i)
2789 c write(4,'(5e15.5)')
2790 c & vlogs(1),xlogs(1),HNose,eki,e
2793 c-----------------------------------------------------------------
2794 SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8)
2795 implicit real*8 (a-h,o-z)
2796 include 'DIMENSIONS'
2798 double precision akin,gnkt,dt,aa,gkt,scale
2799 double precision wdti(maxyosh),wdti2(maxyosh),
2800 & wdti4(maxyosh),wdti8(maxyosh)
2801 integer i,iresn,iyosh,inos,nnos1
2810 C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE
2811 C INTEGRATION FROM t=0 TO t=DT/2
2812 C GET THE TOTAL KINETIC ENERGY
2814 c CALL GETKINP(MASS,VX,VY,VZ,AKIN)
2816 GLOGS(1) = (AKIN - GNKT)/QMASS(1)
2817 C START THE MULTIPLE TIME STEP PROCEDURE
2820 C UPDATE THE THERMOSTAT VELOCITIES
2821 VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2823 AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) )
2824 VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA
2825 & + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA
2827 C UPDATE THE PARTICLE VELOCITIES
2828 AA = EXP(-WDTI2(IYOSH)*VLOGS(1) )
2831 GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1)
2832 C UPDATE THE THERMOSTAT POSITIONS
2834 XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH)
2836 C UPDATE THE THERMOSTAT VELOCITIES
2838 AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) )
2839 VLOGS(INOS) = VLOGS(INOS)*AA*AA
2840 & + WDTI4(IYOSH)*GLOGS(INOS)*AA
2841 GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS)
2842 & -GKT)/QMASS(INOS+1)
2844 VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2847 C UPDATE THE PARTICLE VELOCITIES
2848 c outside of this subroutine
2850 c VX(I) = VX(I)*SCALE
2851 c VY(I) = VY(I)*SCALE
2852 c VZ(I) = VZ(I)*SCALE
2856 c-----------------------------------------------------------------
2857 subroutine tnp1_respa_i_step1
2858 c Applying Nose-Poincare algorithm - step 1 to coordinates
2859 c JPSJ 70 75 (2001) S. Nose
2861 c d_t is not updated here
2863 implicit real*8 (a-h,o-z)
2864 include 'DIMENSIONS'
2865 include 'COMMON.CONTROL'
2866 include 'COMMON.VAR'
2868 include 'COMMON.CHAIN'
2869 include 'COMMON.DERIV'
2870 include 'COMMON.GEO'
2871 include 'COMMON.LOCAL'
2872 include 'COMMON.INTERACT'
2873 include 'COMMON.IOUNITS'
2874 include 'COMMON.NAMES'
2875 double precision adt,adt2,tmp
2877 tmp=1+pi_np/(2*Q_np)*0.5*d_time
2880 s12_dt=d_time/s12_np
2881 d_time_s12=d_time*0.5*s12_np
2884 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2885 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2889 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2890 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2894 if (itype(i).ne.10) then
2897 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2898 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2904 c---------------------------------------------------------------------
2905 subroutine tnp1_respa_i_step2
2906 c Step 2 of the velocity Verlet algorithm: update velocities
2907 implicit real*8 (a-h,o-z)
2908 include 'DIMENSIONS'
2909 include 'COMMON.CONTROL'
2910 include 'COMMON.VAR'
2912 include 'COMMON.CHAIN'
2913 include 'COMMON.DERIV'
2914 include 'COMMON.GEO'
2915 include 'COMMON.LOCAL'
2916 include 'COMMON.INTERACT'
2917 include 'COMMON.IOUNITS'
2918 include 'COMMON.NAMES'
2920 double precision d_time_s12
2924 d_t(j,i)=d_t_new(j,i)
2931 d_time_s12=0.5d0*s12_np*d_time
2934 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
2938 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
2942 if (itype(i).ne.10) then
2945 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
2950 pistar=pistar+(EK-0.5*(E_old+potE)
2951 & -dimen3*Rb*t_bath*log(s12_np)+Csplit-dimen3*Rb*t_bath)*d_time
2952 tmp=1+pistar/(2*Q_np)*0.5*d_time
2958 c-------------------------------------------------------
2960 subroutine tnp1_step1
2961 c Applying Nose-Poincare algorithm - step 1 to coordinates
2962 c JPSJ 70 75 (2001) S. Nose
2964 c d_t is not updated here
2966 implicit real*8 (a-h,o-z)
2967 include 'DIMENSIONS'
2968 include 'COMMON.CONTROL'
2969 include 'COMMON.VAR'
2971 include 'COMMON.CHAIN'
2972 include 'COMMON.DERIV'
2973 include 'COMMON.GEO'
2974 include 'COMMON.LOCAL'
2975 include 'COMMON.INTERACT'
2976 include 'COMMON.IOUNITS'
2977 include 'COMMON.NAMES'
2978 double precision adt,adt2,tmp
2980 tmp=1+pi_np/(2*Q_np)*0.5*d_time
2983 s12_dt=d_time/s12_np
2984 d_time_s12=d_time*0.5*s12_np
2987 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2988 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2992 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2993 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2997 if (itype(i).ne.10) then
3000 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
3001 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
3007 c---------------------------------------------------------------------
3008 subroutine tnp1_step2
3009 c Step 2 of the velocity Verlet algorithm: update velocities
3010 implicit real*8 (a-h,o-z)
3011 include 'DIMENSIONS'
3012 include 'COMMON.CONTROL'
3013 include 'COMMON.VAR'
3015 include 'COMMON.CHAIN'
3016 include 'COMMON.DERIV'
3017 include 'COMMON.GEO'
3018 include 'COMMON.LOCAL'
3019 include 'COMMON.INTERACT'
3020 include 'COMMON.IOUNITS'
3021 include 'COMMON.NAMES'
3023 double precision d_time_s12
3027 d_t(j,i)=d_t_new(j,i)
3034 d_time_s12=0.5d0*s12_np*d_time
3037 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
3041 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
3045 if (itype(i).ne.10) then
3048 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
3053 cd write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np
3054 pistar=pistar+(EK-0.5*(E_old+potE)
3055 & -dimen3*Rb*t_bath*log(s12_np)+H0-dimen3*Rb*t_bath)*d_time
3056 tmp=1+pistar/(2*Q_np)*0.5*d_time
3063 c-----------------------------------------------------------------
3064 subroutine tnp_respa_i_step1
3065 c Applying Nose-Poincare algorithm - step 1 to coordinates
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
3084 d_time_s=d_time*0.5*s_np
3085 ct2 d_time_s=d_time*0.5*s12_np
3088 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3092 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3096 if (itype(i).ne.10) then
3099 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3106 d_t(j,i)=d_t_new(j,i)
3113 C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-Csplit)
3116 pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3117 tmp=0.5*d_time*pistar/Q_np
3118 s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3120 d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3121 ct2 d_time_ss=d_time/s12_np
3122 c d_time_ss=0.5*d_time*(1.0/sold_np+1.0/s_np)
3125 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3129 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3133 if (itype(i).ne.10) then
3136 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3143 c---------------------------------------------------------------------
3145 subroutine tnp_respa_i_step2
3146 c Step 2 of the velocity Verlet algorithm: update velocities
3147 implicit real*8 (a-h,o-z)
3148 include 'DIMENSIONS'
3149 include 'COMMON.CONTROL'
3150 include 'COMMON.VAR'
3152 include 'COMMON.CHAIN'
3153 include 'COMMON.DERIV'
3154 include 'COMMON.GEO'
3155 include 'COMMON.LOCAL'
3156 include 'COMMON.INTERACT'
3157 include 'COMMON.IOUNITS'
3158 include 'COMMON.NAMES'
3160 double precision d_time_s
3162 EK=EK*(s_np/s12_np)**2
3163 HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3164 pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath
3167 cr print '(a,5f)','i_step2',EK,potE,HNose1,pi_np,E_long
3168 d_time_s=d_time*0.5*s12_np
3169 c d_time_s=d_time*0.5*s_np
3172 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3176 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3180 if (itype(i).ne.10) then
3183 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
3192 c-----------------------------------------------------------------
3193 subroutine tnp_respa_step1
3194 c Applying Nose-Poincare algorithm - step 1 to vel for RESPA
3195 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3197 c d_t is not updated here, it is destroyed
3199 implicit real*8 (a-h,o-z)
3200 include 'DIMENSIONS'
3201 include 'COMMON.CONTROL'
3202 include 'COMMON.VAR'
3204 include 'COMMON.CHAIN'
3205 include 'COMMON.DERIV'
3206 include 'COMMON.GEO'
3207 include 'COMMON.LOCAL'
3208 include 'COMMON.INTERACT'
3209 include 'COMMON.IOUNITS'
3210 include 'COMMON.NAMES'
3211 double precision C_np,d_time_s,tmp,d_time_ss
3212 double precision energia(0:n_ene)
3214 d_time_s=d_time*0.5*s_np
3217 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3221 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3225 if (itype(i).ne.10) then
3228 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3234 c C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3237 c pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3238 c tmp=0.5*d_time*pistar/Q_np
3239 c s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3240 c write(iout,*) 'tnp_respa_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3246 c-------------------------------------
3247 c test of reviewer's comment
3248 pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3249 cr print '(a,3f)','1 pi_np,s_np',pi_np,s_np,E_long
3250 c-------------------------------------
3254 c---------------------------------------------------------------------
3255 subroutine tnp_respa_step2
3256 c Step 2 of the velocity Verlet algorithm: update velocities for RESPA
3257 implicit real*8 (a-h,o-z)
3258 include 'DIMENSIONS'
3259 include 'COMMON.CONTROL'
3260 include 'COMMON.VAR'
3262 include 'COMMON.CHAIN'
3263 include 'COMMON.DERIV'
3264 include 'COMMON.GEO'
3265 include 'COMMON.LOCAL'
3266 include 'COMMON.INTERACT'
3267 include 'COMMON.IOUNITS'
3268 include 'COMMON.NAMES'
3270 double precision d_time_s
3276 ct HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3277 ct pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3278 ct & -0.5*d_time*(HNose1-H0)
3280 c-------------------------------------
3281 c test of reviewer's comment
3282 pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3283 cr print '(a,3f)','2 pi_np,s_np',pi_np,s_np,E_long
3284 c-------------------------------------
3285 d_time_s=d_time*0.5*s_np
3288 d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3292 d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3296 if (itype(i).ne.10) then
3299 d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3308 c---------------------------------------------------------------------
3309 subroutine tnp_step1
3310 c Applying Nose-Poincare algorithm - step 1 to coordinates
3311 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3313 c d_t is not updated here, it is destroyed
3315 implicit real*8 (a-h,o-z)
3316 include 'DIMENSIONS'
3317 include 'COMMON.CONTROL'
3318 include 'COMMON.VAR'
3320 include 'COMMON.CHAIN'
3321 include 'COMMON.DERIV'
3322 include 'COMMON.GEO'
3323 include 'COMMON.LOCAL'
3324 include 'COMMON.INTERACT'
3325 include 'COMMON.IOUNITS'
3326 include 'COMMON.NAMES'
3327 double precision C_np,d_time_s,tmp,d_time_ss
3329 d_time_s=d_time*0.5*s_np
3332 d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3336 d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3340 if (itype(i).ne.10) then
3343 d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3350 d_t(j,i)=d_t_new(j,i)
3357 C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3360 pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3361 tmp=0.5*d_time*pistar/Q_np
3362 s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3363 c write(iout,*) 'tnp_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3365 d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3368 dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3372 dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3376 if (itype(i).ne.10) then
3379 dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3386 c-----------------------------------------------------------------
3387 subroutine tnp_step2
3388 c Step 2 of the velocity Verlet algorithm: update velocities
3389 implicit real*8 (a-h,o-z)
3390 include 'DIMENSIONS'
3391 include 'COMMON.CONTROL'
3392 include 'COMMON.VAR'
3394 include 'COMMON.CHAIN'
3395 include 'COMMON.DERIV'
3396 include 'COMMON.GEO'
3397 include 'COMMON.LOCAL'
3398 include 'COMMON.INTERACT'
3399 include 'COMMON.IOUNITS'
3400 include 'COMMON.NAMES'
3402 double precision d_time_s
3404 EK=EK*(s_np/s12_np)**2
3405 HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3406 pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3407 & -0.5*d_time*(HNose1-H0)
3409 cd write(iout,'(a,4f)') 'mmm',EK,potE,HNose1,pi_np
3410 d_time_s=d_time*0.5*s12_np
3413 d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3417 d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3421 if (itype(i).ne.10) then
3424 d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s