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))) 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(27)
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(27)
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),energia(0:n_ene)
1510 write (iout,*) "init_MD INDPDB",indpdb
1512 c write(iout,*) "d_time", d_time
1513 c Compute the standard deviations of stochastic forces for Langevin dynamics
1514 c if the friction coefficients do not depend on surface area
1515 if (lang.gt.0 .and. .not.surfarea) then
1517 stdforcp(i)=stdfp*dsqrt(gamp)
1520 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1521 & *dsqrt(gamsc(iabs(itype(i))))
1524 c Open the pdb file for snapshotshots
1527 if (ilen(tmpdir).gt.0)
1528 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1529 & liczba(:ilen(liczba))//".pdb")
1531 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1535 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1536 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1537 & liczba(:ilen(liczba))//".x")
1538 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1541 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1542 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1543 & liczba(:ilen(liczba))//".cx")
1544 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1550 if (ilen(tmpdir).gt.0)
1551 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1552 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1554 if (ilen(tmpdir).gt.0)
1555 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1556 cartname=prefix(:ilen(prefix))//"_MD.cx"
1560 write (qstr,'(256(1h ))')
1563 iq = qinfrag(i,iset)*10
1564 iw = wfrag(i,iset)/100
1566 if(me.eq.king.or..not.out1file)
1567 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1568 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1573 iq = qinpair(i,iset)*10
1574 iw = wpair(i,iset)/100
1576 if(me.eq.king.or..not.out1file)
1577 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1578 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1582 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1584 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1586 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1588 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1591 write (iout,*) "REST ",rest
1593 if (restart1file) then
1595 & inquire(file=mremd_rst_name,exist=file_exist)
1597 write (*,*) me," Before broadcast: file_exist",file_exist
1598 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1600 write (*,*) me," After broadcast: file_exist",file_exist
1601 c inquire(file=mremd_rst_name,exist=file_exist)
1603 if(me.eq.king.or..not.out1file)
1604 & write(iout,*) "Initial state read by master and distributed"
1606 if (ilen(tmpdir).gt.0)
1607 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1608 & //liczba(:ilen(liczba))//'.rst')
1609 inquire(file=rest2name,exist=file_exist)
1612 if(.not.restart1file) then
1613 if(me.eq.king.or..not.out1file)
1614 & write(iout,*) "Initial state will be read from file ",
1615 & rest2name(:ilen(rest2name))
1618 call rescale_weights(t_bath)
1621 if(me.eq.king.or..not.out1file)then
1622 if (restart1file) then
1623 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1626 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1629 write(iout,*) "Initial velocities randomly generated"
1636 c Generate initial velocities
1637 if(me.eq.king.or..not.out1file)
1638 & write(iout,*) "Initial velocities randomly generated"
1641 CtotTafm is the variable for AFM time which eclipsed during
1644 c rest2name = prefix(:ilen(prefix))//'.rst'
1645 if(me.eq.king.or..not.out1file)then
1646 write (iout,*) "Initial velocities"
1648 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1649 & (d_t(j,i+nres),j=1,3)
1651 c Zeroing the total angular momentum of the system
1652 write(iout,*) "Calling the zero-angular
1653 & momentum subroutine"
1656 c Getting the potential energy and forces and velocities and accelerations
1658 write (iout,*) "velocity of the center of the mass:"
1659 write (iout,*) (vcm(j),j=1,3)
1662 d_t(j,0)=d_t(j,0)-vcm(j)
1664 c Removing the velocity of the center of mass
1666 if(me.eq.king.or..not.out1file)then
1667 write (iout,*) "vcm right after adjustment:"
1668 write (iout,*) (vcm(j),j=1,3)
1671 write (iout,*) "init_MD before initial structure REST ",rest
1673 if (iranconf.ne.0) then
1674 c 8/22/17 AL Loop to produce a low-energy random conformation
1678 print *, 'Calling OVERLAP_SC'
1679 call overlap_sc(fail)
1683 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1684 print *,'SC_move',nft_sc,etot
1685 if(me.eq.king.or..not.out1file)
1686 & write(iout,*) 'SC_move',nft_sc,etot
1690 if (me.eq.king.or..not.out1file) write(iout,*)
1691 & 'Minimizing random structure: Calling MINIM_DC'
1692 call minim_dc(etot,iretcode,nfun)
1694 call geom_to_var(nvar,varia)
1695 if (me.eq.king.or..not.out1file) write (iout,*)
1696 & 'Minimizing random structure: Calling MINIMIZE.'
1697 call minimize(etot,varia,iretcode,nfun)
1698 call var_to_geom(nvar,varia)
1700 if (me.eq.king.or..not.out1file)
1701 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1702 if (isnan(etot) .or. etot.gt.1.0d4) then
1703 write (iout,*) "Energy too large",etot,
1704 & " trying another random conformation"
1707 call gen_rand_conf(itmp,*30)
1709 30 write (iout,*) 'Failed to generate random conformation',
1710 & ', itrial=',itrial
1711 write (*,*) 'Processor:',me,
1712 & ' Failed to generate random conformation',
1722 write (iout,'(a,i3,a)') 'Processor:',me,
1723 & ' error in generating random conformation.'
1724 write (*,'(a,i3,a)') 'Processor:',me,
1725 & ' error in generating random conformation.'
1728 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1737 write (iout,'(a,i3,a)') 'Processor:',me,
1738 & ' failed to generate a low-energy random conformation.'
1739 write (*,'(a,i3,a)') 'Processor:',me,
1740 & ' failed to generate a low-energy random conformation.'
1743 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1748 else if (preminim) then
1749 if (start_from_model) then
1750 i_model=iran_num(1,constr_homology)
1751 write (iout,*) 'starting from model ',i_model
1754 c(j,i)=chomo(j,i,i_model)
1757 call int_from_cart(.true.,.false.)
1758 call sc_loc_geom(.false.)
1761 dc(j,i)=c(j,i+1)-c(j,i)
1762 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1767 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1768 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1772 ! Remove SC overlaps if requested
1774 write (iout,*) 'Calling OVERLAP_SC'
1775 call overlap_sc(fail)
1777 ! Search for better SC rotamers if requested
1779 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1780 print *,'SC_move',nft_sc,etot
1781 if (me.eq.king.or..not.out1file)
1782 & write(iout,*) 'SC_move',nft_sc,etot
1784 call etotal(energia(0))
1785 C 8/22/17 AL Minimize initial structure
1787 if (me.eq.king.or..not.out1file) write(iout,*)
1788 & 'Minimizing initial PDB structure: Calling MINIM_DC'
1789 call minim_dc(etot,iretcode,nfun)
1791 call geom_to_var(nvar,varia)
1792 if(me.eq.king.or..not.out1file) write (iout,*)
1793 & 'Minimizing initial PDB structure: Calling MINIMIZE.'
1794 call minimize(etot,varia,iretcode,nfun)
1795 call var_to_geom(nvar,varia)
1797 if (me.eq.king.or..not.out1file)
1798 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1799 if(me.eq.king.or..not.out1file)
1800 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1803 call chainbuild_cart
1808 kinetic_T=2.0d0/(dimen3*Rb)*EK
1809 if(me.eq.king.or..not.out1file)then
1819 call etotal(potEcomp)
1820 if (large) call enerprint(potEcomp)
1823 t_etotal=t_etotal+MPI_Wtime()-tt0
1825 t_etotal=t_etotal+tcpu()-tt0
1829 c write (iout,*) "PotE-homology",potE-potEcomp(27)
1833 if (amax*d_time .gt. dvmax) then
1834 d_time=d_time*dvmax/amax
1835 if(me.eq.king.or..not.out1file) write (iout,*)
1836 & "Time step reduced to",d_time,
1837 & " because of too large initial acceleration."
1839 if(me.eq.king.or..not.out1file)then
1840 write(iout,*) "Potential energy and its components"
1841 call enerprint(potEcomp)
1842 c write(iout,*) (potEcomp(i),i=0,n_ene)
1844 potE=potEcomp(0)-potEcomp(27)
1845 c write (iout,*) "PotE-homology",potE
1848 if (ntwe.ne.0) call statout(itime)
1849 if(me.eq.king.or..not.out1file)
1850 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1851 & " Kinetic energy",EK," Potential energy",potE,
1852 & " Total energy",totE," Maximum acceleration ",
1855 write (iout,*) "Initial coordinates"
1857 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1858 & (c(j,i+nres),j=1,3)
1860 write (iout,*) "Initial dC"
1862 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1863 & (dc(j,i+nres),j=1,3)
1865 write (iout,*) "Initial velocities"
1866 write (iout,"(13x,' backbone ',23x,' side chain')")
1868 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1869 & (d_t(j,i+nres),j=1,3)
1871 write (iout,*) "Initial accelerations"
1873 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1874 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1875 & (d_a(j,i+nres),j=1,3)
1881 d_t_old(j,i)=d_t(j,i)
1882 d_a_old(j,i)=d_a(j,i)
1884 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1893 call etotal_short(energia_short)
1894 if (large) call enerprint(potEcomp)
1897 t_eshort=t_eshort+MPI_Wtime()-tt0
1899 t_eshort=t_eshort+tcpu()-tt0
1904 if(.not.out1file .and. large) then
1905 write (iout,*) "energia_long",energia_long(0),
1906 & " energia_short",energia_short(0),
1907 & " total",energia_long(0)+energia_short(0)
1908 write (iout,*) "Initial fast-force accelerations"
1910 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1911 & (d_a(j,i+nres),j=1,3)
1914 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1917 d_a_short(j,i)=d_a(j,i)
1926 call etotal_long(energia_long)
1927 if (large) call enerprint(potEcomp)
1930 t_elong=t_elong+MPI_Wtime()-tt0
1932 t_elong=t_elong+tcpu()-tt0
1937 if(.not.out1file .and. large) then
1938 write (iout,*) "energia_long",energia_long(0)
1939 write (iout,*) "Initial slow-force accelerations"
1941 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1942 & (d_a(j,i+nres),j=1,3)
1946 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1948 t_enegrad=t_enegrad+tcpu()-tt0
1953 c-----------------------------------------------------------
1954 subroutine random_vel
1955 implicit real*8 (a-h,o-z)
1956 include 'DIMENSIONS'
1957 include 'COMMON.CONTROL'
1958 include 'COMMON.VAR'
1961 include 'COMMON.LANGEVIN'
1963 include 'COMMON.LANGEVIN.lang0'
1965 include 'COMMON.CHAIN'
1966 include 'COMMON.DERIV'
1967 include 'COMMON.GEO'
1968 include 'COMMON.LOCAL'
1969 include 'COMMON.INTERACT'
1970 include 'COMMON.IOUNITS'
1971 include 'COMMON.NAMES'
1972 include 'COMMON.TIME1'
1973 double precision xv,sigv,lowb,highb,vec_afm(3)
1974 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1975 c First generate velocities in the eigenspace of the G matrix
1976 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1983 sigv=dsqrt((Rb*t_bath)/geigen(i))
1986 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1988 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1989 c & " d_t_work_new",d_t_work_new(ii)
1992 C if (SELFGUIDE.gt.0) then
1995 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
1996 C distance=distance+vec_afm(j)**2
1998 C distance=dsqrt(distance)
2000 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
2001 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
2002 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
2003 C & d_t_work_new(j+(afmend-1)*3)
2014 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2017 c write (iout,*) "Ek from eigenvectors",Ek1
2019 c Transform velocities to UNRES coordinate space
2025 d_t_work(ind)=d_t_work(ind)
2026 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2028 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2032 c Transfer to the d_t vector
2034 d_t(j,0)=d_t_work(j)
2040 d_t(j,i)=d_t_work(ind)
2044 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2047 d_t(j,i+nres)=d_t_work(ind)
2052 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2053 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2058 c-----------------------------------------------------------
2059 subroutine sd_verlet_p_setup
2060 c Sets up the parameters of stochastic Verlet algorithm
2061 implicit real*8 (a-h,o-z)
2062 include 'DIMENSIONS'
2066 include 'COMMON.CONTROL'
2067 include 'COMMON.VAR'
2070 include 'COMMON.LANGEVIN'
2072 include 'COMMON.LANGEVIN.lang0'
2074 include 'COMMON.CHAIN'
2075 include 'COMMON.DERIV'
2076 include 'COMMON.GEO'
2077 include 'COMMON.LOCAL'
2078 include 'COMMON.INTERACT'
2079 include 'COMMON.IOUNITS'
2080 include 'COMMON.NAMES'
2081 include 'COMMON.TIME1'
2082 double precision emgdt(MAXRES6),
2083 & pterm,vterm,rho,rhoc,vsig,
2084 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2085 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2086 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2087 logical lprn /.false./
2088 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2089 double precision ktm
2096 c AL 8/17/04 Code adapted from tinker
2098 c Get the frictional and random terms for stochastic dynamics in the
2099 c eigenspace of mass-scaled UNRES friction matrix
2102 gdt = fricgam(i) * d_time
2104 c Stochastic dynamics reduces to simple MD for zero friction
2106 if (gdt .le. zero) then
2107 pfric_vec(i) = 1.0d0
2108 vfric_vec(i) = d_time
2109 afric_vec(i) = 0.5d0 * d_time * d_time
2110 prand_vec(i) = 0.0d0
2111 vrand_vec1(i) = 0.0d0
2112 vrand_vec2(i) = 0.0d0
2114 c Analytical expressions when friction coefficient is large
2117 if (gdt .ge. gdt_radius) then
2120 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2121 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2122 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2123 vterm = 1.0d0 - egdt**2
2124 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2126 c Use series expansions when friction coefficient is small
2137 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2138 & - gdt5/120.0d0 + gdt6/720.0d0
2139 & - gdt7/5040.0d0 + gdt8/40320.0d0
2140 & - gdt9/362880.0d0) / fricgam(i)**2
2141 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2142 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2143 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2144 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2145 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2146 & + 127.0d0*gdt9/90720.0d0
2147 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2148 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2149 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2150 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2151 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2152 & - 17.0d0*gdt2/1280.0d0
2153 & + 17.0d0*gdt3/6144.0d0
2154 & + 40967.0d0*gdt4/34406400.0d0
2155 & - 57203.0d0*gdt5/275251200.0d0
2156 & - 1429487.0d0*gdt6/13212057600.0d0)
2159 c Compute the scaling factors of random terms for the nonzero friction case
2161 ktm = 0.5d0*d_time/fricgam(i)
2162 psig = dsqrt(ktm*pterm) / fricgam(i)
2163 vsig = dsqrt(ktm*vterm)
2164 rhoc = dsqrt(1.0d0 - rho*rho)
2166 vrand_vec1(i) = vsig * rho
2167 vrand_vec2(i) = vsig * rhoc
2172 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2175 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2176 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2180 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2183 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2184 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2185 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2186 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2187 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2188 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2191 t_sdsetup=t_sdsetup+MPI_Wtime()
2193 t_sdsetup=t_sdsetup+tcpu()-tt0
2197 c-------------------------------------------------------------
2198 subroutine eigtransf1(n,ndim,ab,d,c)
2201 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2207 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2213 c-------------------------------------------------------------
2214 subroutine eigtransf(n,ndim,a,b,d,c)
2217 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2223 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2229 c-------------------------------------------------------------
2230 subroutine sd_verlet1
2231 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2232 implicit real*8 (a-h,o-z)
2233 include 'DIMENSIONS'
2234 include 'COMMON.CONTROL'
2235 include 'COMMON.VAR'
2238 include 'COMMON.LANGEVIN'
2240 include 'COMMON.LANGEVIN.lang0'
2242 include 'COMMON.CHAIN'
2243 include 'COMMON.DERIV'
2244 include 'COMMON.GEO'
2245 include 'COMMON.LOCAL'
2246 include 'COMMON.INTERACT'
2247 include 'COMMON.IOUNITS'
2248 include 'COMMON.NAMES'
2249 double precision stochforcvec(MAXRES6)
2250 common /stochcalc/ stochforcvec
2251 logical lprn /.false./
2253 c write (iout,*) "dc_old"
2255 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2256 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2259 dc_work(j)=dc_old(j,0)
2260 d_t_work(j)=d_t_old(j,0)
2261 d_a_work(j)=d_a_old(j,0)
2266 dc_work(ind+j)=dc_old(j,i)
2267 d_t_work(ind+j)=d_t_old(j,i)
2268 d_a_work(ind+j)=d_a_old(j,i)
2273 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2275 dc_work(ind+j)=dc_old(j,i+nres)
2276 d_t_work(ind+j)=d_t_old(j,i+nres)
2277 d_a_work(ind+j)=d_a_old(j,i+nres)
2285 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2289 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2290 & vfric_mat(i,j),afric_mat(i,j),
2291 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2299 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2300 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2301 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2302 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2304 d_t_work_new(i)=ddt1+0.5d0*ddt2
2305 d_t_work(i)=ddt1+ddt2
2310 d_t(j,0)=d_t_work(j)
2315 dc(j,i)=dc_work(ind+j)
2316 d_t(j,i)=d_t_work(ind+j)
2321 if (itype(i).ne.10) then
2324 dc(j,inres)=dc_work(ind+j)
2325 d_t(j,inres)=d_t_work(ind+j)
2332 c--------------------------------------------------------------------------
2333 subroutine sd_verlet2
2334 c Calculating the adjusted velocities for accelerations
2335 implicit real*8 (a-h,o-z)
2336 include 'DIMENSIONS'
2337 include 'COMMON.CONTROL'
2338 include 'COMMON.VAR'
2341 include 'COMMON.LANGEVIN'
2343 include 'COMMON.LANGEVIN.lang0'
2345 include 'COMMON.CHAIN'
2346 include 'COMMON.DERIV'
2347 include 'COMMON.GEO'
2348 include 'COMMON.LOCAL'
2349 include 'COMMON.INTERACT'
2350 include 'COMMON.IOUNITS'
2351 include 'COMMON.NAMES'
2352 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2353 common /stochcalc/ stochforcvec
2355 c Compute the stochastic forces which contribute to velocity change
2357 call stochastic_force(stochforcvecV)
2364 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2365 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2366 & vrand_mat2(i,j)*stochforcvecV(j)
2368 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2372 d_t(j,0)=d_t_work(j)
2377 d_t(j,i)=d_t_work(ind+j)
2382 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2385 d_t(j,inres)=d_t_work(ind+j)
2392 c-----------------------------------------------------------
2393 subroutine sd_verlet_ciccotti_setup
2394 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2396 implicit real*8 (a-h,o-z)
2397 include 'DIMENSIONS'
2401 include 'COMMON.CONTROL'
2402 include 'COMMON.VAR'
2405 include 'COMMON.LANGEVIN'
2407 include 'COMMON.LANGEVIN.lang0'
2409 include 'COMMON.CHAIN'
2410 include 'COMMON.DERIV'
2411 include 'COMMON.GEO'
2412 include 'COMMON.LOCAL'
2413 include 'COMMON.INTERACT'
2414 include 'COMMON.IOUNITS'
2415 include 'COMMON.NAMES'
2416 include 'COMMON.TIME1'
2417 double precision emgdt(MAXRES6),
2418 & pterm,vterm,rho,rhoc,vsig,
2419 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2420 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2421 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2422 logical lprn /.false./
2423 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2424 double precision ktm
2431 c AL 8/17/04 Code adapted from tinker
2433 c Get the frictional and random terms for stochastic dynamics in the
2434 c eigenspace of mass-scaled UNRES friction matrix
2437 write (iout,*) "i",i," fricgam",fricgam(i)
2438 gdt = fricgam(i) * d_time
2440 c Stochastic dynamics reduces to simple MD for zero friction
2442 if (gdt .le. zero) then
2443 pfric_vec(i) = 1.0d0
2444 vfric_vec(i) = d_time
2445 afric_vec(i) = 0.5d0*d_time*d_time
2446 prand_vec(i) = afric_vec(i)
2447 vrand_vec2(i) = vfric_vec(i)
2449 c Analytical expressions when friction coefficient is large
2454 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2455 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2456 prand_vec(i) = afric_vec(i)
2457 vrand_vec2(i) = vfric_vec(i)
2459 c Compute the scaling factors of random terms for the nonzero friction case
2461 c ktm = 0.5d0*d_time/fricgam(i)
2462 c psig = dsqrt(ktm*pterm) / fricgam(i)
2463 c vsig = dsqrt(ktm*vterm)
2464 c prand_vec(i) = psig*afric_vec(i)
2465 c vrand_vec2(i) = vsig*vfric_vec(i)
2470 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2473 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2474 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2478 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2480 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2481 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2482 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2483 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2484 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2486 t_sdsetup=t_sdsetup+MPI_Wtime()
2488 t_sdsetup=t_sdsetup+tcpu()-tt0
2492 c-------------------------------------------------------------
2493 subroutine sd_verlet1_ciccotti
2494 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2495 implicit real*8 (a-h,o-z)
2496 include 'DIMENSIONS'
2500 include 'COMMON.CONTROL'
2501 include 'COMMON.VAR'
2504 include 'COMMON.LANGEVIN'
2506 include 'COMMON.LANGEVIN.lang0'
2508 include 'COMMON.CHAIN'
2509 include 'COMMON.DERIV'
2510 include 'COMMON.GEO'
2511 include 'COMMON.LOCAL'
2512 include 'COMMON.INTERACT'
2513 include 'COMMON.IOUNITS'
2514 include 'COMMON.NAMES'
2515 double precision stochforcvec(MAXRES6)
2516 common /stochcalc/ stochforcvec
2517 logical lprn /.false./
2519 c write (iout,*) "dc_old"
2521 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2522 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2525 dc_work(j)=dc_old(j,0)
2526 d_t_work(j)=d_t_old(j,0)
2527 d_a_work(j)=d_a_old(j,0)
2532 dc_work(ind+j)=dc_old(j,i)
2533 d_t_work(ind+j)=d_t_old(j,i)
2534 d_a_work(ind+j)=d_a_old(j,i)
2539 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2541 dc_work(ind+j)=dc_old(j,i+nres)
2542 d_t_work(ind+j)=d_t_old(j,i+nres)
2543 d_a_work(ind+j)=d_a_old(j,i+nres)
2552 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2556 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2557 & vfric_mat(i,j),afric_mat(i,j),
2558 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2566 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2567 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2568 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2569 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2571 d_t_work_new(i)=ddt1+0.5d0*ddt2
2572 d_t_work(i)=ddt1+ddt2
2577 d_t(j,0)=d_t_work(j)
2582 dc(j,i)=dc_work(ind+j)
2583 d_t(j,i)=d_t_work(ind+j)
2588 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2591 dc(j,inres)=dc_work(ind+j)
2592 d_t(j,inres)=d_t_work(ind+j)
2599 c--------------------------------------------------------------------------
2600 subroutine sd_verlet2_ciccotti
2601 c Calculating the adjusted velocities for accelerations
2602 implicit real*8 (a-h,o-z)
2603 include 'DIMENSIONS'
2604 include 'COMMON.CONTROL'
2605 include 'COMMON.VAR'
2608 include 'COMMON.LANGEVIN'
2610 include 'COMMON.LANGEVIN.lang0'
2612 include 'COMMON.CHAIN'
2613 include 'COMMON.DERIV'
2614 include 'COMMON.GEO'
2615 include 'COMMON.LOCAL'
2616 include 'COMMON.INTERACT'
2617 include 'COMMON.IOUNITS'
2618 include 'COMMON.NAMES'
2619 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2620 common /stochcalc/ stochforcvec
2622 c Compute the stochastic forces which contribute to velocity change
2624 call stochastic_force(stochforcvecV)
2631 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2632 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2633 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2635 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2639 d_t(j,0)=d_t_work(j)
2644 d_t(j,i)=d_t_work(ind+j)
2649 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2652 d_t(j,inres)=d_t_work(ind+j)