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 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"
1655 c write (iout,*) "After inertia"
1657 c Getting the potential energy and forces and velocities and accelerations
1658 if(me.eq.king.or..not.out1file)then
1659 write(iout,*) "Calling the vcm_vel"
1663 write (iout,*) "velocity of the center of the mass:"
1664 write (iout,*) (vcm(j),j=1,3)
1667 d_t(j,0)=d_t(j,0)-vcm(j)
1669 c Removing the velocity of the center of mass
1670 if(me.eq.king.or..not.out1file)then
1671 write(iout,*) "Calling the vcm_vel"
1675 if(me.eq.king.or..not.out1file)then
1676 write (iout,*) "vcm right after adjustment:"
1677 write (iout,*) (vcm(j),j=1,3)
1681 if(iranconf.ne.0) then
1682 c 8/22/17 AL Loop to produce a low-energy random conformation
1684 write (iout,*) "iranmin",iranmin
1687 print *, 'Calling OVERLAP_SC'
1688 call overlap_sc(fail)
1692 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1693 print *,'SC_move',nft_sc,etot
1694 if(me.eq.king.or..not.out1file)
1695 & write(iout,*) 'SC_move',nft_sc,etot
1699 print *, 'Calling MINIM_DC'
1700 call minim_dc(etot,iretcode,nfun)
1702 call geom_to_var(nvar,varia)
1703 print *,'Calling MINIMIZE.'
1704 call minimize(etot,varia,iretcode,nfun)
1705 call var_to_geom(nvar,varia)
1707 if(me.eq.king.or..not.out1file)
1708 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1709 if (isnan(etot) .or. etot.gt.1.0d4) then
1710 write (iout,*) "Energy too large",etot,
1711 & " trying another random conformation"
1715 call gen_rand_conf(itmp,nrestmp,*30)
1717 30 write (iout,*) 'Failed to generate random conformation',
1718 & ', itrial=',itrial
1719 write (*,*) 'Processor:',me,
1720 & ' Failed to generate random conformation',
1730 write (iout,'(a,i3,a)') 'Processor:',me,
1731 & ' error in generating random conformation.'
1732 write (*,'(a,i3,a)') 'Processor:',me,
1733 & ' error in generating random conformation.'
1736 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1745 write (iout,'(a,i3,a)') 'Processor:',me,
1746 & ' failed to generate a low-energy random conformation.'
1747 write (*,'(a,i3,a)') 'Processor:',me,
1748 & ' failed to generate a low-energy random conformation.'
1751 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1756 else if (indpdb.gt.0) then
1757 C 8/22/17 AL Minimize initial PDB structure
1759 print *, 'Calling MINIM_DC'
1760 call minim_dc(etot,iretcode,nfun)
1762 call geom_to_var(nvar,varia)
1763 print *,'Calling MINIMIZE.'
1764 call minimize(etot,varia,iretcode,nfun)
1765 call var_to_geom(nvar,varia)
1767 if(me.eq.king.or..not.out1file)
1768 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1771 call chainbuild_cart
1776 kinetic_T=2.0d0/(dimen3*Rb)*EK
1777 if(me.eq.king.or..not.out1file)then
1787 call etotal(potEcomp)
1788 if (large) call enerprint(potEcomp)
1791 t_etotal=t_etotal+MPI_Wtime()-tt0
1793 t_etotal=t_etotal+tcpu()-tt0
1800 if (amax*d_time .gt. dvmax) then
1801 d_time=d_time*dvmax/amax
1802 if(me.eq.king.or..not.out1file) write (iout,*)
1803 & "Time step reduced to",d_time,
1804 & " because of too large initial acceleration."
1806 if(me.eq.king.or..not.out1file)then
1807 write(iout,*) "Potential energy and its components"
1808 call enerprint(potEcomp)
1809 write(iout,*) (potEcomp(i),i=0,n_ene)
1811 potE=potEcomp(0)-potEcomp(20)
1814 if (ntwe.ne.0) call statout(itime)
1815 if(me.eq.king.or..not.out1file)
1816 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1817 & " Kinetic energy",EK," Potential energy",potE,
1818 & " Total energy",totE," Maximum acceleration ",
1821 write (iout,*) "Initial coordinates"
1823 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1824 & (c(j,i+nres),j=1,3)
1826 write (iout,*) "Initial dC"
1828 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1829 & (dc(j,i+nres),j=1,3)
1831 write (iout,*) "Initial velocities"
1832 write (iout,"(13x,' backbone ',23x,' side chain')")
1834 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1835 & (d_t(j,i+nres),j=1,3)
1837 write (iout,*) "Initial accelerations"
1839 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1840 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1841 & (d_a(j,i+nres),j=1,3)
1847 d_t_old(j,i)=d_t(j,i)
1848 d_a_old(j,i)=d_a(j,i)
1850 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1859 call etotal_short(energia_short)
1860 if (large) call enerprint(potEcomp)
1863 t_eshort=t_eshort+MPI_Wtime()-tt0
1865 t_eshort=t_eshort+tcpu()-tt0
1870 if(.not.out1file .and. large) then
1871 write (iout,*) "energia_long",energia_long(0),
1872 & " energia_short",energia_short(0),
1873 & " total",energia_long(0)+energia_short(0)
1874 write (iout,*) "Initial fast-force accelerations"
1876 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1877 & (d_a(j,i+nres),j=1,3)
1880 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1883 d_a_short(j,i)=d_a(j,i)
1892 call etotal_long(energia_long)
1893 if (large) call enerprint(potEcomp)
1896 t_elong=t_elong+MPI_Wtime()-tt0
1898 t_elong=t_elong+tcpu()-tt0
1903 if(.not.out1file .and. large) then
1904 write (iout,*) "energia_long",energia_long(0)
1905 write (iout,*) "Initial slow-force accelerations"
1907 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1908 & (d_a(j,i+nres),j=1,3)
1912 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1914 t_enegrad=t_enegrad+tcpu()-tt0
1919 c-----------------------------------------------------------
1920 subroutine random_vel
1921 implicit real*8 (a-h,o-z)
1922 include 'DIMENSIONS'
1923 include 'COMMON.CONTROL'
1924 include 'COMMON.VAR'
1927 include 'COMMON.LANGEVIN'
1929 include 'COMMON.LANGEVIN.lang0'
1931 include 'COMMON.CHAIN'
1932 include 'COMMON.DERIV'
1933 include 'COMMON.GEO'
1934 include 'COMMON.LOCAL'
1935 include 'COMMON.INTERACT'
1936 include 'COMMON.IOUNITS'
1937 include 'COMMON.NAMES'
1938 include 'COMMON.TIME1'
1939 double precision xv,sigv,lowb,highb,vec_afm(3)
1940 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1941 c First generate velocities in the eigenspace of the G matrix
1942 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1949 sigv=dsqrt((Rb*t_bath)/geigen(i))
1952 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1954 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1955 c & " d_t_work_new",d_t_work_new(ii)
1958 C if (SELFGUIDE.gt.0) then
1961 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
1962 C distance=distance+vec_afm(j)**2
1964 C distance=dsqrt(distance)
1966 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
1967 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
1968 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
1969 C & d_t_work_new(j+(afmend-1)*3)
1980 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1983 c write (iout,*) "Ek from eigenvectors",Ek1
1985 c Transform velocities to UNRES coordinate space
1991 d_t_work(ind)=d_t_work(ind)
1992 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1994 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1998 c Transfer to the d_t vector
2000 d_t(j,0)=d_t_work(j)
2006 d_t(j,i)=d_t_work(ind)
2010 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2013 d_t(j,i+nres)=d_t_work(ind)
2018 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2019 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2024 c-----------------------------------------------------------
2025 subroutine sd_verlet_p_setup
2026 c Sets up the parameters of stochastic Verlet algorithm
2027 implicit real*8 (a-h,o-z)
2028 include 'DIMENSIONS'
2032 include 'COMMON.CONTROL'
2033 include 'COMMON.VAR'
2036 include 'COMMON.LANGEVIN'
2038 include 'COMMON.LANGEVIN.lang0'
2040 include 'COMMON.CHAIN'
2041 include 'COMMON.DERIV'
2042 include 'COMMON.GEO'
2043 include 'COMMON.LOCAL'
2044 include 'COMMON.INTERACT'
2045 include 'COMMON.IOUNITS'
2046 include 'COMMON.NAMES'
2047 include 'COMMON.TIME1'
2048 double precision emgdt(MAXRES6),
2049 & pterm,vterm,rho,rhoc,vsig,
2050 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2051 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2052 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2053 logical lprn /.false./
2054 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2055 double precision ktm
2062 c AL 8/17/04 Code adapted from tinker
2064 c Get the frictional and random terms for stochastic dynamics in the
2065 c eigenspace of mass-scaled UNRES friction matrix
2068 gdt = fricgam(i) * d_time
2070 c Stochastic dynamics reduces to simple MD for zero friction
2072 if (gdt .le. zero) then
2073 pfric_vec(i) = 1.0d0
2074 vfric_vec(i) = d_time
2075 afric_vec(i) = 0.5d0 * d_time * d_time
2076 prand_vec(i) = 0.0d0
2077 vrand_vec1(i) = 0.0d0
2078 vrand_vec2(i) = 0.0d0
2080 c Analytical expressions when friction coefficient is large
2083 if (gdt .ge. gdt_radius) then
2086 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2087 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2088 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2089 vterm = 1.0d0 - egdt**2
2090 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2092 c Use series expansions when friction coefficient is small
2103 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2104 & - gdt5/120.0d0 + gdt6/720.0d0
2105 & - gdt7/5040.0d0 + gdt8/40320.0d0
2106 & - gdt9/362880.0d0) / fricgam(i)**2
2107 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2108 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2109 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2110 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2111 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2112 & + 127.0d0*gdt9/90720.0d0
2113 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2114 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2115 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2116 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2117 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2118 & - 17.0d0*gdt2/1280.0d0
2119 & + 17.0d0*gdt3/6144.0d0
2120 & + 40967.0d0*gdt4/34406400.0d0
2121 & - 57203.0d0*gdt5/275251200.0d0
2122 & - 1429487.0d0*gdt6/13212057600.0d0)
2125 c Compute the scaling factors of random terms for the nonzero friction case
2127 ktm = 0.5d0*d_time/fricgam(i)
2128 psig = dsqrt(ktm*pterm) / fricgam(i)
2129 vsig = dsqrt(ktm*vterm)
2130 rhoc = dsqrt(1.0d0 - rho*rho)
2132 vrand_vec1(i) = vsig * rho
2133 vrand_vec2(i) = vsig * rhoc
2138 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2141 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2142 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2146 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2149 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2150 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2151 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2152 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2153 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2154 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2157 t_sdsetup=t_sdsetup+MPI_Wtime()
2159 t_sdsetup=t_sdsetup+tcpu()-tt0
2163 c-------------------------------------------------------------
2164 subroutine eigtransf1(n,ndim,ab,d,c)
2167 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2173 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2179 c-------------------------------------------------------------
2180 subroutine eigtransf(n,ndim,a,b,d,c)
2183 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2189 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2195 c-------------------------------------------------------------
2196 subroutine sd_verlet1
2197 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2198 implicit real*8 (a-h,o-z)
2199 include 'DIMENSIONS'
2200 include 'COMMON.CONTROL'
2201 include 'COMMON.VAR'
2204 include 'COMMON.LANGEVIN'
2206 include 'COMMON.LANGEVIN.lang0'
2208 include 'COMMON.CHAIN'
2209 include 'COMMON.DERIV'
2210 include 'COMMON.GEO'
2211 include 'COMMON.LOCAL'
2212 include 'COMMON.INTERACT'
2213 include 'COMMON.IOUNITS'
2214 include 'COMMON.NAMES'
2215 double precision stochforcvec(MAXRES6)
2216 common /stochcalc/ stochforcvec
2217 logical lprn /.false./
2219 c write (iout,*) "dc_old"
2221 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2222 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2225 dc_work(j)=dc_old(j,0)
2226 d_t_work(j)=d_t_old(j,0)
2227 d_a_work(j)=d_a_old(j,0)
2232 dc_work(ind+j)=dc_old(j,i)
2233 d_t_work(ind+j)=d_t_old(j,i)
2234 d_a_work(ind+j)=d_a_old(j,i)
2239 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2241 dc_work(ind+j)=dc_old(j,i+nres)
2242 d_t_work(ind+j)=d_t_old(j,i+nres)
2243 d_a_work(ind+j)=d_a_old(j,i+nres)
2251 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2255 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2256 & vfric_mat(i,j),afric_mat(i,j),
2257 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2265 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2266 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2267 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2268 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2270 d_t_work_new(i)=ddt1+0.5d0*ddt2
2271 d_t_work(i)=ddt1+ddt2
2276 d_t(j,0)=d_t_work(j)
2281 dc(j,i)=dc_work(ind+j)
2282 d_t(j,i)=d_t_work(ind+j)
2287 if (itype(i).ne.10) then
2290 dc(j,inres)=dc_work(ind+j)
2291 d_t(j,inres)=d_t_work(ind+j)
2298 c--------------------------------------------------------------------------
2299 subroutine sd_verlet2
2300 c Calculating the adjusted velocities for accelerations
2301 implicit real*8 (a-h,o-z)
2302 include 'DIMENSIONS'
2303 include 'COMMON.CONTROL'
2304 include 'COMMON.VAR'
2307 include 'COMMON.LANGEVIN'
2309 include 'COMMON.LANGEVIN.lang0'
2311 include 'COMMON.CHAIN'
2312 include 'COMMON.DERIV'
2313 include 'COMMON.GEO'
2314 include 'COMMON.LOCAL'
2315 include 'COMMON.INTERACT'
2316 include 'COMMON.IOUNITS'
2317 include 'COMMON.NAMES'
2318 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2319 common /stochcalc/ stochforcvec
2321 c Compute the stochastic forces which contribute to velocity change
2323 call stochastic_force(stochforcvecV)
2330 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2331 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2332 & vrand_mat2(i,j)*stochforcvecV(j)
2334 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2338 d_t(j,0)=d_t_work(j)
2343 d_t(j,i)=d_t_work(ind+j)
2348 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2351 d_t(j,inres)=d_t_work(ind+j)
2358 c-----------------------------------------------------------
2359 subroutine sd_verlet_ciccotti_setup
2360 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2362 implicit real*8 (a-h,o-z)
2363 include 'DIMENSIONS'
2367 include 'COMMON.CONTROL'
2368 include 'COMMON.VAR'
2371 include 'COMMON.LANGEVIN'
2373 include 'COMMON.LANGEVIN.lang0'
2375 include 'COMMON.CHAIN'
2376 include 'COMMON.DERIV'
2377 include 'COMMON.GEO'
2378 include 'COMMON.LOCAL'
2379 include 'COMMON.INTERACT'
2380 include 'COMMON.IOUNITS'
2381 include 'COMMON.NAMES'
2382 include 'COMMON.TIME1'
2383 double precision emgdt(MAXRES6),
2384 & pterm,vterm,rho,rhoc,vsig,
2385 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2386 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2387 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2388 logical lprn /.false./
2389 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2390 double precision ktm
2397 c AL 8/17/04 Code adapted from tinker
2399 c Get the frictional and random terms for stochastic dynamics in the
2400 c eigenspace of mass-scaled UNRES friction matrix
2403 write (iout,*) "i",i," fricgam",fricgam(i)
2404 gdt = fricgam(i) * d_time
2406 c Stochastic dynamics reduces to simple MD for zero friction
2408 if (gdt .le. zero) then
2409 pfric_vec(i) = 1.0d0
2410 vfric_vec(i) = d_time
2411 afric_vec(i) = 0.5d0*d_time*d_time
2412 prand_vec(i) = afric_vec(i)
2413 vrand_vec2(i) = vfric_vec(i)
2415 c Analytical expressions when friction coefficient is large
2420 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2421 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2422 prand_vec(i) = afric_vec(i)
2423 vrand_vec2(i) = vfric_vec(i)
2425 c Compute the scaling factors of random terms for the nonzero friction case
2427 c ktm = 0.5d0*d_time/fricgam(i)
2428 c psig = dsqrt(ktm*pterm) / fricgam(i)
2429 c vsig = dsqrt(ktm*vterm)
2430 c prand_vec(i) = psig*afric_vec(i)
2431 c vrand_vec2(i) = vsig*vfric_vec(i)
2436 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2439 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2440 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2444 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2446 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2447 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2448 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2449 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2450 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2452 t_sdsetup=t_sdsetup+MPI_Wtime()
2454 t_sdsetup=t_sdsetup+tcpu()-tt0
2458 c-------------------------------------------------------------
2459 subroutine sd_verlet1_ciccotti
2460 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2461 implicit real*8 (a-h,o-z)
2462 include 'DIMENSIONS'
2466 include 'COMMON.CONTROL'
2467 include 'COMMON.VAR'
2470 include 'COMMON.LANGEVIN'
2472 include 'COMMON.LANGEVIN.lang0'
2474 include 'COMMON.CHAIN'
2475 include 'COMMON.DERIV'
2476 include 'COMMON.GEO'
2477 include 'COMMON.LOCAL'
2478 include 'COMMON.INTERACT'
2479 include 'COMMON.IOUNITS'
2480 include 'COMMON.NAMES'
2481 double precision stochforcvec(MAXRES6)
2482 common /stochcalc/ stochforcvec
2483 logical lprn /.false./
2485 c write (iout,*) "dc_old"
2487 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2488 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2491 dc_work(j)=dc_old(j,0)
2492 d_t_work(j)=d_t_old(j,0)
2493 d_a_work(j)=d_a_old(j,0)
2498 dc_work(ind+j)=dc_old(j,i)
2499 d_t_work(ind+j)=d_t_old(j,i)
2500 d_a_work(ind+j)=d_a_old(j,i)
2505 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2507 dc_work(ind+j)=dc_old(j,i+nres)
2508 d_t_work(ind+j)=d_t_old(j,i+nres)
2509 d_a_work(ind+j)=d_a_old(j,i+nres)
2518 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2522 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2523 & vfric_mat(i,j),afric_mat(i,j),
2524 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2532 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2533 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2534 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2535 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2537 d_t_work_new(i)=ddt1+0.5d0*ddt2
2538 d_t_work(i)=ddt1+ddt2
2543 d_t(j,0)=d_t_work(j)
2548 dc(j,i)=dc_work(ind+j)
2549 d_t(j,i)=d_t_work(ind+j)
2554 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2557 dc(j,inres)=dc_work(ind+j)
2558 d_t(j,inres)=d_t_work(ind+j)
2565 c--------------------------------------------------------------------------
2566 subroutine sd_verlet2_ciccotti
2567 c Calculating the adjusted velocities for accelerations
2568 implicit real*8 (a-h,o-z)
2569 include 'DIMENSIONS'
2570 include 'COMMON.CONTROL'
2571 include 'COMMON.VAR'
2574 include 'COMMON.LANGEVIN'
2576 include 'COMMON.LANGEVIN.lang0'
2578 include 'COMMON.CHAIN'
2579 include 'COMMON.DERIV'
2580 include 'COMMON.GEO'
2581 include 'COMMON.LOCAL'
2582 include 'COMMON.INTERACT'
2583 include 'COMMON.IOUNITS'
2584 include 'COMMON.NAMES'
2585 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2586 common /stochcalc/ stochforcvec
2588 c Compute the stochastic forces which contribute to velocity change
2590 call stochastic_force(stochforcvecV)
2597 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2598 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2599 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2601 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2605 d_t(j,0)=d_t_work(j)
2610 d_t(j,i)=d_t_work(ind+j)
2615 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2618 d_t(j,inres)=d_t_work(ind+j)