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 c 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 ! AL 4/17/17: Reduce the steps if NaNs occurred.
406 if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
411 if (large.and. mod(itime,ntwe).eq.0)
412 & call enerprint(potEcomp)
415 t_etotal=t_etotal+MPI_Wtime()-tt0
417 t_etotal=t_etotal+tcpu()-tt0
420 potE=potEcomp(0)-potEcomp(20)
422 c Get the new accelerations
425 t_enegrad=t_enegrad+MPI_Wtime()-tt0
427 t_enegrad=t_enegrad+tcpu()-tt0
429 c Determine maximum acceleration and scale down the timestep if needed
431 amax=amax/(itime_scal+1)**2
432 call predict_edrift(epdrift)
433 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
434 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
436 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
438 itime_scal=itime_scal+ifac_time
439 c fac_time=dmin1(damax/amax,0.5d0)
440 fac_time=0.5d0**ifac_time
441 d_time=d_time*fac_time
442 if (lang.eq.2 .or. lang.eq.3) then
444 c write (iout,*) "Calling sd_verlet_setup: 1"
445 c Rescale the stochastic forces and recalculate or restore
446 c the matrices of tinker integrator
447 if (itime_scal.gt.maxflag_stoch) then
448 if (large) write (iout,'(a,i5,a)')
449 & "Calculate matrices for stochastic step;",
450 & " itime_scal ",itime_scal
452 call sd_verlet_p_setup
454 call sd_verlet_ciccotti_setup
456 write (iout,'(2a,i3,a,i3,1h.)')
457 & "Warning: cannot store matrices for stochastic",
458 & " integration because the index",itime_scal,
459 & " is greater than",maxflag_stoch
460 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
461 & " integration Langevin algorithm for better efficiency."
462 else if (flag_stoch(itime_scal)) then
463 if (large) write (iout,'(a,i5,a,l1)')
464 & "Restore matrices for stochastic step; itime_scal ",
465 & itime_scal," flag ",flag_stoch(itime_scal)
468 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
469 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
470 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
471 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
472 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
473 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
477 if (large) write (iout,'(2a,i5,a,l1)')
478 & "Calculate & store matrices for stochastic step;",
479 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
481 call sd_verlet_p_setup
483 call sd_verlet_ciccotti_setup
485 flag_stoch(ifac_time)=.true.
488 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
489 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
490 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
491 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
492 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
493 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
497 fac_time=1.0d0/dsqrt(fac_time)
499 stochforcvec(i)=fac_time*stochforcvec(i)
502 else if (lang.eq.1) then
503 c Rescale the accelerations due to stochastic forces
504 fac_time=1.0d0/dsqrt(fac_time)
506 d_as_work(i)=d_as_work(i)*fac_time
509 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
510 & "itime",itime," Timestep scaled down to ",
511 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
513 c Second step of the velocity Verlet algorithm
518 else if (lang.eq.3) then
520 call sd_verlet2_ciccotti
522 else if (lang.eq.1) then
527 if (rattle) call rattle2
530 C print *,totTafm,"TU?"
531 if (d_time.ne.d_time0) then
534 if (lang.eq.2 .or. lang.eq.3) then
535 if (large) write (iout,'(a)')
536 & "Restore original matrices for stochastic step"
537 c write (iout,*) "Calling sd_verlet_setup: 2"
538 c Restore the matrices of tinker integrator if the time step has been restored
541 pfric_mat(i,j)=pfric0_mat(i,j,0)
542 afric_mat(i,j)=afric0_mat(i,j,0)
543 vfric_mat(i,j)=vfric0_mat(i,j,0)
544 prand_mat(i,j)=prand0_mat(i,j,0)
545 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
546 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
555 c Calculate the kinetic and the total energy and the kinetic temperature
560 c write (iout,*) "step",itime," EK",EK," EK1",EK1
562 c Couple the system to Berendsen bath if needed
563 if (tbf .and. lang.eq.0) then
566 kinetic_T=2.0d0/(dimen3*Rb)*EK
567 c Backup the coordinates, velocities, and accelerations
571 d_t_old(j,i)=d_t(j,i)
572 d_a_old(j,i)=d_a(j,i)
576 if (mod(itime,ntwe).eq.0 .and. large) then
577 write (iout,*) "Velocities, step 2"
579 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
580 & (d_t(j,i+nres),j=1,3)
586 c-------------------------------------------------------------------------------
587 subroutine RESPA_step(itime)
588 c-------------------------------------------------------------------------------
589 c Perform a single RESPA step.
590 c-------------------------------------------------------------------------------
591 implicit real*8 (a-h,o-z)
595 integer IERROR,ERRCODE
597 include 'COMMON.SETUP'
598 include 'COMMON.CONTROL'
602 include 'COMMON.LANGEVIN'
604 include 'COMMON.LANGEVIN.lang0'
606 include 'COMMON.CHAIN'
607 include 'COMMON.DERIV'
609 include 'COMMON.LOCAL'
610 include 'COMMON.INTERACT'
611 include 'COMMON.IOUNITS'
612 include 'COMMON.NAMES'
613 include 'COMMON.TIME1'
615 common /shortcheck/ lprint_short
616 double precision energia_short(0:n_ene),
617 & energia_long(0:n_ene)
618 double precision cm(3),L(3),vcm(3),incr(3)
619 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
620 & d_a_old0(3,0:maxres2)
621 logical PRINT_AMTS_MSG /.false./
622 integer ilen,count,rstcount
625 integer maxcount_scale /10/
627 double precision stochforcvec(MAXRES6)
628 common /stochcalc/ stochforcvec
631 common /cipiszcze/ itt
634 if (large.and. mod(itime,ntwe).eq.0) then
635 write (iout,*) "***************** RESPA itime",itime
636 write (iout,*) "Cartesian and internal coordinates: step 0"
641 write (iout,*) "Accelerations from long-range forces"
643 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
644 & (d_a(j,i+nres),j=1,3)
646 write (iout,*) "Velocities, step 0"
648 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
649 & (d_t(j,i+nres),j=1,3)
654 c Perform the initial RESPA step (increment velocities)
655 c write (iout,*) "*********************** RESPA ini"
658 if (mod(itime,ntwe).eq.0 .and. large) then
659 write (iout,*) "Velocities, end"
661 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
662 & (d_t(j,i+nres),j=1,3)
666 c Compute the short-range forces
672 C 7/2/2009 commented out
674 c call etotal_short(energia_short)
677 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
680 d_a(j,i)=d_a_short(j,i)
684 if (large.and. mod(itime,ntwe).eq.0) then
685 write (iout,*) "energia_short",energia_short(0)
686 write (iout,*) "Accelerations from short-range forces"
688 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
689 & (d_a(j,i+nres),j=1,3)
694 t_enegrad=t_enegrad+MPI_Wtime()-tt0
696 t_enegrad=t_enegrad+tcpu()-tt0
701 d_t_old(j,i)=d_t(j,i)
702 d_a_old(j,i)=d_a(j,i)
705 c 6/30/08 A-MTS: attempt at increasing the split number
708 dc_old0(j,i)=dc_old(j,i)
709 d_t_old0(j,i)=d_t_old(j,i)
710 d_a_old0(j,i)=d_a_old(j,i)
713 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
714 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
721 c write (iout,*) "itime",itime," ntime_split",ntime_split
722 c Split the time step
723 d_time=d_time0/ntime_split
724 c Perform the short-range RESPA steps (velocity Verlet increments of
725 c positions and velocities using short-range forces)
726 c write (iout,*) "*********************** RESPA split"
727 do itsplit=1,ntime_split
730 else if (lang.eq.2 .or. lang.eq.3) then
732 call stochastic_force(stochforcvec)
735 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
737 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
742 c First step of the velocity Verlet algorithm
747 else if (lang.eq.3) then
749 call sd_verlet1_ciccotti
751 else if (lang.eq.1) then
756 c Build the chain from the newly calculated coordinates
758 if (rattle) call rattle1
760 if (large.and. mod(itime,ntwe).eq.0) then
761 write (iout,*) "***** ITSPLIT",itsplit
762 write (iout,*) "Cartesian and internal coordinates: step 1"
767 write (iout,*) "Velocities, step 1"
769 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
770 & (d_t(j,i+nres),j=1,3)
779 c Calculate energy and forces
780 c if (large.and. mod(itime,ntwe).eq.0) lprint_short=.true.
782 call etotal_short(energia_short)
783 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
784 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) )
786 if (PRINT_AMTS_MSG) write (iout,*)
787 & "Infinities/NaNs in energia_short",
788 & energia_short(0),"; increasing ntime_split to",ntime_split
789 ntime_split=ntime_split*2
790 if (ntime_split.gt.maxtime_split) then
792 write (iout,*) "Cannot rescue the run; aborting job.",
793 & " Retry with a smaller time step"
795 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
797 write (iout,*) "Cannot rescue the run; terminating.",
798 & " Retry with a smaller time step"
804 if (large.and. mod(itime,ntwe).eq.0)
805 & call enerprint(energia_short)
808 t_eshort=t_eshort+MPI_Wtime()-tt0
810 t_eshort=t_eshort+tcpu()-tt0
814 c Get the new accelerations
816 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
819 d_a_short(j,i)=d_a(j,i)
823 if (large.and. mod(itime,ntwe).eq.0) then
824 write (iout,*)"energia_short",energia_short(0)
825 write (iout,*) "Accelerations from short-range forces"
827 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
828 & (d_a(j,i+nres),j=1,3)
833 c Determine maximum acceleration and scale down the timestep if needed
835 amax=amax/ntime_split**2
836 call predict_edrift(epdrift)
837 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
838 & write (iout,*) "amax",amax," damax",damax,
839 & " epdrift",epdrift," epdriftmax",epdriftmax
840 c Exit loop and try with increased split number if the change of
841 c acceleration is too big
842 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
843 if (ntime_split.lt.maxtime_split) then
845 ntime_split=ntime_split*2
848 dc_old(j,i)=dc_old0(j,i)
849 d_t_old(j,i)=d_t_old0(j,i)
850 d_a_old(j,i)=d_a_old0(j,i)
853 if (PRINT_AMTS_MSG) then
854 write (iout,*) "acceleration/energy drift too large",amax,
855 & epdrift," split increased to ",ntime_split," itime",itime,
861 & "Uh-hu. Bumpy landscape. Maximum splitting number",
863 & " already reached!!! Trying to carry on!"
867 t_enegrad=t_enegrad+MPI_Wtime()-tt0
869 t_enegrad=t_enegrad+tcpu()-tt0
871 c Second step of the velocity Verlet algorithm
876 else if (lang.eq.3) then
878 call sd_verlet2_ciccotti
880 else if (lang.eq.1) then
885 if (rattle) call rattle2
886 c Backup the coordinates, velocities, and accelerations
890 d_t_old(j,i)=d_t(j,i)
891 d_a_old(j,i)=d_a(j,i)
898 c Restore the time step
900 c Compute long-range forces
907 call etotal_long(energia_long)
908 ! AL 4/17/2017 Handling NaNs
909 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
912 & "Infinitied/NaNs in energia_long, Aborting MPI job."
913 call enerprint(energia_long)
915 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
917 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
918 call enerprint(energia_long)
923 c lprint_short=.false.
924 if (large.and. mod(itime,ntwe).eq.0)
925 & call enerprint(energia_long)
928 t_elong=t_elong+MPI_Wtime()-tt0
930 t_elong=t_elong+tcpu()-tt0
936 t_enegrad=t_enegrad+MPI_Wtime()-tt0
938 t_enegrad=t_enegrad+tcpu()-tt0
940 c Compute accelerations from long-range forces
942 if (large.and. mod(itime,ntwe).eq.0) then
943 write (iout,*) "energia_long",energia_long(0)
944 write (iout,*) "Cartesian and internal coordinates: step 2"
949 write (iout,*) "Accelerations from long-range forces"
951 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
952 & (d_a(j,i+nres),j=1,3)
954 write (iout,*) "Velocities, step 2"
956 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
957 & (d_t(j,i+nres),j=1,3)
961 c Compute the final RESPA step (increment velocities)
962 c write (iout,*) "*********************** RESPA fin"
964 c Compute the complete potential energy
966 potEcomp(i)=energia_short(i)+energia_long(i)
968 potE=potEcomp(0)-potEcomp(20)
970 if (large.and. mod(itime,ntwe).eq.0) then
971 call enerprint(potEcomp)
972 write (iout,*) "potE",potD
975 c potE=energia_short(0)+energia_long(0)
978 c Calculate the kinetic and the total energy and the kinetic temperature
981 c Couple the system to Berendsen bath if needed
982 if (tbf .and. lang.eq.0) then
985 kinetic_T=2.0d0/(dimen3*Rb)*EK
986 c Backup the coordinates, velocities, and accelerations
988 if (mod(itime,ntwe).eq.0 .and. large) then
989 write (iout,*) "Velocities, end"
991 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
992 & (d_t(j,i+nres),j=1,3)
998 c---------------------------------------------------------------------
1000 c First and last RESPA step (incrementing velocities using long-range
1002 implicit real*8 (a-h,o-z)
1003 include 'DIMENSIONS'
1004 include 'COMMON.CONTROL'
1005 include 'COMMON.VAR'
1007 include 'COMMON.CHAIN'
1008 include 'COMMON.DERIV'
1009 include 'COMMON.GEO'
1010 include 'COMMON.LOCAL'
1011 include 'COMMON.INTERACT'
1012 include 'COMMON.IOUNITS'
1013 include 'COMMON.NAMES'
1015 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1019 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1023 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1026 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1032 c-----------------------------------------------------------------
1034 c Applying velocity Verlet algorithm - step 1 to coordinates
1035 implicit real*8 (a-h,o-z)
1036 include 'DIMENSIONS'
1037 include 'COMMON.CONTROL'
1038 include 'COMMON.VAR'
1040 include 'COMMON.CHAIN'
1041 include 'COMMON.DERIV'
1042 include 'COMMON.GEO'
1043 include 'COMMON.LOCAL'
1044 include 'COMMON.INTERACT'
1045 include 'COMMON.IOUNITS'
1046 include 'COMMON.NAMES'
1047 double precision adt,adt2
1050 write (iout,*) "VELVERLET1 START: DC"
1052 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1053 & (dc(j,i+nres),j=1,3)
1057 adt=d_a_old(j,0)*d_time
1059 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1060 d_t_new(j,0)=d_t_old(j,0)+adt2
1061 d_t(j,0)=d_t_old(j,0)+adt
1067 adt=d_a_old(j,i)*d_time
1069 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1070 d_t_new(j,i)=d_t_old(j,i)+adt2
1071 d_t(j,i)=d_t_old(j,i)+adt
1076 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1079 adt=d_a_old(j,inres)*d_time
1081 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1082 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1083 d_t(j,inres)=d_t_old(j,inres)+adt
1088 write (iout,*) "VELVERLET1 END: DC"
1090 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1091 & (dc(j,i+nres),j=1,3)
1096 c---------------------------------------------------------------------
1098 c Step 2 of the velocity Verlet algorithm: update velocities
1099 implicit real*8 (a-h,o-z)
1100 include 'DIMENSIONS'
1101 include 'COMMON.CONTROL'
1102 include 'COMMON.VAR'
1104 include 'COMMON.CHAIN'
1105 include 'COMMON.DERIV'
1106 include 'COMMON.GEO'
1107 include 'COMMON.LOCAL'
1108 include 'COMMON.INTERACT'
1109 include 'COMMON.IOUNITS'
1110 include 'COMMON.NAMES'
1112 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1116 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1120 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1123 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1129 c-----------------------------------------------------------------
1130 subroutine sddir_precalc
1131 c Applying velocity Verlet algorithm - step 1 to coordinates
1132 implicit real*8 (a-h,o-z)
1133 include 'DIMENSIONS'
1137 include 'COMMON.CONTROL'
1138 include 'COMMON.VAR'
1141 include 'COMMON.LANGEVIN'
1143 include 'COMMON.LANGEVIN.lang0'
1145 include 'COMMON.CHAIN'
1146 include 'COMMON.DERIV'
1147 include 'COMMON.GEO'
1148 include 'COMMON.LOCAL'
1149 include 'COMMON.INTERACT'
1150 include 'COMMON.IOUNITS'
1151 include 'COMMON.NAMES'
1152 include 'COMMON.TIME1'
1153 double precision stochforcvec(MAXRES6)
1154 common /stochcalc/ stochforcvec
1156 c Compute friction and stochastic forces
1165 time_fric=time_fric+MPI_Wtime()-time00
1168 time_fric=time_fric+tcpu()-time00
1171 call stochastic_force(stochforcvec)
1173 time_stoch=time_stoch+MPI_Wtime()-time00
1175 time_stoch=time_stoch+tcpu()-time00
1178 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1179 c forces (d_as_work)
1181 call ginv_mult(fric_work, d_af_work)
1182 call ginv_mult(stochforcvec, d_as_work)
1185 c---------------------------------------------------------------------
1186 subroutine sddir_verlet1
1187 c Applying velocity Verlet algorithm - step 1 to velocities
1188 implicit real*8 (a-h,o-z)
1189 include 'DIMENSIONS'
1190 include 'COMMON.CONTROL'
1191 include 'COMMON.VAR'
1194 include 'COMMON.LANGEVIN'
1196 include 'COMMON.LANGEVIN.lang0'
1198 include 'COMMON.CHAIN'
1199 include 'COMMON.DERIV'
1200 include 'COMMON.GEO'
1201 include 'COMMON.LOCAL'
1202 include 'COMMON.INTERACT'
1203 include 'COMMON.IOUNITS'
1204 include 'COMMON.NAMES'
1205 c Revised 3/31/05 AL: correlation between random contributions to
1206 c position and velocity increments included.
1207 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1208 double precision adt,adt2
1210 c Add the contribution from BOTH friction and stochastic force to the
1211 c coordinates, but ONLY the contribution from the friction forces to velocities
1214 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1215 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1216 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1217 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1218 d_t(j,0)=d_t_old(j,0)+adt
1223 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1224 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1225 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1226 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1227 d_t(j,i)=d_t_old(j,i)+adt
1232 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1235 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1236 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1237 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1238 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1239 d_t(j,inres)=d_t_old(j,inres)+adt
1246 c---------------------------------------------------------------------
1247 subroutine sddir_verlet2
1248 c Calculating the adjusted velocities for accelerations
1249 implicit real*8 (a-h,o-z)
1250 include 'DIMENSIONS'
1251 include 'COMMON.CONTROL'
1252 include 'COMMON.VAR'
1255 include 'COMMON.LANGEVIN'
1257 include 'COMMON.LANGEVIN.lang0'
1259 include 'COMMON.CHAIN'
1260 include 'COMMON.DERIV'
1261 include 'COMMON.GEO'
1262 include 'COMMON.LOCAL'
1263 include 'COMMON.INTERACT'
1264 include 'COMMON.IOUNITS'
1265 include 'COMMON.NAMES'
1266 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1267 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1268 c Revised 3/31/05 AL: correlation between random contributions to
1269 c position and velocity increments included.
1270 c The correlation coefficients are calculated at low-friction limit.
1271 c Also, friction forces are now not calculated with new velocities.
1273 c call friction_force
1274 call stochastic_force(stochforcvec)
1276 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1277 c forces (d_as_work)
1279 call ginv_mult(stochforcvec, d_as_work1)
1285 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1286 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1291 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1292 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1297 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1300 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1301 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1302 & +cos60*d_as_work1(ind+j))*d_time
1309 c---------------------------------------------------------------------
1310 subroutine max_accel
1312 c Find the maximum difference in the accelerations of the the sites
1313 c at the beginning and the end of the time step.
1315 implicit real*8 (a-h,o-z)
1316 include 'DIMENSIONS'
1317 include 'COMMON.CONTROL'
1318 include 'COMMON.VAR'
1320 include 'COMMON.CHAIN'
1321 include 'COMMON.DERIV'
1322 include 'COMMON.GEO'
1323 include 'COMMON.LOCAL'
1324 include 'COMMON.INTERACT'
1325 include 'COMMON.IOUNITS'
1326 double precision aux(3),accel(3),accel_old(3),dacc
1328 c aux(j)=d_a(j,0)-d_a_old(j,0)
1329 accel_old(j)=d_a_old(j,0)
1336 c 7/3/08 changed to asymmetric difference
1338 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1339 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1340 accel(j)=accel(j)+0.5d0*d_a(j,i)
1341 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1342 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1343 dacc=dabs(accel(j)-accel_old(j))
1344 c write (iout,*) i,dacc
1345 if (dacc.gt.amax) amax=dacc
1353 accel_old(j)=d_a_old(j,0)
1358 accel_old(j)=accel_old(j)+d_a_old(j,1)
1359 accel(j)=accel(j)+d_a(j,1)
1363 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1365 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1366 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1367 accel(j)=accel(j)+d_a(j,i+nres)
1371 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1372 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1373 dacc=dabs(accel(j)-accel_old(j))
1374 c write (iout,*) "side-chain",i,dacc
1375 if (dacc.gt.amax) amax=dacc
1379 accel_old(j)=accel_old(j)+d_a_old(j,i)
1380 accel(j)=accel(j)+d_a(j,i)
1381 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1386 c---------------------------------------------------------------------
1387 subroutine predict_edrift(epdrift)
1389 c Predict the drift of the potential energy
1391 implicit real*8 (a-h,o-z)
1392 include 'DIMENSIONS'
1393 include 'COMMON.CONTROL'
1394 include 'COMMON.VAR'
1396 include 'COMMON.CHAIN'
1397 include 'COMMON.DERIV'
1398 include 'COMMON.GEO'
1399 include 'COMMON.LOCAL'
1400 include 'COMMON.INTERACT'
1401 include 'COMMON.IOUNITS'
1402 include 'COMMON.MUCA'
1403 double precision epdrift,epdriftij
1404 c Drift of the potential energy
1410 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1411 if (lmuca) epdriftij=epdriftij*factor
1412 c write (iout,*) "back",i,j,epdriftij
1413 if (epdriftij.gt.epdrift) epdrift=epdriftij
1417 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1420 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1421 if (lmuca) epdriftij=epdriftij*factor
1422 c write (iout,*) "side",i,j,epdriftij
1423 if (epdriftij.gt.epdrift) epdrift=epdriftij
1427 epdrift=0.5d0*epdrift*d_time*d_time
1428 c write (iout,*) "epdrift",epdrift
1431 c-----------------------------------------------------------------------
1432 subroutine verlet_bath
1434 c Coupling to the thermostat by using the Berendsen algorithm
1436 implicit real*8 (a-h,o-z)
1437 include 'DIMENSIONS'
1438 include 'COMMON.CONTROL'
1439 include 'COMMON.VAR'
1441 include 'COMMON.CHAIN'
1442 include 'COMMON.DERIV'
1443 include 'COMMON.GEO'
1444 include 'COMMON.LOCAL'
1445 include 'COMMON.INTERACT'
1446 include 'COMMON.IOUNITS'
1447 include 'COMMON.NAMES'
1448 double precision T_half,fact
1450 T_half=2.0d0/(dimen3*Rb)*EK
1451 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1452 c write(iout,*) "T_half", T_half
1453 c write(iout,*) "EK", EK
1454 c write(iout,*) "fact", fact
1456 d_t(j,0)=fact*d_t(j,0)
1460 d_t(j,i)=fact*d_t(j,i)
1464 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1467 d_t(j,inres)=fact*d_t(j,inres)
1473 c---------------------------------------------------------
1475 c Set up the initial conditions of a MD simulation
1476 implicit real*8 (a-h,o-z)
1477 include 'DIMENSIONS'
1481 integer IERROR,ERRCODE
1483 include 'COMMON.SETUP'
1484 include 'COMMON.CONTROL'
1485 include 'COMMON.VAR'
1488 include 'COMMON.LANGEVIN'
1490 include 'COMMON.LANGEVIN.lang0'
1492 include 'COMMON.CHAIN'
1493 include 'COMMON.DERIV'
1494 include 'COMMON.GEO'
1495 include 'COMMON.LOCAL'
1496 include 'COMMON.INTERACT'
1497 include 'COMMON.IOUNITS'
1498 include 'COMMON.NAMES'
1499 include 'COMMON.REMD'
1500 real*8 energia_long(0:n_ene),
1501 & energia_short(0:n_ene),vcm(3),incr(3)
1502 double precision cm(3),L(3),xv,sigv,lowb,highb
1503 double precision varia(maxvar)
1511 c write(iout,*) "d_time", d_time
1512 c Compute the standard deviations of stochastic forces for Langevin dynamics
1513 c if the friction coefficients do not depend on surface area
1514 if (lang.gt.0 .and. .not.surfarea) then
1516 stdforcp(i)=stdfp*dsqrt(gamp)
1519 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1520 & *dsqrt(gamsc(iabs(itype(i))))
1523 c Open the pdb file for snapshotshots
1526 if (ilen(tmpdir).gt.0)
1527 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1528 & liczba(:ilen(liczba))//".pdb")
1530 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1534 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1535 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1536 & liczba(:ilen(liczba))//".x")
1537 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1540 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1541 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1542 & liczba(:ilen(liczba))//".cx")
1543 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1549 if (ilen(tmpdir).gt.0)
1550 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1551 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1553 if (ilen(tmpdir).gt.0)
1554 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1555 cartname=prefix(:ilen(prefix))//"_MD.cx"
1559 write (qstr,'(256(1h ))')
1562 iq = qinfrag(i,iset)*10
1563 iw = wfrag(i,iset)/100
1565 if(me.eq.king.or..not.out1file)
1566 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1567 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1572 iq = qinpair(i,iset)*10
1573 iw = wpair(i,iset)/100
1575 if(me.eq.king.or..not.out1file)
1576 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1577 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1581 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1583 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1585 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1587 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1591 if (restart1file) then
1593 & inquire(file=mremd_rst_name,exist=file_exist)
1595 write (*,*) me," Before broadcast: file_exist",file_exist
1596 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1598 write (*,*) me," After broadcast: file_exist",file_exist
1599 c inquire(file=mremd_rst_name,exist=file_exist)
1601 if(me.eq.king.or..not.out1file)
1602 & write(iout,*) "Initial state read by master and distributed"
1604 if (ilen(tmpdir).gt.0)
1605 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1606 & //liczba(:ilen(liczba))//'.rst')
1607 inquire(file=rest2name,exist=file_exist)
1610 if(.not.restart1file) then
1611 if(me.eq.king.or..not.out1file)
1612 & write(iout,*) "Initial state will be read from file ",
1613 & rest2name(:ilen(rest2name))
1616 call rescale_weights(t_bath)
1619 if(me.eq.king.or..not.out1file)then
1620 if (restart1file) then
1621 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1624 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1627 write(iout,*) "Initial velocities randomly generated"
1634 c Generate initial velocities
1635 if(me.eq.king.or..not.out1file)
1636 & write(iout,*) "Initial velocities randomly generated"
1639 CtotTafm is the variable for AFM time which eclipsed during
1642 c rest2name = prefix(:ilen(prefix))//'.rst'
1643 if(me.eq.king.or..not.out1file)then
1644 write (iout,*) "Initial velocities"
1646 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1647 & (d_t(j,i+nres),j=1,3)
1649 c Zeroing the total angular momentum of the system
1650 write(iout,*) "Calling the zero-angular
1651 & momentum subroutine"
1654 c Getting the potential energy and forces and velocities and accelerations
1656 write (iout,*) "velocity of the center of the mass:"
1657 write (iout,*) (vcm(j),j=1,3)
1660 d_t(j,0)=d_t(j,0)-vcm(j)
1662 c Removing the velocity of the center of mass
1664 if(me.eq.king.or..not.out1file)then
1665 write (iout,*) "vcm right after adjustment:"
1666 write (iout,*) (vcm(j),j=1,3)
1670 if(iranconf.ne.0) then
1671 c 8/22/17 AL Loop to produce a low-energy random conformation
1675 print *, 'Calling OVERLAP_SC'
1676 call overlap_sc(fail)
1680 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1681 print *,'SC_move',nft_sc,etot
1682 if(me.eq.king.or..not.out1file)
1683 & write(iout,*) 'SC_move',nft_sc,etot
1687 print *, 'Calling MINIM_DC'
1688 call minim_dc(etot,iretcode,nfun)
1690 call geom_to_var(nvar,varia)
1691 print *,'Calling MINIMIZE.'
1692 call minimize(etot,varia,iretcode,nfun)
1693 call var_to_geom(nvar,varia)
1695 if(me.eq.king.or..not.out1file)
1696 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1697 if (isnan(etot) .or. etot.gt.1.0d4) then
1698 write (iout,*) "Energy too large",etot,
1699 & " trying another random conformation"
1702 call gen_rand_conf(itmp,*30)
1704 30 write (iout,*) 'Failed to generate random conformation',
1705 & ', itrial=',itrial
1706 write (*,*) 'Processor:',me,
1707 & ' Failed to generate random conformation',
1717 write (iout,'(a,i3,a)') 'Processor:',me,
1718 & ' error in generating random conformation.'
1719 write (*,'(a,i3,a)') 'Processor:',me,
1720 & ' error in generating random conformation.'
1723 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1732 write (iout,'(a,i3,a)') 'Processor:',me,
1733 & ' failed to generate a low-energy random conformation.'
1734 write (*,'(a,i3,a)') 'Processor:',me,
1735 & ' failed to generate a low-energy random conformation.'
1738 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1743 else if (indpdb.gt.0) then
1744 C 8/22/17 AL Minimize initial PDB structure
1746 print *, 'Calling MINIM_DC'
1747 call minim_dc(etot,iretcode,nfun)
1749 call geom_to_var(nvar,varia)
1750 print *,'Calling MINIMIZE.'
1751 call minimize(etot,varia,iretcode,nfun)
1752 call var_to_geom(nvar,varia)
1754 if(me.eq.king.or..not.out1file)
1755 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1758 call chainbuild_cart
1763 kinetic_T=2.0d0/(dimen3*Rb)*EK
1764 if(me.eq.king.or..not.out1file)then
1774 call etotal(potEcomp)
1775 if (large) call enerprint(potEcomp)
1778 t_etotal=t_etotal+MPI_Wtime()-tt0
1780 t_etotal=t_etotal+tcpu()-tt0
1787 if (amax*d_time .gt. dvmax) then
1788 d_time=d_time*dvmax/amax
1789 if(me.eq.king.or..not.out1file) write (iout,*)
1790 & "Time step reduced to",d_time,
1791 & " because of too large initial acceleration."
1793 if(me.eq.king.or..not.out1file)then
1794 write(iout,*) "Potential energy and its components"
1795 call enerprint(potEcomp)
1796 write(iout,*) (potEcomp(i),i=0,n_ene)
1798 potE=potEcomp(0)-potEcomp(20)
1801 if (ntwe.ne.0) call statout(itime)
1802 if(me.eq.king.or..not.out1file)
1803 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1804 & " Kinetic energy",EK," Potential energy",potE,
1805 & " Total energy",totE," Maximum acceleration ",
1808 write (iout,*) "Initial coordinates"
1810 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1811 & (c(j,i+nres),j=1,3)
1813 write (iout,*) "Initial dC"
1815 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1816 & (dc(j,i+nres),j=1,3)
1818 write (iout,*) "Initial velocities"
1819 write (iout,"(13x,' backbone ',23x,' side chain')")
1821 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1822 & (d_t(j,i+nres),j=1,3)
1824 write (iout,*) "Initial accelerations"
1826 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1827 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1828 & (d_a(j,i+nres),j=1,3)
1834 d_t_old(j,i)=d_t(j,i)
1835 d_a_old(j,i)=d_a(j,i)
1837 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1846 call etotal_short(energia_short)
1847 if (large) call enerprint(potEcomp)
1850 t_eshort=t_eshort+MPI_Wtime()-tt0
1852 t_eshort=t_eshort+tcpu()-tt0
1857 if(.not.out1file .and. large) then
1858 write (iout,*) "energia_long",energia_long(0),
1859 & " energia_short",energia_short(0),
1860 & " total",energia_long(0)+energia_short(0)
1861 write (iout,*) "Initial fast-force accelerations"
1863 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1864 & (d_a(j,i+nres),j=1,3)
1867 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1870 d_a_short(j,i)=d_a(j,i)
1879 call etotal_long(energia_long)
1880 if (large) call enerprint(potEcomp)
1883 t_elong=t_elong+MPI_Wtime()-tt0
1885 t_elong=t_elong+tcpu()-tt0
1890 if(.not.out1file .and. large) then
1891 write (iout,*) "energia_long",energia_long(0)
1892 write (iout,*) "Initial slow-force accelerations"
1894 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1895 & (d_a(j,i+nres),j=1,3)
1899 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1901 t_enegrad=t_enegrad+tcpu()-tt0
1906 c-----------------------------------------------------------
1907 subroutine random_vel
1908 implicit real*8 (a-h,o-z)
1909 include 'DIMENSIONS'
1910 include 'COMMON.CONTROL'
1911 include 'COMMON.VAR'
1914 include 'COMMON.LANGEVIN'
1916 include 'COMMON.LANGEVIN.lang0'
1918 include 'COMMON.CHAIN'
1919 include 'COMMON.DERIV'
1920 include 'COMMON.GEO'
1921 include 'COMMON.LOCAL'
1922 include 'COMMON.INTERACT'
1923 include 'COMMON.IOUNITS'
1924 include 'COMMON.NAMES'
1925 include 'COMMON.TIME1'
1926 double precision xv,sigv,lowb,highb,vec_afm(3)
1927 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1928 c First generate velocities in the eigenspace of the G matrix
1929 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1936 sigv=dsqrt((Rb*t_bath)/geigen(i))
1939 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1941 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1942 c & " d_t_work_new",d_t_work_new(ii)
1945 C if (SELFGUIDE.gt.0) then
1948 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
1949 C distance=distance+vec_afm(j)**2
1951 C distance=dsqrt(distance)
1953 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
1954 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
1955 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
1956 C & d_t_work_new(j+(afmend-1)*3)
1967 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1970 c write (iout,*) "Ek from eigenvectors",Ek1
1972 c Transform velocities to UNRES coordinate space
1978 d_t_work(ind)=d_t_work(ind)
1979 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1981 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1985 c Transfer to the d_t vector
1987 d_t(j,0)=d_t_work(j)
1993 d_t(j,i)=d_t_work(ind)
1997 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2000 d_t(j,i+nres)=d_t_work(ind)
2005 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2006 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2011 c-----------------------------------------------------------
2012 subroutine sd_verlet_p_setup
2013 c Sets up the parameters of stochastic Verlet algorithm
2014 implicit real*8 (a-h,o-z)
2015 include 'DIMENSIONS'
2019 include 'COMMON.CONTROL'
2020 include 'COMMON.VAR'
2023 include 'COMMON.LANGEVIN'
2025 include 'COMMON.LANGEVIN.lang0'
2027 include 'COMMON.CHAIN'
2028 include 'COMMON.DERIV'
2029 include 'COMMON.GEO'
2030 include 'COMMON.LOCAL'
2031 include 'COMMON.INTERACT'
2032 include 'COMMON.IOUNITS'
2033 include 'COMMON.NAMES'
2034 include 'COMMON.TIME1'
2035 double precision emgdt(MAXRES6),
2036 & pterm,vterm,rho,rhoc,vsig,
2037 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2038 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2039 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2040 logical lprn /.false./
2041 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2042 double precision ktm
2049 c AL 8/17/04 Code adapted from tinker
2051 c Get the frictional and random terms for stochastic dynamics in the
2052 c eigenspace of mass-scaled UNRES friction matrix
2055 gdt = fricgam(i) * d_time
2057 c Stochastic dynamics reduces to simple MD for zero friction
2059 if (gdt .le. zero) then
2060 pfric_vec(i) = 1.0d0
2061 vfric_vec(i) = d_time
2062 afric_vec(i) = 0.5d0 * d_time * d_time
2063 prand_vec(i) = 0.0d0
2064 vrand_vec1(i) = 0.0d0
2065 vrand_vec2(i) = 0.0d0
2067 c Analytical expressions when friction coefficient is large
2070 if (gdt .ge. gdt_radius) then
2073 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2074 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2075 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2076 vterm = 1.0d0 - egdt**2
2077 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2079 c Use series expansions when friction coefficient is small
2090 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2091 & - gdt5/120.0d0 + gdt6/720.0d0
2092 & - gdt7/5040.0d0 + gdt8/40320.0d0
2093 & - gdt9/362880.0d0) / fricgam(i)**2
2094 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2095 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2096 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2097 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2098 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2099 & + 127.0d0*gdt9/90720.0d0
2100 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2101 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2102 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2103 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2104 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2105 & - 17.0d0*gdt2/1280.0d0
2106 & + 17.0d0*gdt3/6144.0d0
2107 & + 40967.0d0*gdt4/34406400.0d0
2108 & - 57203.0d0*gdt5/275251200.0d0
2109 & - 1429487.0d0*gdt6/13212057600.0d0)
2112 c Compute the scaling factors of random terms for the nonzero friction case
2114 ktm = 0.5d0*d_time/fricgam(i)
2115 psig = dsqrt(ktm*pterm) / fricgam(i)
2116 vsig = dsqrt(ktm*vterm)
2117 rhoc = dsqrt(1.0d0 - rho*rho)
2119 vrand_vec1(i) = vsig * rho
2120 vrand_vec2(i) = vsig * rhoc
2125 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2128 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2129 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2133 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2136 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2137 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2138 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2139 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2140 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2141 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2144 t_sdsetup=t_sdsetup+MPI_Wtime()
2146 t_sdsetup=t_sdsetup+tcpu()-tt0
2150 c-------------------------------------------------------------
2151 subroutine eigtransf1(n,ndim,ab,d,c)
2154 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2160 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2166 c-------------------------------------------------------------
2167 subroutine eigtransf(n,ndim,a,b,d,c)
2170 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2176 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2182 c-------------------------------------------------------------
2183 subroutine sd_verlet1
2184 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2185 implicit real*8 (a-h,o-z)
2186 include 'DIMENSIONS'
2187 include 'COMMON.CONTROL'
2188 include 'COMMON.VAR'
2191 include 'COMMON.LANGEVIN'
2193 include 'COMMON.LANGEVIN.lang0'
2195 include 'COMMON.CHAIN'
2196 include 'COMMON.DERIV'
2197 include 'COMMON.GEO'
2198 include 'COMMON.LOCAL'
2199 include 'COMMON.INTERACT'
2200 include 'COMMON.IOUNITS'
2201 include 'COMMON.NAMES'
2202 double precision stochforcvec(MAXRES6)
2203 common /stochcalc/ stochforcvec
2204 logical lprn /.false./
2206 c write (iout,*) "dc_old"
2208 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2209 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2212 dc_work(j)=dc_old(j,0)
2213 d_t_work(j)=d_t_old(j,0)
2214 d_a_work(j)=d_a_old(j,0)
2219 dc_work(ind+j)=dc_old(j,i)
2220 d_t_work(ind+j)=d_t_old(j,i)
2221 d_a_work(ind+j)=d_a_old(j,i)
2226 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2228 dc_work(ind+j)=dc_old(j,i+nres)
2229 d_t_work(ind+j)=d_t_old(j,i+nres)
2230 d_a_work(ind+j)=d_a_old(j,i+nres)
2238 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2242 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2243 & vfric_mat(i,j),afric_mat(i,j),
2244 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2252 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2253 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2254 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2255 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2257 d_t_work_new(i)=ddt1+0.5d0*ddt2
2258 d_t_work(i)=ddt1+ddt2
2263 d_t(j,0)=d_t_work(j)
2268 dc(j,i)=dc_work(ind+j)
2269 d_t(j,i)=d_t_work(ind+j)
2274 if (itype(i).ne.10) then
2277 dc(j,inres)=dc_work(ind+j)
2278 d_t(j,inres)=d_t_work(ind+j)
2285 c--------------------------------------------------------------------------
2286 subroutine sd_verlet2
2287 c Calculating the adjusted velocities for accelerations
2288 implicit real*8 (a-h,o-z)
2289 include 'DIMENSIONS'
2290 include 'COMMON.CONTROL'
2291 include 'COMMON.VAR'
2294 include 'COMMON.LANGEVIN'
2296 include 'COMMON.LANGEVIN.lang0'
2298 include 'COMMON.CHAIN'
2299 include 'COMMON.DERIV'
2300 include 'COMMON.GEO'
2301 include 'COMMON.LOCAL'
2302 include 'COMMON.INTERACT'
2303 include 'COMMON.IOUNITS'
2304 include 'COMMON.NAMES'
2305 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2306 common /stochcalc/ stochforcvec
2308 c Compute the stochastic forces which contribute to velocity change
2310 call stochastic_force(stochforcvecV)
2317 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2318 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2319 & vrand_mat2(i,j)*stochforcvecV(j)
2321 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2325 d_t(j,0)=d_t_work(j)
2330 d_t(j,i)=d_t_work(ind+j)
2335 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2338 d_t(j,inres)=d_t_work(ind+j)
2345 c-----------------------------------------------------------
2346 subroutine sd_verlet_ciccotti_setup
2347 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2349 implicit real*8 (a-h,o-z)
2350 include 'DIMENSIONS'
2354 include 'COMMON.CONTROL'
2355 include 'COMMON.VAR'
2358 include 'COMMON.LANGEVIN'
2360 include 'COMMON.LANGEVIN.lang0'
2362 include 'COMMON.CHAIN'
2363 include 'COMMON.DERIV'
2364 include 'COMMON.GEO'
2365 include 'COMMON.LOCAL'
2366 include 'COMMON.INTERACT'
2367 include 'COMMON.IOUNITS'
2368 include 'COMMON.NAMES'
2369 include 'COMMON.TIME1'
2370 double precision emgdt(MAXRES6),
2371 & pterm,vterm,rho,rhoc,vsig,
2372 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2373 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2374 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2375 logical lprn /.false./
2376 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2377 double precision ktm
2384 c AL 8/17/04 Code adapted from tinker
2386 c Get the frictional and random terms for stochastic dynamics in the
2387 c eigenspace of mass-scaled UNRES friction matrix
2390 write (iout,*) "i",i," fricgam",fricgam(i)
2391 gdt = fricgam(i) * d_time
2393 c Stochastic dynamics reduces to simple MD for zero friction
2395 if (gdt .le. zero) then
2396 pfric_vec(i) = 1.0d0
2397 vfric_vec(i) = d_time
2398 afric_vec(i) = 0.5d0*d_time*d_time
2399 prand_vec(i) = afric_vec(i)
2400 vrand_vec2(i) = vfric_vec(i)
2402 c Analytical expressions when friction coefficient is large
2407 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2408 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2409 prand_vec(i) = afric_vec(i)
2410 vrand_vec2(i) = vfric_vec(i)
2412 c Compute the scaling factors of random terms for the nonzero friction case
2414 c ktm = 0.5d0*d_time/fricgam(i)
2415 c psig = dsqrt(ktm*pterm) / fricgam(i)
2416 c vsig = dsqrt(ktm*vterm)
2417 c prand_vec(i) = psig*afric_vec(i)
2418 c vrand_vec2(i) = vsig*vfric_vec(i)
2423 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2426 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2427 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2431 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2433 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2434 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2435 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2436 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2437 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2439 t_sdsetup=t_sdsetup+MPI_Wtime()
2441 t_sdsetup=t_sdsetup+tcpu()-tt0
2445 c-------------------------------------------------------------
2446 subroutine sd_verlet1_ciccotti
2447 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2448 implicit real*8 (a-h,o-z)
2449 include 'DIMENSIONS'
2453 include 'COMMON.CONTROL'
2454 include 'COMMON.VAR'
2457 include 'COMMON.LANGEVIN'
2459 include 'COMMON.LANGEVIN.lang0'
2461 include 'COMMON.CHAIN'
2462 include 'COMMON.DERIV'
2463 include 'COMMON.GEO'
2464 include 'COMMON.LOCAL'
2465 include 'COMMON.INTERACT'
2466 include 'COMMON.IOUNITS'
2467 include 'COMMON.NAMES'
2468 double precision stochforcvec(MAXRES6)
2469 common /stochcalc/ stochforcvec
2470 logical lprn /.false./
2472 c write (iout,*) "dc_old"
2474 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2475 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2478 dc_work(j)=dc_old(j,0)
2479 d_t_work(j)=d_t_old(j,0)
2480 d_a_work(j)=d_a_old(j,0)
2485 dc_work(ind+j)=dc_old(j,i)
2486 d_t_work(ind+j)=d_t_old(j,i)
2487 d_a_work(ind+j)=d_a_old(j,i)
2492 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2494 dc_work(ind+j)=dc_old(j,i+nres)
2495 d_t_work(ind+j)=d_t_old(j,i+nres)
2496 d_a_work(ind+j)=d_a_old(j,i+nres)
2505 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2509 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2510 & vfric_mat(i,j),afric_mat(i,j),
2511 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2519 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2520 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2521 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2522 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2524 d_t_work_new(i)=ddt1+0.5d0*ddt2
2525 d_t_work(i)=ddt1+ddt2
2530 d_t(j,0)=d_t_work(j)
2535 dc(j,i)=dc_work(ind+j)
2536 d_t(j,i)=d_t_work(ind+j)
2541 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2544 dc(j,inres)=dc_work(ind+j)
2545 d_t(j,inres)=d_t_work(ind+j)
2552 c--------------------------------------------------------------------------
2553 subroutine sd_verlet2_ciccotti
2554 c Calculating the adjusted velocities for accelerations
2555 implicit real*8 (a-h,o-z)
2556 include 'DIMENSIONS'
2557 include 'COMMON.CONTROL'
2558 include 'COMMON.VAR'
2561 include 'COMMON.LANGEVIN'
2563 include 'COMMON.LANGEVIN.lang0'
2565 include 'COMMON.CHAIN'
2566 include 'COMMON.DERIV'
2567 include 'COMMON.GEO'
2568 include 'COMMON.LOCAL'
2569 include 'COMMON.INTERACT'
2570 include 'COMMON.IOUNITS'
2571 include 'COMMON.NAMES'
2572 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2573 common /stochcalc/ stochforcvec
2575 c Compute the stochastic forces which contribute to velocity change
2577 call stochastic_force(stochforcvecV)
2584 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2585 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2586 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2588 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2592 d_t(j,0)=d_t_work(j)
2597 d_t(j,i)=d_t_work(ind+j)
2602 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2605 d_t(j,inres)=d_t_work(ind+j)