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
331 else if (lang.eq.2 .or. lang.eq.3) then
333 call stochastic_force(stochforcvec)
336 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
338 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
345 icount_scale=icount_scale+1
346 if (icount_scale.gt.maxcount_scale) then
348 & "ERROR: too many attempts at scaling down the time step. ",
349 & "amax=",amax,"epdrift=",epdrift,
350 & "damax=",damax,"edriftmax=",edriftmax,
354 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
358 c First step of the velocity Verlet algorithm
363 else if (lang.eq.3) then
365 call sd_verlet1_ciccotti
367 else if (lang.eq.1) then
372 c Build the chain from the newly calculated coordinates
374 if (rattle) call rattle1
376 if (large.and. mod(itime,ntwe).eq.0) then
377 write (iout,*) "Cartesian and internal coordinates: step 1"
382 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
383 & (dc(j,i+nres),j=1,3)
385 write (iout,*) "Accelerations"
387 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
388 & (d_a(j,i+nres),j=1,3)
390 write (iout,*) "Velocities, step 1"
392 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
393 & (d_t(j,i+nres),j=1,3)
402 c Calculate energy and forces
404 call etotal(potEcomp)
405 if (large.and. mod(itime,ntwe).eq.0)
406 & call enerprint(potEcomp)
409 t_etotal=t_etotal+MPI_Wtime()-tt0
411 t_etotal=t_etotal+tcpu()-tt0
414 potE=potEcomp(0)-potEcomp(20)
416 c Get the new accelerations
419 t_enegrad=t_enegrad+MPI_Wtime()-tt0
421 t_enegrad=t_enegrad+tcpu()-tt0
423 c Determine maximum acceleration and scale down the timestep if needed
425 amax=amax/(itime_scal+1)**2
426 call predict_edrift(epdrift)
427 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
428 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
430 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
432 itime_scal=itime_scal+ifac_time
433 c fac_time=dmin1(damax/amax,0.5d0)
434 fac_time=0.5d0**ifac_time
435 d_time=d_time*fac_time
436 if (lang.eq.2 .or. lang.eq.3) then
438 c write (iout,*) "Calling sd_verlet_setup: 1"
439 c Rescale the stochastic forces and recalculate or restore
440 c the matrices of tinker integrator
441 if (itime_scal.gt.maxflag_stoch) then
442 if (large) write (iout,'(a,i5,a)')
443 & "Calculate matrices for stochastic step;",
444 & " itime_scal ",itime_scal
446 call sd_verlet_p_setup
448 call sd_verlet_ciccotti_setup
450 write (iout,'(2a,i3,a,i3,1h.)')
451 & "Warning: cannot store matrices for stochastic",
452 & " integration because the index",itime_scal,
453 & " is greater than",maxflag_stoch
454 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
455 & " integration Langevin algorithm for better efficiency."
456 else if (flag_stoch(itime_scal)) then
457 if (large) write (iout,'(a,i5,a,l1)')
458 & "Restore matrices for stochastic step; itime_scal ",
459 & itime_scal," flag ",flag_stoch(itime_scal)
462 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
463 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
464 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
465 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
466 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
467 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
471 if (large) write (iout,'(2a,i5,a,l1)')
472 & "Calculate & store matrices for stochastic step;",
473 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
475 call sd_verlet_p_setup
477 call sd_verlet_ciccotti_setup
479 flag_stoch(ifac_time)=.true.
482 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
483 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
484 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
485 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
486 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
487 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
491 fac_time=1.0d0/dsqrt(fac_time)
493 stochforcvec(i)=fac_time*stochforcvec(i)
496 else if (lang.eq.1) then
497 c Rescale the accelerations due to stochastic forces
498 fac_time=1.0d0/dsqrt(fac_time)
500 d_as_work(i)=d_as_work(i)*fac_time
503 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
504 & "itime",itime," Timestep scaled down to ",
505 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
507 c Second step of the velocity Verlet algorithm
512 else if (lang.eq.3) then
514 call sd_verlet2_ciccotti
516 else if (lang.eq.1) then
521 if (rattle) call rattle2
524 C print *,totTafm,"TU?"
525 if (d_time.ne.d_time0) then
528 if (lang.eq.2 .or. lang.eq.3) then
529 if (large) write (iout,'(a)')
530 & "Restore original matrices for stochastic step"
531 c write (iout,*) "Calling sd_verlet_setup: 2"
532 c Restore the matrices of tinker integrator if the time step has been restored
535 pfric_mat(i,j)=pfric0_mat(i,j,0)
536 afric_mat(i,j)=afric0_mat(i,j,0)
537 vfric_mat(i,j)=vfric0_mat(i,j,0)
538 prand_mat(i,j)=prand0_mat(i,j,0)
539 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
540 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
549 c Calculate the kinetic and the total energy and the kinetic temperature
554 c write (iout,*) "step",itime," EK",EK," EK1",EK1
556 c Couple the system to Berendsen bath if needed
557 if (tbf .and. lang.eq.0) then
560 kinetic_T=2.0d0/(dimen3*Rb)*EK
561 c Backup the coordinates, velocities, and accelerations
565 d_t_old(j,i)=d_t(j,i)
566 d_a_old(j,i)=d_a(j,i)
570 if (mod(itime,ntwe).eq.0 .and. large) then
571 write (iout,*) "Velocities, step 2"
573 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
574 & (d_t(j,i+nres),j=1,3)
580 c-------------------------------------------------------------------------------
581 subroutine RESPA_step(itime)
582 c-------------------------------------------------------------------------------
583 c Perform a single RESPA step.
584 c-------------------------------------------------------------------------------
585 implicit real*8 (a-h,o-z)
589 integer IERROR,ERRCODE
591 include 'COMMON.SETUP'
592 include 'COMMON.CONTROL'
596 include 'COMMON.LANGEVIN'
598 include 'COMMON.LANGEVIN.lang0'
600 include 'COMMON.CHAIN'
601 include 'COMMON.DERIV'
603 include 'COMMON.LOCAL'
604 include 'COMMON.INTERACT'
605 include 'COMMON.IOUNITS'
606 include 'COMMON.NAMES'
607 include 'COMMON.TIME1'
608 double precision energia_short(0:n_ene),
609 & energia_long(0:n_ene)
610 double precision cm(3),L(3),vcm(3),incr(3)
611 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
612 & d_a_old0(3,0:maxres2)
613 logical PRINT_AMTS_MSG /.false./
614 integer ilen,count,rstcount
617 integer maxcount_scale /10/
619 double precision stochforcvec(MAXRES6)
620 common /stochcalc/ stochforcvec
623 common /cipiszcze/ itt
626 if (large.and. mod(itime,ntwe).eq.0) then
627 write (iout,*) "***************** RESPA itime",itime
628 write (iout,*) "Cartesian and internal coordinates: step 0"
630 call pdbout(0.0d0,"cipiszcze",iout)
632 write (iout,*) "Accelerations from long-range forces"
634 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
635 & (d_a(j,i+nres),j=1,3)
637 write (iout,*) "Velocities, step 0"
639 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
640 & (d_t(j,i+nres),j=1,3)
645 c Perform the initial RESPA step (increment velocities)
646 c write (iout,*) "*********************** RESPA ini"
649 if (mod(itime,ntwe).eq.0 .and. large) then
650 write (iout,*) "Velocities, end"
652 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
653 & (d_t(j,i+nres),j=1,3)
657 c Compute the short-range forces
663 C 7/2/2009 commented out
665 c call etotal_short(energia_short)
668 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
671 d_a(j,i)=d_a_short(j,i)
675 if (large.and. mod(itime,ntwe).eq.0) then
676 write (iout,*) "energia_short",energia_short(0)
677 write (iout,*) "Accelerations from short-range forces"
679 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
680 & (d_a(j,i+nres),j=1,3)
685 t_enegrad=t_enegrad+MPI_Wtime()-tt0
687 t_enegrad=t_enegrad+tcpu()-tt0
692 d_t_old(j,i)=d_t(j,i)
693 d_a_old(j,i)=d_a(j,i)
696 c 6/30/08 A-MTS: attempt at increasing the split number
699 dc_old0(j,i)=dc_old(j,i)
700 d_t_old0(j,i)=d_t_old(j,i)
701 d_a_old0(j,i)=d_a_old(j,i)
704 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
705 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
712 c write (iout,*) "itime",itime," ntime_split",ntime_split
713 c Split the time step
714 d_time=d_time0/ntime_split
715 c Perform the short-range RESPA steps (velocity Verlet increments of
716 c positions and velocities using short-range forces)
717 c write (iout,*) "*********************** RESPA split"
718 do itsplit=1,ntime_split
721 else if (lang.eq.2 .or. lang.eq.3) then
723 call stochastic_force(stochforcvec)
726 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
728 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
733 c First step of the velocity Verlet algorithm
738 else if (lang.eq.3) then
740 call sd_verlet1_ciccotti
742 else if (lang.eq.1) then
747 c Build the chain from the newly calculated coordinates
749 if (rattle) call rattle1
751 if (large.and. mod(itime,ntwe).eq.0) then
752 write (iout,*) "***** ITSPLIT",itsplit
753 write (iout,*) "Cartesian and internal coordinates: step 1"
754 call pdbout(0.0d0,"cipiszcze",iout)
757 write (iout,*) "Velocities, step 1"
759 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
760 & (d_t(j,i+nres),j=1,3)
769 c Calculate energy and forces
771 call etotal_short(energia_short)
772 if (large.and. mod(itime,ntwe).eq.0)
773 & call enerprint(energia_short)
776 t_eshort=t_eshort+MPI_Wtime()-tt0
778 t_eshort=t_eshort+tcpu()-tt0
782 c Get the new accelerations
784 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
787 d_a_short(j,i)=d_a(j,i)
791 if (large.and. mod(itime,ntwe).eq.0) then
792 write (iout,*)"energia_short",energia_short(0)
793 write (iout,*) "Accelerations from short-range forces"
795 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
796 & (d_a(j,i+nres),j=1,3)
801 c Determine maximum acceleration and scale down the timestep if needed
803 amax=amax/ntime_split**2
804 call predict_edrift(epdrift)
805 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
806 & write (iout,*) "amax",amax," damax",damax,
807 & " epdrift",epdrift," epdriftmax",epdriftmax
808 c Exit loop and try with increased split number if the change of
809 c acceleration is too big
810 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
811 if (ntime_split.lt.maxtime_split) then
813 ntime_split=ntime_split*2
816 dc_old(j,i)=dc_old0(j,i)
817 d_t_old(j,i)=d_t_old0(j,i)
818 d_a_old(j,i)=d_a_old0(j,i)
821 if (PRINT_AMTS_MSG) then
822 write (iout,*) "acceleration/energy drift too large",amax,
823 & epdrift," split increased to ",ntime_split," itime",itime,
829 & "Uh-hu. Bumpy landscape. Maximum splitting number",
831 & " already reached!!! Trying to carry on!"
835 t_enegrad=t_enegrad+MPI_Wtime()-tt0
837 t_enegrad=t_enegrad+tcpu()-tt0
839 c Second step of the velocity Verlet algorithm
844 else if (lang.eq.3) then
846 call sd_verlet2_ciccotti
848 else if (lang.eq.1) then
853 if (rattle) call rattle2
854 c Backup the coordinates, velocities, and accelerations
858 d_t_old(j,i)=d_t(j,i)
859 d_a_old(j,i)=d_a(j,i)
866 c Restore the time step
868 c Compute long-range forces
875 call etotal_long(energia_long)
876 if (large.and. mod(itime,ntwe).eq.0)
877 & call enerprint(energia_long)
880 t_elong=t_elong+MPI_Wtime()-tt0
882 t_elong=t_elong+tcpu()-tt0
888 t_enegrad=t_enegrad+MPI_Wtime()-tt0
890 t_enegrad=t_enegrad+tcpu()-tt0
892 c Compute accelerations from long-range forces
894 if (large.and. mod(itime,ntwe).eq.0) then
895 write (iout,*) "energia_long",energia_long(0)
896 write (iout,*) "Cartesian and internal coordinates: step 2"
898 call pdbout(0.0d0,"cipiszcze",iout)
900 write (iout,*) "Accelerations from long-range forces"
902 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
903 & (d_a(j,i+nres),j=1,3)
905 write (iout,*) "Velocities, step 2"
907 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
908 & (d_t(j,i+nres),j=1,3)
912 c Compute the final RESPA step (increment velocities)
913 c write (iout,*) "*********************** RESPA fin"
915 c Compute the complete potential energy
917 potEcomp(i)=energia_short(i)+energia_long(i)
919 potE=potEcomp(0)-potEcomp(20)
920 c potE=energia_short(0)+energia_long(0)
923 c Calculate the kinetic and the total energy and the kinetic temperature
926 c Couple the system to Berendsen bath if needed
927 if (tbf .and. lang.eq.0) then
930 kinetic_T=2.0d0/(dimen3*Rb)*EK
931 c Backup the coordinates, velocities, and accelerations
933 if (mod(itime,ntwe).eq.0 .and. large) then
934 write (iout,*) "Velocities, end"
936 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
937 & (d_t(j,i+nres),j=1,3)
943 c---------------------------------------------------------------------
945 c First and last RESPA step (incrementing velocities using long-range
947 implicit real*8 (a-h,o-z)
949 include 'COMMON.CONTROL'
952 include 'COMMON.CHAIN'
953 include 'COMMON.DERIV'
955 include 'COMMON.LOCAL'
956 include 'COMMON.INTERACT'
957 include 'COMMON.IOUNITS'
958 include 'COMMON.NAMES'
960 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
964 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
968 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
971 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
977 c-----------------------------------------------------------------
979 c Applying velocity Verlet algorithm - step 1 to coordinates
980 implicit real*8 (a-h,o-z)
982 include 'COMMON.CONTROL'
985 include 'COMMON.CHAIN'
986 include 'COMMON.DERIV'
988 include 'COMMON.LOCAL'
989 include 'COMMON.INTERACT'
990 include 'COMMON.IOUNITS'
991 include 'COMMON.NAMES'
992 double precision adt,adt2
995 write (iout,*) "VELVERLET1 START: DC"
997 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
998 & (dc(j,i+nres),j=1,3)
1002 adt=d_a_old(j,0)*d_time
1004 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1005 d_t_new(j,0)=d_t_old(j,0)+adt2
1006 d_t(j,0)=d_t_old(j,0)+adt
1012 adt=d_a_old(j,i)*d_time
1014 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1015 d_t_new(j,i)=d_t_old(j,i)+adt2
1016 d_t(j,i)=d_t_old(j,i)+adt
1021 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1024 adt=d_a_old(j,inres)*d_time
1026 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1027 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1028 d_t(j,inres)=d_t_old(j,inres)+adt
1033 write (iout,*) "VELVERLET1 END: DC"
1035 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1036 & (dc(j,i+nres),j=1,3)
1041 c---------------------------------------------------------------------
1043 c Step 2 of the velocity Verlet algorithm: update velocities
1044 implicit real*8 (a-h,o-z)
1045 include 'DIMENSIONS'
1046 include 'COMMON.CONTROL'
1047 include 'COMMON.VAR'
1049 include 'COMMON.CHAIN'
1050 include 'COMMON.DERIV'
1051 include 'COMMON.GEO'
1052 include 'COMMON.LOCAL'
1053 include 'COMMON.INTERACT'
1054 include 'COMMON.IOUNITS'
1055 include 'COMMON.NAMES'
1057 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1061 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1065 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1068 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1074 c-----------------------------------------------------------------
1075 subroutine sddir_precalc
1076 c Applying velocity Verlet algorithm - step 1 to coordinates
1077 implicit real*8 (a-h,o-z)
1078 include 'DIMENSIONS'
1082 include 'COMMON.CONTROL'
1083 include 'COMMON.VAR'
1086 include 'COMMON.LANGEVIN'
1088 include 'COMMON.LANGEVIN.lang0'
1090 include 'COMMON.CHAIN'
1091 include 'COMMON.DERIV'
1092 include 'COMMON.GEO'
1093 include 'COMMON.LOCAL'
1094 include 'COMMON.INTERACT'
1095 include 'COMMON.IOUNITS'
1096 include 'COMMON.NAMES'
1097 include 'COMMON.TIME1'
1098 double precision stochforcvec(MAXRES6)
1099 common /stochcalc/ stochforcvec
1101 c Compute friction and stochastic forces
1110 time_fric=time_fric+MPI_Wtime()-time00
1113 time_fric=time_fric+tcpu()-time00
1116 call stochastic_force(stochforcvec)
1118 time_stoch=time_stoch+MPI_Wtime()-time00
1120 time_stoch=time_stoch+tcpu()-time00
1123 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1124 c forces (d_as_work)
1126 call ginv_mult(fric_work, d_af_work)
1127 call ginv_mult(stochforcvec, d_as_work)
1130 c---------------------------------------------------------------------
1131 subroutine sddir_verlet1
1132 c Applying velocity Verlet algorithm - step 1 to velocities
1133 implicit real*8 (a-h,o-z)
1134 include 'DIMENSIONS'
1135 include 'COMMON.CONTROL'
1136 include 'COMMON.VAR'
1139 include 'COMMON.LANGEVIN'
1141 include 'COMMON.LANGEVIN.lang0'
1143 include 'COMMON.CHAIN'
1144 include 'COMMON.DERIV'
1145 include 'COMMON.GEO'
1146 include 'COMMON.LOCAL'
1147 include 'COMMON.INTERACT'
1148 include 'COMMON.IOUNITS'
1149 include 'COMMON.NAMES'
1150 c Revised 3/31/05 AL: correlation between random contributions to
1151 c position and velocity increments included.
1152 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1153 double precision adt,adt2
1155 c Add the contribution from BOTH friction and stochastic force to the
1156 c coordinates, but ONLY the contribution from the friction forces to velocities
1159 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1160 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1161 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1162 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1163 d_t(j,0)=d_t_old(j,0)+adt
1168 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1169 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1170 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1171 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1172 d_t(j,i)=d_t_old(j,i)+adt
1177 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1180 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1181 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1182 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1183 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1184 d_t(j,inres)=d_t_old(j,inres)+adt
1191 c---------------------------------------------------------------------
1192 subroutine sddir_verlet2
1193 c Calculating the adjusted velocities for accelerations
1194 implicit real*8 (a-h,o-z)
1195 include 'DIMENSIONS'
1196 include 'COMMON.CONTROL'
1197 include 'COMMON.VAR'
1200 include 'COMMON.LANGEVIN'
1202 include 'COMMON.LANGEVIN.lang0'
1204 include 'COMMON.CHAIN'
1205 include 'COMMON.DERIV'
1206 include 'COMMON.GEO'
1207 include 'COMMON.LOCAL'
1208 include 'COMMON.INTERACT'
1209 include 'COMMON.IOUNITS'
1210 include 'COMMON.NAMES'
1211 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1212 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1213 c Revised 3/31/05 AL: correlation between random contributions to
1214 c position and velocity increments included.
1215 c The correlation coefficients are calculated at low-friction limit.
1216 c Also, friction forces are now not calculated with new velocities.
1218 c call friction_force
1219 call stochastic_force(stochforcvec)
1221 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1222 c forces (d_as_work)
1224 call ginv_mult(stochforcvec, d_as_work1)
1230 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1231 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1236 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1237 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1242 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1245 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1246 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1247 & +cos60*d_as_work1(ind+j))*d_time
1254 c---------------------------------------------------------------------
1255 subroutine max_accel
1257 c Find the maximum difference in the accelerations of the the sites
1258 c at the beginning and the end of the time step.
1260 implicit real*8 (a-h,o-z)
1261 include 'DIMENSIONS'
1262 include 'COMMON.CONTROL'
1263 include 'COMMON.VAR'
1265 include 'COMMON.CHAIN'
1266 include 'COMMON.DERIV'
1267 include 'COMMON.GEO'
1268 include 'COMMON.LOCAL'
1269 include 'COMMON.INTERACT'
1270 include 'COMMON.IOUNITS'
1271 double precision aux(3),accel(3),accel_old(3),dacc
1273 c aux(j)=d_a(j,0)-d_a_old(j,0)
1274 accel_old(j)=d_a_old(j,0)
1281 c 7/3/08 changed to asymmetric difference
1283 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1284 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1285 accel(j)=accel(j)+0.5d0*d_a(j,i)
1286 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1287 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1288 dacc=dabs(accel(j)-accel_old(j))
1289 c write (iout,*) i,dacc
1290 if (dacc.gt.amax) amax=dacc
1298 accel_old(j)=d_a_old(j,0)
1303 accel_old(j)=accel_old(j)+d_a_old(j,1)
1304 accel(j)=accel(j)+d_a(j,1)
1308 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1310 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1311 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1312 accel(j)=accel(j)+d_a(j,i+nres)
1316 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1317 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1318 dacc=dabs(accel(j)-accel_old(j))
1319 c write (iout,*) "side-chain",i,dacc
1320 if (dacc.gt.amax) amax=dacc
1324 accel_old(j)=accel_old(j)+d_a_old(j,i)
1325 accel(j)=accel(j)+d_a(j,i)
1326 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1331 c---------------------------------------------------------------------
1332 subroutine predict_edrift(epdrift)
1334 c Predict the drift of the potential energy
1336 implicit real*8 (a-h,o-z)
1337 include 'DIMENSIONS'
1338 include 'COMMON.CONTROL'
1339 include 'COMMON.VAR'
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.MUCA'
1348 double precision epdrift,epdriftij
1349 c Drift of the potential energy
1355 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1356 if (lmuca) epdriftij=epdriftij*factor
1357 c write (iout,*) "back",i,j,epdriftij
1358 if (epdriftij.gt.epdrift) epdrift=epdriftij
1362 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1365 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1366 if (lmuca) epdriftij=epdriftij*factor
1367 c write (iout,*) "side",i,j,epdriftij
1368 if (epdriftij.gt.epdrift) epdrift=epdriftij
1372 epdrift=0.5d0*epdrift*d_time*d_time
1373 c write (iout,*) "epdrift",epdrift
1376 c-----------------------------------------------------------------------
1377 subroutine verlet_bath
1379 c Coupling to the thermostat by using the Berendsen algorithm
1381 implicit real*8 (a-h,o-z)
1382 include 'DIMENSIONS'
1383 include 'COMMON.CONTROL'
1384 include 'COMMON.VAR'
1386 include 'COMMON.CHAIN'
1387 include 'COMMON.DERIV'
1388 include 'COMMON.GEO'
1389 include 'COMMON.LOCAL'
1390 include 'COMMON.INTERACT'
1391 include 'COMMON.IOUNITS'
1392 include 'COMMON.NAMES'
1393 double precision T_half,fact
1395 T_half=2.0d0/(dimen3*Rb)*EK
1396 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1397 c write(iout,*) "T_half", T_half
1398 c write(iout,*) "EK", EK
1399 c write(iout,*) "fact", fact
1401 d_t(j,0)=fact*d_t(j,0)
1405 d_t(j,i)=fact*d_t(j,i)
1409 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1412 d_t(j,inres)=fact*d_t(j,inres)
1418 c---------------------------------------------------------
1420 c Set up the initial conditions of a MD simulation
1421 implicit real*8 (a-h,o-z)
1422 include 'DIMENSIONS'
1426 integer IERROR,ERRCODE
1428 include 'COMMON.SETUP'
1429 include 'COMMON.CONTROL'
1430 include 'COMMON.VAR'
1433 include 'COMMON.LANGEVIN'
1435 include 'COMMON.LANGEVIN.lang0'
1437 include 'COMMON.CHAIN'
1438 include 'COMMON.DERIV'
1439 include 'COMMON.GEO'
1440 include 'COMMON.LOCAL'
1441 include 'COMMON.INTERACT'
1442 include 'COMMON.IOUNITS'
1443 include 'COMMON.NAMES'
1444 include 'COMMON.REMD'
1445 real*8 energia_long(0:n_ene),
1446 & energia_short(0:n_ene),vcm(3),incr(3)
1447 double precision cm(3),L(3),xv,sigv,lowb,highb
1448 double precision varia(maxvar)
1456 c write(iout,*) "d_time", d_time
1457 c Compute the standard deviations of stochastic forces for Langevin dynamics
1458 c if the friction coefficients do not depend on surface area
1459 if (lang.gt.0 .and. .not.surfarea) then
1461 stdforcp(i)=stdfp*dsqrt(gamp)
1464 stdforcsc(i)=stdfsc(iabs(itype(i)))
1465 & *dsqrt(gamsc(iabs(itype(i))))
1468 c Open the pdb file for snapshotshots
1471 if (ilen(tmpdir).gt.0)
1472 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1473 & liczba(:ilen(liczba))//".pdb")
1475 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1479 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1480 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1481 & liczba(:ilen(liczba))//".x")
1482 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1485 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1486 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1487 & liczba(:ilen(liczba))//".cx")
1488 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1494 if (ilen(tmpdir).gt.0)
1495 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1496 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1498 if (ilen(tmpdir).gt.0)
1499 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1500 cartname=prefix(:ilen(prefix))//"_MD.cx"
1504 write (qstr,'(256(1h ))')
1507 iq = qinfrag(i,iset)*10
1508 iw = wfrag(i,iset)/100
1510 if(me.eq.king.or..not.out1file)
1511 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1512 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1517 iq = qinpair(i,iset)*10
1518 iw = wpair(i,iset)/100
1520 if(me.eq.king.or..not.out1file)
1521 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1522 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1526 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1528 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1530 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1532 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1536 if (restart1file) then
1538 & inquire(file=mremd_rst_name,exist=file_exist)
1540 write (*,*) me," Before broadcast: file_exist",file_exist
1541 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1543 write (*,*) me," After broadcast: file_exist",file_exist
1544 c inquire(file=mremd_rst_name,exist=file_exist)
1546 if(me.eq.king.or..not.out1file)
1547 & write(iout,*) "Initial state read by master and distributed"
1549 if (ilen(tmpdir).gt.0)
1550 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1551 & //liczba(:ilen(liczba))//'.rst')
1552 inquire(file=rest2name,exist=file_exist)
1555 if(.not.restart1file) then
1556 if(me.eq.king.or..not.out1file)
1557 & write(iout,*) "Initial state will be read from file ",
1558 & rest2name(:ilen(rest2name))
1561 call rescale_weights(t_bath)
1563 if(me.eq.king.or..not.out1file)then
1564 if (restart1file) then
1565 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1568 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1571 write(iout,*) "Initial velocities randomly generated"
1578 c Generate initial velocities
1579 if(me.eq.king.or..not.out1file)
1580 & write(iout,*) "Initial velocities randomly generated"
1583 CtotTafm is the variable for AFM time which eclipsed during
1586 c rest2name = prefix(:ilen(prefix))//'.rst'
1587 if(me.eq.king.or..not.out1file)then
1588 write (iout,*) "Initial velocities"
1590 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1591 & (d_t(j,i+nres),j=1,3)
1593 c Zeroing the total angular momentum of the system
1594 write(iout,*) "Calling the zero-angular
1595 & momentum subroutine"
1598 c Getting the potential energy and forces and velocities and accelerations
1600 c write (iout,*) "velocity of the center of the mass:"
1601 c write (iout,*) (vcm(j),j=1,3)
1603 d_t(j,0)=d_t(j,0)-vcm(j)
1605 c Removing the velocity of the center of mass
1607 if(me.eq.king.or..not.out1file)then
1608 write (iout,*) "vcm right after adjustment:"
1609 write (iout,*) (vcm(j),j=1,3)
1612 if (start_from_model) then
1613 i_model=iran_num(1,constr_homology)
1614 write (iout,*) 'starting from model ',i_model
1617 c(j,i)=chomo(j,i,i_model)
1621 call int_from_cart(.true.,.false.)
1622 call sc_loc_geom(.false.)
1625 dc(j,i)=c(j,i+1)-c(j,i)
1626 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1631 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1632 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1636 c call chainbuild_cart
1637 write (iout,*) "PREMINIM ",preminim
1638 if(iranconf.eq.0 .and. preminim) then
1640 write (iout,*) 'Calling OVERLAP_SC'
1641 call overlap_sc(fail)
1645 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1646 print *,'SC_move',nft_sc,etot
1647 if(me.eq.king.or..not.out1file)
1648 & write(iout,*) 'SC_move',nft_sc,etot
1652 print *, 'Calling MINIM_DC'
1653 call minim_dc(etot,iretcode,nfun)
1655 call geom_to_var(nvar,varia)
1656 print *,'Calling MINIMIZE.'
1657 call minimize(etot,varia,iretcode,nfun)
1658 call var_to_geom(nvar,varia)
1660 if(me.eq.king.or..not.out1file) then
1661 write(iout,*) "Minimized energy is",etot
1662 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1663 call etotal(potEcomp)
1664 call enerprint(potEcomp)
1667 if(iranconf.ne.0) then
1668 c 8/22/17 AL Loop to produce a low-energy random conformation
1671 if(me.eq.king.or..not.out1file)
1672 & write (iout,*) 'Calling OVERLAP_SC'
1673 call overlap_sc(fail)
1677 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1678 print *,'SC_move',nft_sc,etot
1679 if(me.eq.king.or..not.out1file)
1680 & write(iout,*) 'SC_move',nft_sc,etot
1684 print *, 'Calling MINIM_DC'
1685 call minim_dc(etot,iretcode,nfun)
1687 call geom_to_var(nvar,varia)
1688 print *,'Calling MINIMIZE.'
1689 call minimize(etot,varia,iretcode,nfun)
1690 call var_to_geom(nvar,varia)
1692 if(me.eq.king.or..not.out1file)
1693 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1695 if (isnan(etot) .or. etot.gt.1.0d4) then
1696 write (iout,*) "Energy too large",etot,
1697 & " trying another random conformation"
1700 call gen_rand_conf(itmp,*30)
1702 30 write (iout,*) 'Failed to generate random conformation',
1703 & ', itrial=',itrial
1704 write (*,*) 'Processor:',me,
1705 & ' Failed to generate random conformation',
1714 write (iout,'(a,i3,a)') 'Processor:',me,
1715 & ' error in generating random conformation.'
1716 write (*,'(a,i3,a)') 'Processor:',me,
1717 & ' error in generating random conformation.'
1720 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1730 write (iout,'(a,i3,a)') 'Processor:',me,
1731 & ' failed to generate a low-energy random conformation.'
1732 write (*,'(a,i3,a)') 'Processor:',me,
1733 & ' failed to generate a low-energy random conformation.'
1736 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1743 call chainbuild_cart
1748 kinetic_T=2.0d0/(dimen3*Rb)*EK
1749 if(me.eq.king.or..not.out1file)then
1759 call etotal(potEcomp)
1760 if (large) call enerprint(potEcomp)
1763 t_etotal=t_etotal+MPI_Wtime()-tt0
1765 t_etotal=t_etotal+tcpu()-tt0
1772 if (amax*d_time .gt. dvmax) then
1773 d_time=d_time*dvmax/amax
1774 if(me.eq.king.or..not.out1file) write (iout,*)
1775 & "Time step reduced to",d_time,
1776 & " because of too large initial acceleration."
1778 if(me.eq.king.or..not.out1file)then
1779 write(iout,*) "Potential energy and its components"
1780 call enerprint(potEcomp)
1781 c write(iout,*) (potEcomp(i),i=0,n_ene)
1783 potE=potEcomp(0)-potEcomp(20)
1786 if (ntwe.ne.0) call statout(itime)
1787 if(me.eq.king.or..not.out1file)
1788 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1789 & " Kinetic energy",EK," Potential energy",potE,
1790 & " Total energy",totE," Maximum acceleration ",
1793 write (iout,*) "Initial coordinates"
1795 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1796 & (c(j,i+nres),j=1,3)
1798 write (iout,*) "Initial dC"
1800 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1801 & (dc(j,i+nres),j=1,3)
1803 write (iout,*) "Initial velocities"
1804 write (iout,"(13x,' backbone ',23x,' side chain')")
1806 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1807 & (d_t(j,i+nres),j=1,3)
1809 write (iout,*) "Initial accelerations"
1811 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1812 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1813 & (d_a(j,i+nres),j=1,3)
1819 d_t_old(j,i)=d_t(j,i)
1820 d_a_old(j,i)=d_a(j,i)
1822 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1831 call etotal_short(energia_short)
1832 if (large) call enerprint(potEcomp)
1835 t_eshort=t_eshort+MPI_Wtime()-tt0
1837 t_eshort=t_eshort+tcpu()-tt0
1842 if(.not.out1file .and. large) then
1843 write (iout,*) "energia_long",energia_long(0),
1844 & " energia_short",energia_short(0),
1845 & " total",energia_long(0)+energia_short(0)
1846 write (iout,*) "Initial fast-force accelerations"
1848 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1849 & (d_a(j,i+nres),j=1,3)
1852 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1855 d_a_short(j,i)=d_a(j,i)
1864 call etotal_long(energia_long)
1865 if (large) call enerprint(potEcomp)
1868 t_elong=t_elong+MPI_Wtime()-tt0
1870 t_elong=t_elong+tcpu()-tt0
1875 if(.not.out1file .and. large) then
1876 write (iout,*) "energia_long",energia_long(0)
1877 write (iout,*) "Initial slow-force accelerations"
1879 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1880 & (d_a(j,i+nres),j=1,3)
1884 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1886 t_enegrad=t_enegrad+tcpu()-tt0
1891 c-----------------------------------------------------------
1892 subroutine random_vel
1893 implicit real*8 (a-h,o-z)
1894 include 'DIMENSIONS'
1895 include 'COMMON.CONTROL'
1896 include 'COMMON.VAR'
1899 include 'COMMON.LANGEVIN'
1901 include 'COMMON.LANGEVIN.lang0'
1903 include 'COMMON.CHAIN'
1904 include 'COMMON.DERIV'
1905 include 'COMMON.GEO'
1906 include 'COMMON.LOCAL'
1907 include 'COMMON.INTERACT'
1908 include 'COMMON.IOUNITS'
1909 include 'COMMON.NAMES'
1910 include 'COMMON.TIME1'
1911 double precision xv,sigv,lowb,highb,vec_afm(3)
1912 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1913 c First generate velocities in the eigenspace of the G matrix
1914 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1921 sigv=dsqrt((Rb*t_bath)/geigen(i))
1924 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1926 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1927 c & " d_t_work_new",d_t_work_new(ii)
1930 C if (SELFGUIDE.gt.0) then
1933 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
1934 C distance=distance+vec_afm(j)**2
1936 C distance=dsqrt(distance)
1938 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
1939 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
1940 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
1941 C & d_t_work_new(j+(afmend-1)*3)
1952 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1955 c write (iout,*) "Ek from eigenvectors",Ek1
1957 c Transform velocities to UNRES coordinate space
1963 d_t_work(ind)=d_t_work(ind)
1964 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1966 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1970 c Transfer to the d_t vector
1972 d_t(j,0)=d_t_work(j)
1978 d_t(j,i)=d_t_work(ind)
1982 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1985 d_t(j,i+nres)=d_t_work(ind)
1990 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1991 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1996 c-----------------------------------------------------------
1997 subroutine sd_verlet_p_setup
1998 c Sets up the parameters of stochastic Verlet algorithm
1999 implicit real*8 (a-h,o-z)
2000 include 'DIMENSIONS'
2004 include 'COMMON.CONTROL'
2005 include 'COMMON.VAR'
2008 include 'COMMON.LANGEVIN'
2010 include 'COMMON.LANGEVIN.lang0'
2012 include 'COMMON.CHAIN'
2013 include 'COMMON.DERIV'
2014 include 'COMMON.GEO'
2015 include 'COMMON.LOCAL'
2016 include 'COMMON.INTERACT'
2017 include 'COMMON.IOUNITS'
2018 include 'COMMON.NAMES'
2019 include 'COMMON.TIME1'
2020 double precision emgdt(MAXRES6),
2021 & pterm,vterm,rho,rhoc,vsig,
2022 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2023 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2024 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2025 logical lprn /.false./
2026 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2027 double precision ktm
2034 c AL 8/17/04 Code adapted from tinker
2036 c Get the frictional and random terms for stochastic dynamics in the
2037 c eigenspace of mass-scaled UNRES friction matrix
2040 gdt = fricgam(i) * d_time
2042 c Stochastic dynamics reduces to simple MD for zero friction
2044 if (gdt .le. zero) then
2045 pfric_vec(i) = 1.0d0
2046 vfric_vec(i) = d_time
2047 afric_vec(i) = 0.5d0 * d_time * d_time
2048 prand_vec(i) = 0.0d0
2049 vrand_vec1(i) = 0.0d0
2050 vrand_vec2(i) = 0.0d0
2052 c Analytical expressions when friction coefficient is large
2055 if (gdt .ge. gdt_radius) then
2058 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2059 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2060 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2061 vterm = 1.0d0 - egdt**2
2062 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2064 c Use series expansions when friction coefficient is small
2075 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2076 & - gdt5/120.0d0 + gdt6/720.0d0
2077 & - gdt7/5040.0d0 + gdt8/40320.0d0
2078 & - gdt9/362880.0d0) / fricgam(i)**2
2079 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2080 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2081 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2082 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2083 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2084 & + 127.0d0*gdt9/90720.0d0
2085 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2086 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2087 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2088 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2089 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2090 & - 17.0d0*gdt2/1280.0d0
2091 & + 17.0d0*gdt3/6144.0d0
2092 & + 40967.0d0*gdt4/34406400.0d0
2093 & - 57203.0d0*gdt5/275251200.0d0
2094 & - 1429487.0d0*gdt6/13212057600.0d0)
2097 c Compute the scaling factors of random terms for the nonzero friction case
2099 ktm = 0.5d0*d_time/fricgam(i)
2100 psig = dsqrt(ktm*pterm) / fricgam(i)
2101 vsig = dsqrt(ktm*vterm)
2102 rhoc = dsqrt(1.0d0 - rho*rho)
2104 vrand_vec1(i) = vsig * rho
2105 vrand_vec2(i) = vsig * rhoc
2110 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2113 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2114 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2118 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2121 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2122 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2123 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2124 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2125 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2126 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2129 t_sdsetup=t_sdsetup+MPI_Wtime()
2131 t_sdsetup=t_sdsetup+tcpu()-tt0
2135 c-------------------------------------------------------------
2136 subroutine eigtransf1(n,ndim,ab,d,c)
2139 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2145 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2151 c-------------------------------------------------------------
2152 subroutine eigtransf(n,ndim,a,b,d,c)
2155 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2161 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2167 c-------------------------------------------------------------
2168 subroutine sd_verlet1
2169 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2170 implicit real*8 (a-h,o-z)
2171 include 'DIMENSIONS'
2172 include 'COMMON.CONTROL'
2173 include 'COMMON.VAR'
2176 include 'COMMON.LANGEVIN'
2178 include 'COMMON.LANGEVIN.lang0'
2180 include 'COMMON.CHAIN'
2181 include 'COMMON.DERIV'
2182 include 'COMMON.GEO'
2183 include 'COMMON.LOCAL'
2184 include 'COMMON.INTERACT'
2185 include 'COMMON.IOUNITS'
2186 include 'COMMON.NAMES'
2187 double precision stochforcvec(MAXRES6)
2188 common /stochcalc/ stochforcvec
2189 logical lprn /.false./
2191 c write (iout,*) "dc_old"
2193 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2194 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2197 dc_work(j)=dc_old(j,0)
2198 d_t_work(j)=d_t_old(j,0)
2199 d_a_work(j)=d_a_old(j,0)
2204 dc_work(ind+j)=dc_old(j,i)
2205 d_t_work(ind+j)=d_t_old(j,i)
2206 d_a_work(ind+j)=d_a_old(j,i)
2211 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2213 dc_work(ind+j)=dc_old(j,i+nres)
2214 d_t_work(ind+j)=d_t_old(j,i+nres)
2215 d_a_work(ind+j)=d_a_old(j,i+nres)
2223 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2227 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2228 & vfric_mat(i,j),afric_mat(i,j),
2229 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2237 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2238 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2239 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2240 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2242 d_t_work_new(i)=ddt1+0.5d0*ddt2
2243 d_t_work(i)=ddt1+ddt2
2248 d_t(j,0)=d_t_work(j)
2253 dc(j,i)=dc_work(ind+j)
2254 d_t(j,i)=d_t_work(ind+j)
2259 if (itype(i).ne.10) then
2262 dc(j,inres)=dc_work(ind+j)
2263 d_t(j,inres)=d_t_work(ind+j)
2270 c--------------------------------------------------------------------------
2271 subroutine sd_verlet2
2272 c Calculating the adjusted velocities for accelerations
2273 implicit real*8 (a-h,o-z)
2274 include 'DIMENSIONS'
2275 include 'COMMON.CONTROL'
2276 include 'COMMON.VAR'
2279 include 'COMMON.LANGEVIN'
2281 include 'COMMON.LANGEVIN.lang0'
2283 include 'COMMON.CHAIN'
2284 include 'COMMON.DERIV'
2285 include 'COMMON.GEO'
2286 include 'COMMON.LOCAL'
2287 include 'COMMON.INTERACT'
2288 include 'COMMON.IOUNITS'
2289 include 'COMMON.NAMES'
2290 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2291 common /stochcalc/ stochforcvec
2293 c Compute the stochastic forces which contribute to velocity change
2295 call stochastic_force(stochforcvecV)
2302 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2303 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2304 & vrand_mat2(i,j)*stochforcvecV(j)
2306 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2310 d_t(j,0)=d_t_work(j)
2315 d_t(j,i)=d_t_work(ind+j)
2320 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2323 d_t(j,inres)=d_t_work(ind+j)
2330 c-----------------------------------------------------------
2331 subroutine sd_verlet_ciccotti_setup
2332 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2334 implicit real*8 (a-h,o-z)
2335 include 'DIMENSIONS'
2339 include 'COMMON.CONTROL'
2340 include 'COMMON.VAR'
2343 include 'COMMON.LANGEVIN'
2345 include 'COMMON.LANGEVIN.lang0'
2347 include 'COMMON.CHAIN'
2348 include 'COMMON.DERIV'
2349 include 'COMMON.GEO'
2350 include 'COMMON.LOCAL'
2351 include 'COMMON.INTERACT'
2352 include 'COMMON.IOUNITS'
2353 include 'COMMON.NAMES'
2354 include 'COMMON.TIME1'
2355 double precision emgdt(MAXRES6),
2356 & pterm,vterm,rho,rhoc,vsig,
2357 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2358 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2359 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2360 logical lprn /.false./
2361 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2362 double precision ktm
2369 c AL 8/17/04 Code adapted from tinker
2371 c Get the frictional and random terms for stochastic dynamics in the
2372 c eigenspace of mass-scaled UNRES friction matrix
2375 write (iout,*) "i",i," fricgam",fricgam(i)
2376 gdt = fricgam(i) * d_time
2378 c Stochastic dynamics reduces to simple MD for zero friction
2380 if (gdt .le. zero) then
2381 pfric_vec(i) = 1.0d0
2382 vfric_vec(i) = d_time
2383 afric_vec(i) = 0.5d0*d_time*d_time
2384 prand_vec(i) = afric_vec(i)
2385 vrand_vec2(i) = vfric_vec(i)
2387 c Analytical expressions when friction coefficient is large
2392 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2393 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2394 prand_vec(i) = afric_vec(i)
2395 vrand_vec2(i) = vfric_vec(i)
2397 c Compute the scaling factors of random terms for the nonzero friction case
2399 c ktm = 0.5d0*d_time/fricgam(i)
2400 c psig = dsqrt(ktm*pterm) / fricgam(i)
2401 c vsig = dsqrt(ktm*vterm)
2402 c prand_vec(i) = psig*afric_vec(i)
2403 c vrand_vec2(i) = vsig*vfric_vec(i)
2408 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2411 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2412 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2416 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2418 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2419 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2420 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2421 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2422 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2424 t_sdsetup=t_sdsetup+MPI_Wtime()
2426 t_sdsetup=t_sdsetup+tcpu()-tt0
2430 c-------------------------------------------------------------
2431 subroutine sd_verlet1_ciccotti
2432 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2433 implicit real*8 (a-h,o-z)
2434 include 'DIMENSIONS'
2438 include 'COMMON.CONTROL'
2439 include 'COMMON.VAR'
2442 include 'COMMON.LANGEVIN'
2444 include 'COMMON.LANGEVIN.lang0'
2446 include 'COMMON.CHAIN'
2447 include 'COMMON.DERIV'
2448 include 'COMMON.GEO'
2449 include 'COMMON.LOCAL'
2450 include 'COMMON.INTERACT'
2451 include 'COMMON.IOUNITS'
2452 include 'COMMON.NAMES'
2453 double precision stochforcvec(MAXRES6)
2454 common /stochcalc/ stochforcvec
2455 logical lprn /.false./
2457 c write (iout,*) "dc_old"
2459 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2460 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2463 dc_work(j)=dc_old(j,0)
2464 d_t_work(j)=d_t_old(j,0)
2465 d_a_work(j)=d_a_old(j,0)
2470 dc_work(ind+j)=dc_old(j,i)
2471 d_t_work(ind+j)=d_t_old(j,i)
2472 d_a_work(ind+j)=d_a_old(j,i)
2477 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2479 dc_work(ind+j)=dc_old(j,i+nres)
2480 d_t_work(ind+j)=d_t_old(j,i+nres)
2481 d_a_work(ind+j)=d_a_old(j,i+nres)
2490 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2494 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2495 & vfric_mat(i,j),afric_mat(i,j),
2496 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2504 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2505 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2506 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2507 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2509 d_t_work_new(i)=ddt1+0.5d0*ddt2
2510 d_t_work(i)=ddt1+ddt2
2515 d_t(j,0)=d_t_work(j)
2520 dc(j,i)=dc_work(ind+j)
2521 d_t(j,i)=d_t_work(ind+j)
2526 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2529 dc(j,inres)=dc_work(ind+j)
2530 d_t(j,inres)=d_t_work(ind+j)
2537 c--------------------------------------------------------------------------
2538 subroutine sd_verlet2_ciccotti
2539 c Calculating the adjusted velocities for accelerations
2540 implicit real*8 (a-h,o-z)
2541 include 'DIMENSIONS'
2542 include 'COMMON.CONTROL'
2543 include 'COMMON.VAR'
2546 include 'COMMON.LANGEVIN'
2548 include 'COMMON.LANGEVIN.lang0'
2550 include 'COMMON.CHAIN'
2551 include 'COMMON.DERIV'
2552 include 'COMMON.GEO'
2553 include 'COMMON.LOCAL'
2554 include 'COMMON.INTERACT'
2555 include 'COMMON.IOUNITS'
2556 include 'COMMON.NAMES'
2557 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2558 common /stochcalc/ stochforcvec
2560 c Compute the stochastic forces which contribute to velocity change
2562 call stochastic_force(stochforcvecV)
2569 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2570 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2571 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2573 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2577 d_t(j,0)=d_t_work(j)
2582 d_t(j,i)=d_t_work(ind+j)
2587 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2590 d_t(j,inres)=d_t_work(ind+j)