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)
39 integer nharp,iharp(4,maxres/3)
42 if (ilen(tmpdir).gt.0)
43 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
44 & //liczba(:ilen(liczba))//'.rst')
46 if (ilen(tmpdir).gt.0)
47 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
54 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
56 & "Cartesian coordinates of the initial structure"
57 write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
58 & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
60 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
61 & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
62 & (c(j,ires+nres),j=1,3)
65 & "Initial dC vectors of the chain"
66 write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
67 & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
69 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
70 & restyp(itype(ires)),ires,(dc(j,ires),j=1,3),
71 & (dc(j,ires+nres),j=1,3)
74 & "Initial dC_norm vectors of the chain"
75 write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
76 & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
78 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
79 & restyp(itype(ires)),ires,(dc_norm(j,ires),j=1,3),
80 & (dc_norm(j,ires+nres),j=1,3)
87 c Determine the inverse of the inertia matrix.
88 call setup_MD_matrices
92 t_MDsetup = MPI_Wtime()-tt0
94 t_MDsetup = tcpu()-tt0
97 c Entering the MD loop
103 if (lang.eq.2 .or. lang.eq.3) then
107 call sd_verlet_p_setup
109 call sd_verlet_ciccotti_setup
113 pfric0_mat(i,j,0)=pfric_mat(i,j)
114 afric0_mat(i,j,0)=afric_mat(i,j)
115 vfric0_mat(i,j,0)=vfric_mat(i,j)
116 prand0_mat(i,j,0)=prand_mat(i,j)
117 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
118 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
123 flag_stoch(i)=.false.
127 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
129 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
133 else if (lang.eq.1 .or. lang.eq.4) then
137 t_langsetup=MPI_Wtime()-tt0
140 t_langsetup=tcpu()-tt0
143 do itime=1,n_timestep
145 if (large.and. mod(itime,ntwe).eq.0)
146 & write (iout,*) "itime",itime
148 if (lang.gt.0 .and. surfarea .and.
149 & mod(itime,reset_fricmat).eq.0) then
150 if (lang.eq.2 .or. lang.eq.3) then
154 call sd_verlet_p_setup
156 call sd_verlet_ciccotti_setup
160 pfric0_mat(i,j,0)=pfric_mat(i,j)
161 afric0_mat(i,j,0)=afric_mat(i,j)
162 vfric0_mat(i,j,0)=vfric_mat(i,j)
163 prand0_mat(i,j,0)=prand_mat(i,j)
164 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
165 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
170 flag_stoch(i)=.false.
173 else if (lang.eq.1 .or. lang.eq.4) then
176 write (iout,'(a,i10)')
177 & "Friction matrix reset based on surface area, itime",itime
179 if (reset_vel .and. tbf .and. lang.eq.0
180 & .and. mod(itime,count_reset_vel).eq.0) then
182 write(iout,'(a,f20.2)')
183 & "Velocities reset to random values, time",totT
186 d_t_old(j,i)=d_t(j,i)
190 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
194 d_t(j,0)=d_t(j,0)-vcm(j)
197 kinetic_T=2.0d0/(dimen3*Rb)*EK
198 scalfac=dsqrt(T_bath/kinetic_T)
199 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
202 d_t_old(j,i)=scalfac*d_t(j,i)
208 c Time-reversible RESPA algorithm
209 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
210 call RESPA_step(itime)
212 c Variable time step algorithm.
213 call velverlet_step(itime)
217 call brown_step(itime)
219 print *,"Brown dynamics not here!"
221 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
227 if (mod(itime,ntwe).eq.0) then
229 C call enerprint(potEcomp)
230 C print *,itime,'AFM',Eafmforc,etot
244 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
247 v_work(ind)=d_t(j,i+nres)
252 write (66,'(80f10.5)')
253 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
257 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
259 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
261 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
264 if (mod(itime,ntwx).eq.0) then
265 write(iout,*) 'time=',itime
266 C call check_ecartint
268 write (tytul,'("time",f8.2)') totT
270 call hairpin(.true.,nharp,iharp)
271 call secondary2(.true.)
272 call pdbout(potE,tytul,ipdb)
277 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
278 open(irest2,file=rest2name,status='unknown')
279 write(irest2,*) totT,EK,potE,totE,t_bath
281 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
284 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
295 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
297 & 'MD calculations setup:',t_MDsetup,
298 & 'Energy & gradient evaluation:',t_enegrad,
299 & 'Stochastic MD setup:',t_langsetup,
300 & 'Stochastic MD step setup:',t_sdsetup,
302 write (iout,'(/28(1h=),a25,27(1h=))')
303 & ' End of MD calculation '
305 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
307 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
308 & " time_fricmatmult",time_fricmatmult," time_fsample ",
313 c-------------------------------------------------------------------------------
314 subroutine velverlet_step(itime)
315 c-------------------------------------------------------------------------------
316 c Perform a single velocity Verlet step; the time step can be rescaled if
317 c increments in accelerations exceed the threshold
318 c-------------------------------------------------------------------------------
319 implicit real*8 (a-h,o-z)
323 integer ierror,ierrcode
325 include 'COMMON.SETUP'
326 include 'COMMON.CONTROL'
330 include 'COMMON.LANGEVIN'
332 include 'COMMON.LANGEVIN.lang0'
334 include 'COMMON.CHAIN'
335 include 'COMMON.DERIV'
337 include 'COMMON.LOCAL'
338 include 'COMMON.INTERACT'
339 include 'COMMON.IOUNITS'
340 include 'COMMON.NAMES'
341 include 'COMMON.TIME1'
342 include 'COMMON.MUCA'
343 double precision vcm(3),incr(3)
344 double precision cm(3),L(3)
345 integer ilen,count,rstcount
348 integer maxcount_scale /20/
350 double precision stochforcvec(MAXRES6)
351 common /stochcalc/ stochforcvec
359 else if (lang.eq.2 .or. lang.eq.3) then
361 call stochastic_force(stochforcvec)
364 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
366 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
373 icount_scale=icount_scale+1
374 if (icount_scale.gt.maxcount_scale) then
376 & "ERROR: too many attempts at scaling down the time step. ",
377 & "amax=",amax,"epdrift=",epdrift,
378 & "damax=",damax,"edriftmax=",edriftmax,
382 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
386 c First step of the velocity Verlet algorithm
391 else if (lang.eq.3) then
393 call sd_verlet1_ciccotti
395 else if (lang.eq.1) then
400 c Build the chain from the newly calculated coordinates
402 if (rattle) call rattle1
404 if (large.and. mod(itime,ntwe).eq.0) then
405 write (iout,*) "Cartesian and internal coordinates: step 1"
410 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
411 & (dc(j,i+nres),j=1,3)
413 write (iout,*) "Accelerations"
415 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
416 & (d_a(j,i+nres),j=1,3)
418 write (iout,*) "Velocities, step 1"
420 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
421 & (d_t(j,i+nres),j=1,3)
430 c Calculate energy and forces
432 call etotal(potEcomp)
433 if (large.and. mod(itime,ntwe).eq.0)
434 & call enerprint(potEcomp)
437 t_etotal=t_etotal+MPI_Wtime()-tt0
439 t_etotal=t_etotal+tcpu()-tt0
442 potE=potEcomp(0)-potEcomp(20)
444 c Get the new accelerations
447 t_enegrad=t_enegrad+MPI_Wtime()-tt0
449 t_enegrad=t_enegrad+tcpu()-tt0
451 c Determine maximum acceleration and scale down the timestep if needed
453 amax=amax/(itime_scal+1)**2
454 call predict_edrift(epdrift)
455 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
456 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
458 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
460 itime_scal=itime_scal+ifac_time
461 c fac_time=dmin1(damax/amax,0.5d0)
462 fac_time=0.5d0**ifac_time
463 d_time=d_time*fac_time
464 if (lang.eq.2 .or. lang.eq.3) then
466 c write (iout,*) "Calling sd_verlet_setup: 1"
467 c Rescale the stochastic forces and recalculate or restore
468 c the matrices of tinker integrator
469 if (itime_scal.gt.maxflag_stoch) then
470 if (large) write (iout,'(a,i5,a)')
471 & "Calculate matrices for stochastic step;",
472 & " itime_scal ",itime_scal
474 call sd_verlet_p_setup
476 call sd_verlet_ciccotti_setup
478 write (iout,'(2a,i3,a,i3,1h.)')
479 & "Warning: cannot store matrices for stochastic",
480 & " integration because the index",itime_scal,
481 & " is greater than",maxflag_stoch
482 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
483 & " integration Langevin algorithm for better efficiency."
484 else if (flag_stoch(itime_scal)) then
485 if (large) write (iout,'(a,i5,a,l1)')
486 & "Restore matrices for stochastic step; itime_scal ",
487 & itime_scal," flag ",flag_stoch(itime_scal)
490 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
491 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
492 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
493 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
494 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
495 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
499 if (large) write (iout,'(2a,i5,a,l1)')
500 & "Calculate & store matrices for stochastic step;",
501 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
503 call sd_verlet_p_setup
505 call sd_verlet_ciccotti_setup
507 flag_stoch(ifac_time)=.true.
510 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
511 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
512 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
513 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
514 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
515 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
519 fac_time=1.0d0/dsqrt(fac_time)
521 stochforcvec(i)=fac_time*stochforcvec(i)
524 else if (lang.eq.1) then
525 c Rescale the accelerations due to stochastic forces
526 fac_time=1.0d0/dsqrt(fac_time)
528 d_as_work(i)=d_as_work(i)*fac_time
531 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
532 & "itime",itime," Timestep scaled down to ",
533 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
535 c Second step of the velocity Verlet algorithm
540 else if (lang.eq.3) then
542 call sd_verlet2_ciccotti
544 else if (lang.eq.1) then
549 if (rattle) call rattle2
552 C print *,totTafm,"TU?"
553 if (d_time.ne.d_time0) then
556 if (lang.eq.2 .or. lang.eq.3) then
557 if (large) write (iout,'(a)')
558 & "Restore original matrices for stochastic step"
559 c write (iout,*) "Calling sd_verlet_setup: 2"
560 c Restore the matrices of tinker integrator if the time step has been restored
563 pfric_mat(i,j)=pfric0_mat(i,j,0)
564 afric_mat(i,j)=afric0_mat(i,j,0)
565 vfric_mat(i,j)=vfric0_mat(i,j,0)
566 prand_mat(i,j)=prand0_mat(i,j,0)
567 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
568 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
577 c Calculate the kinetic and the total energy and the kinetic temperature
582 c write (iout,*) "step",itime," EK",EK," EK1",EK1
584 c Couple the system to Berendsen bath if needed
585 if (tbf .and. lang.eq.0) then
588 kinetic_T=2.0d0/(dimen3*Rb)*EK
589 c Backup the coordinates, velocities, and accelerations
593 d_t_old(j,i)=d_t(j,i)
594 d_a_old(j,i)=d_a(j,i)
598 if (mod(itime,ntwe).eq.0 .and. large) then
599 write (iout,*) "Velocities, step 2"
601 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
602 & (d_t(j,i+nres),j=1,3)
608 c-------------------------------------------------------------------------------
609 subroutine RESPA_step(itime)
610 c-------------------------------------------------------------------------------
611 c Perform a single RESPA step.
612 c-------------------------------------------------------------------------------
613 implicit real*8 (a-h,o-z)
617 integer IERROR,ERRCODE
619 include 'COMMON.SETUP'
620 include 'COMMON.CONTROL'
624 include 'COMMON.LANGEVIN'
626 include 'COMMON.LANGEVIN.lang0'
628 include 'COMMON.CHAIN'
629 include 'COMMON.DERIV'
631 include 'COMMON.LOCAL'
632 include 'COMMON.INTERACT'
633 include 'COMMON.IOUNITS'
634 include 'COMMON.NAMES'
635 include 'COMMON.TIME1'
636 double precision energia_short(0:n_ene),
637 & energia_long(0:n_ene)
638 double precision cm(3),L(3),vcm(3),incr(3)
639 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
640 & d_a_old0(3,0:maxres2)
641 logical PRINT_AMTS_MSG /.false./
642 integer ilen,count,rstcount
645 integer maxcount_scale /10/
647 double precision stochforcvec(MAXRES6)
648 common /stochcalc/ stochforcvec
651 common /cipiszcze/ itt
654 if (large.and. mod(itime,ntwe).eq.0) then
655 write (iout,*) "***************** RESPA itime",itime
656 write (iout,*) "Cartesian and internal coordinates: step 0"
658 call pdbout(0.0d0,"cipiszcze",iout)
660 write (iout,*) "Accelerations from long-range forces"
662 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
663 & (d_a(j,i+nres),j=1,3)
665 write (iout,*) "Velocities, step 0"
667 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
668 & (d_t(j,i+nres),j=1,3)
673 c Perform the initial RESPA step (increment velocities)
674 c write (iout,*) "*********************** RESPA ini"
677 if (mod(itime,ntwe).eq.0 .and. large) then
678 write (iout,*) "Velocities, end"
680 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
681 & (d_t(j,i+nres),j=1,3)
685 c Compute the short-range forces
691 C 7/2/2009 commented out
693 c call etotal_short(energia_short)
696 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
699 d_a(j,i)=d_a_short(j,i)
703 if (large.and. mod(itime,ntwe).eq.0) then
704 write (iout,*) "energia_short",energia_short(0)
705 write (iout,*) "Accelerations from short-range forces"
707 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
708 & (d_a(j,i+nres),j=1,3)
713 t_enegrad=t_enegrad+MPI_Wtime()-tt0
715 t_enegrad=t_enegrad+tcpu()-tt0
720 d_t_old(j,i)=d_t(j,i)
721 d_a_old(j,i)=d_a(j,i)
724 c 6/30/08 A-MTS: attempt at increasing the split number
727 dc_old0(j,i)=dc_old(j,i)
728 d_t_old0(j,i)=d_t_old(j,i)
729 d_a_old0(j,i)=d_a_old(j,i)
732 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
733 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
740 c write (iout,*) "itime",itime," ntime_split",ntime_split
741 c Split the time step
742 d_time=d_time0/ntime_split
743 c Perform the short-range RESPA steps (velocity Verlet increments of
744 c positions and velocities using short-range forces)
745 c write (iout,*) "*********************** RESPA split"
746 do itsplit=1,ntime_split
749 else if (lang.eq.2 .or. lang.eq.3) then
751 call stochastic_force(stochforcvec)
754 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
756 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
761 c First step of the velocity Verlet algorithm
766 else if (lang.eq.3) then
768 call sd_verlet1_ciccotti
770 else if (lang.eq.1) then
775 c Build the chain from the newly calculated coordinates
777 if (rattle) call rattle1
779 if (large.and. mod(itime,ntwe).eq.0) then
780 write (iout,*) "***** ITSPLIT",itsplit
781 write (iout,*) "Cartesian and internal coordinates: step 1"
782 call pdbout(0.0d0,"cipiszcze",iout)
785 write (iout,*) "Velocities, step 1"
787 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
788 & (d_t(j,i+nres),j=1,3)
797 c Calculate energy and forces
799 call etotal_short(energia_short)
800 if (large.and. mod(itime,ntwe).eq.0)
801 & call enerprint(energia_short)
804 t_eshort=t_eshort+MPI_Wtime()-tt0
806 t_eshort=t_eshort+tcpu()-tt0
810 c Get the new accelerations
812 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
815 d_a_short(j,i)=d_a(j,i)
819 if (large.and. mod(itime,ntwe).eq.0) then
820 write (iout,*)"energia_short",energia_short(0)
821 write (iout,*) "Accelerations from short-range forces"
823 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
824 & (d_a(j,i+nres),j=1,3)
829 c Determine maximum acceleration and scale down the timestep if needed
831 amax=amax/ntime_split**2
832 call predict_edrift(epdrift)
833 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
834 & write (iout,*) "amax",amax," damax",damax,
835 & " epdrift",epdrift," epdriftmax",epdriftmax
836 c Exit loop and try with increased split number if the change of
837 c acceleration is too big
838 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
839 if (ntime_split.lt.maxtime_split) then
841 ntime_split=ntime_split*2
844 dc_old(j,i)=dc_old0(j,i)
845 d_t_old(j,i)=d_t_old0(j,i)
846 d_a_old(j,i)=d_a_old0(j,i)
849 if (PRINT_AMTS_MSG) then
850 write (iout,*) "acceleration/energy drift too large",amax,
851 & epdrift," split increased to ",ntime_split," itime",itime,
857 & "Uh-hu. Bumpy landscape. Maximum splitting number",
859 & " already reached!!! Trying to carry on!"
863 t_enegrad=t_enegrad+MPI_Wtime()-tt0
865 t_enegrad=t_enegrad+tcpu()-tt0
867 c Second step of the velocity Verlet algorithm
872 else if (lang.eq.3) then
874 call sd_verlet2_ciccotti
876 else if (lang.eq.1) then
881 if (rattle) call rattle2
882 c Backup the coordinates, velocities, and accelerations
886 d_t_old(j,i)=d_t(j,i)
887 d_a_old(j,i)=d_a(j,i)
894 c Restore the time step
896 c Compute long-range forces
903 call etotal_long(energia_long)
904 if (large.and. mod(itime,ntwe).eq.0)
905 & call enerprint(energia_long)
908 t_elong=t_elong+MPI_Wtime()-tt0
910 t_elong=t_elong+tcpu()-tt0
916 t_enegrad=t_enegrad+MPI_Wtime()-tt0
918 t_enegrad=t_enegrad+tcpu()-tt0
920 c Compute accelerations from long-range forces
922 if (large.and. mod(itime,ntwe).eq.0) then
923 write (iout,*) "energia_long",energia_long(0)
924 write (iout,*) "Cartesian and internal coordinates: step 2"
926 call pdbout(0.0d0,"cipiszcze",iout)
928 write (iout,*) "Accelerations from long-range forces"
930 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
931 & (d_a(j,i+nres),j=1,3)
933 write (iout,*) "Velocities, step 2"
935 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
936 & (d_t(j,i+nres),j=1,3)
940 c Compute the final RESPA step (increment velocities)
941 c write (iout,*) "*********************** RESPA fin"
943 c Compute the complete potential energy
945 potEcomp(i)=energia_short(i)+energia_long(i)
947 potE=potEcomp(0)-potEcomp(20)
948 c potE=energia_short(0)+energia_long(0)
951 c Calculate the kinetic and the total energy and the kinetic temperature
954 c Couple the system to Berendsen bath if needed
955 if (tbf .and. lang.eq.0) then
958 kinetic_T=2.0d0/(dimen3*Rb)*EK
959 c Backup the coordinates, velocities, and accelerations
961 if (mod(itime,ntwe).eq.0 .and. large) then
962 write (iout,*) "Velocities, end"
964 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
965 & (d_t(j,i+nres),j=1,3)
971 c---------------------------------------------------------------------
973 c First and last RESPA step (incrementing velocities using long-range
975 implicit real*8 (a-h,o-z)
977 include 'COMMON.CONTROL'
980 include 'COMMON.CHAIN'
981 include 'COMMON.DERIV'
983 include 'COMMON.LOCAL'
984 include 'COMMON.INTERACT'
985 include 'COMMON.IOUNITS'
986 include 'COMMON.NAMES'
988 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
992 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
996 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
999 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1005 c-----------------------------------------------------------------
1007 c Applying velocity Verlet algorithm - step 1 to coordinates
1008 implicit real*8 (a-h,o-z)
1009 include 'DIMENSIONS'
1010 include 'COMMON.CONTROL'
1011 include 'COMMON.VAR'
1013 include 'COMMON.CHAIN'
1014 include 'COMMON.DERIV'
1015 include 'COMMON.GEO'
1016 include 'COMMON.LOCAL'
1017 include 'COMMON.INTERACT'
1018 include 'COMMON.IOUNITS'
1019 include 'COMMON.NAMES'
1020 double precision adt,adt2
1023 write (iout,*) "VELVERLET1 START: DC"
1025 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1026 & (dc(j,i+nres),j=1,3)
1030 adt=d_a_old(j,0)*d_time
1032 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1033 d_t_new(j,0)=d_t_old(j,0)+adt2
1034 d_t(j,0)=d_t_old(j,0)+adt
1040 adt=d_a_old(j,i)*d_time
1042 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1043 d_t_new(j,i)=d_t_old(j,i)+adt2
1044 d_t(j,i)=d_t_old(j,i)+adt
1049 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1052 adt=d_a_old(j,inres)*d_time
1054 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1055 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1056 d_t(j,inres)=d_t_old(j,inres)+adt
1061 write (iout,*) "VELVERLET1 END: DC"
1063 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1064 & (dc(j,i+nres),j=1,3)
1069 c---------------------------------------------------------------------
1071 c Step 2 of the velocity Verlet algorithm: update velocities
1072 implicit real*8 (a-h,o-z)
1073 include 'DIMENSIONS'
1074 include 'COMMON.CONTROL'
1075 include 'COMMON.VAR'
1077 include 'COMMON.CHAIN'
1078 include 'COMMON.DERIV'
1079 include 'COMMON.GEO'
1080 include 'COMMON.LOCAL'
1081 include 'COMMON.INTERACT'
1082 include 'COMMON.IOUNITS'
1083 include 'COMMON.NAMES'
1085 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1089 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1093 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1096 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1102 c-----------------------------------------------------------------
1103 subroutine sddir_precalc
1104 c Applying velocity Verlet algorithm - step 1 to coordinates
1105 implicit real*8 (a-h,o-z)
1106 include 'DIMENSIONS'
1110 include 'COMMON.CONTROL'
1111 include 'COMMON.VAR'
1114 include 'COMMON.LANGEVIN'
1116 include 'COMMON.LANGEVIN.lang0'
1118 include 'COMMON.CHAIN'
1119 include 'COMMON.DERIV'
1120 include 'COMMON.GEO'
1121 include 'COMMON.LOCAL'
1122 include 'COMMON.INTERACT'
1123 include 'COMMON.IOUNITS'
1124 include 'COMMON.NAMES'
1125 include 'COMMON.TIME1'
1126 double precision stochforcvec(MAXRES6)
1127 common /stochcalc/ stochforcvec
1129 c Compute friction and stochastic forces
1138 time_fric=time_fric+MPI_Wtime()-time00
1141 time_fric=time_fric+tcpu()-time00
1144 call stochastic_force(stochforcvec)
1146 time_stoch=time_stoch+MPI_Wtime()-time00
1148 time_stoch=time_stoch+tcpu()-time00
1151 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1152 c forces (d_as_work)
1154 call ginv_mult(fric_work, d_af_work)
1155 call ginv_mult(stochforcvec, d_as_work)
1158 c---------------------------------------------------------------------
1159 subroutine sddir_verlet1
1160 c Applying velocity Verlet algorithm - step 1 to velocities
1161 implicit real*8 (a-h,o-z)
1162 include 'DIMENSIONS'
1163 include 'COMMON.CONTROL'
1164 include 'COMMON.VAR'
1167 include 'COMMON.LANGEVIN'
1169 include 'COMMON.LANGEVIN.lang0'
1171 include 'COMMON.CHAIN'
1172 include 'COMMON.DERIV'
1173 include 'COMMON.GEO'
1174 include 'COMMON.LOCAL'
1175 include 'COMMON.INTERACT'
1176 include 'COMMON.IOUNITS'
1177 include 'COMMON.NAMES'
1178 c Revised 3/31/05 AL: correlation between random contributions to
1179 c position and velocity increments included.
1180 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1181 double precision adt,adt2
1183 c Add the contribution from BOTH friction and stochastic force to the
1184 c coordinates, but ONLY the contribution from the friction forces to velocities
1187 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1188 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1189 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1190 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1191 d_t(j,0)=d_t_old(j,0)+adt
1196 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1197 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1198 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1199 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1200 d_t(j,i)=d_t_old(j,i)+adt
1205 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1208 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1209 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1210 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1211 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1212 d_t(j,inres)=d_t_old(j,inres)+adt
1219 c---------------------------------------------------------------------
1220 subroutine sddir_verlet2
1221 c Calculating the adjusted velocities for accelerations
1222 implicit real*8 (a-h,o-z)
1223 include 'DIMENSIONS'
1224 include 'COMMON.CONTROL'
1225 include 'COMMON.VAR'
1228 include 'COMMON.LANGEVIN'
1230 include 'COMMON.LANGEVIN.lang0'
1232 include 'COMMON.CHAIN'
1233 include 'COMMON.DERIV'
1234 include 'COMMON.GEO'
1235 include 'COMMON.LOCAL'
1236 include 'COMMON.INTERACT'
1237 include 'COMMON.IOUNITS'
1238 include 'COMMON.NAMES'
1239 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1240 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1241 c Revised 3/31/05 AL: correlation between random contributions to
1242 c position and velocity increments included.
1243 c The correlation coefficients are calculated at low-friction limit.
1244 c Also, friction forces are now not calculated with new velocities.
1246 c call friction_force
1247 call stochastic_force(stochforcvec)
1249 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1250 c forces (d_as_work)
1252 call ginv_mult(stochforcvec, d_as_work1)
1258 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1259 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1264 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1265 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1270 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1273 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1274 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1275 & +cos60*d_as_work1(ind+j))*d_time
1282 c---------------------------------------------------------------------
1283 subroutine max_accel
1285 c Find the maximum difference in the accelerations of the the sites
1286 c at the beginning and the end of the time step.
1288 implicit real*8 (a-h,o-z)
1289 include 'DIMENSIONS'
1290 include 'COMMON.CONTROL'
1291 include 'COMMON.VAR'
1293 include 'COMMON.CHAIN'
1294 include 'COMMON.DERIV'
1295 include 'COMMON.GEO'
1296 include 'COMMON.LOCAL'
1297 include 'COMMON.INTERACT'
1298 include 'COMMON.IOUNITS'
1299 double precision aux(3),accel(3),accel_old(3),dacc
1301 c aux(j)=d_a(j,0)-d_a_old(j,0)
1302 accel_old(j)=d_a_old(j,0)
1309 c 7/3/08 changed to asymmetric difference
1311 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1312 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1313 accel(j)=accel(j)+0.5d0*d_a(j,i)
1314 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1315 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1316 dacc=dabs(accel(j)-accel_old(j))
1317 c write (iout,*) i,dacc
1318 if (dacc.gt.amax) amax=dacc
1326 accel_old(j)=d_a_old(j,0)
1331 accel_old(j)=accel_old(j)+d_a_old(j,1)
1332 accel(j)=accel(j)+d_a(j,1)
1336 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1338 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1339 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1340 accel(j)=accel(j)+d_a(j,i+nres)
1344 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1345 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1346 dacc=dabs(accel(j)-accel_old(j))
1347 c write (iout,*) "side-chain",i,dacc
1348 if (dacc.gt.amax) amax=dacc
1352 accel_old(j)=accel_old(j)+d_a_old(j,i)
1353 accel(j)=accel(j)+d_a(j,i)
1354 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1359 c---------------------------------------------------------------------
1360 subroutine predict_edrift(epdrift)
1362 c Predict the drift of the potential energy
1364 implicit real*8 (a-h,o-z)
1365 include 'DIMENSIONS'
1366 include 'COMMON.CONTROL'
1367 include 'COMMON.VAR'
1369 include 'COMMON.CHAIN'
1370 include 'COMMON.DERIV'
1371 include 'COMMON.GEO'
1372 include 'COMMON.LOCAL'
1373 include 'COMMON.INTERACT'
1374 include 'COMMON.IOUNITS'
1375 include 'COMMON.MUCA'
1376 double precision epdrift,epdriftij
1377 c Drift of the potential energy
1383 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1384 if (lmuca) epdriftij=epdriftij*factor
1385 c write (iout,*) "back",i,j,epdriftij
1386 if (epdriftij.gt.epdrift) epdrift=epdriftij
1390 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1393 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1394 if (lmuca) epdriftij=epdriftij*factor
1395 c write (iout,*) "side",i,j,epdriftij
1396 if (epdriftij.gt.epdrift) epdrift=epdriftij
1400 epdrift=0.5d0*epdrift*d_time*d_time
1401 c write (iout,*) "epdrift",epdrift
1404 c-----------------------------------------------------------------------
1405 subroutine verlet_bath
1407 c Coupling to the thermostat by using the Berendsen algorithm
1409 implicit real*8 (a-h,o-z)
1410 include 'DIMENSIONS'
1411 include 'COMMON.CONTROL'
1412 include 'COMMON.VAR'
1414 include 'COMMON.CHAIN'
1415 include 'COMMON.DERIV'
1416 include 'COMMON.GEO'
1417 include 'COMMON.LOCAL'
1418 include 'COMMON.INTERACT'
1419 include 'COMMON.IOUNITS'
1420 include 'COMMON.NAMES'
1421 double precision T_half,fact
1423 T_half=2.0d0/(dimen3*Rb)*EK
1424 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1425 c write(iout,*) "T_half", T_half
1426 c write(iout,*) "EK", EK
1427 c write(iout,*) "fact", fact
1429 d_t(j,0)=fact*d_t(j,0)
1433 d_t(j,i)=fact*d_t(j,i)
1437 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1440 d_t(j,inres)=fact*d_t(j,inres)
1446 c---------------------------------------------------------
1448 c Set up the initial conditions of a MD simulation
1449 implicit real*8 (a-h,o-z)
1450 include 'DIMENSIONS'
1454 integer IERROR,ERRCODE
1456 include 'COMMON.SETUP'
1457 include 'COMMON.CONTROL'
1458 include 'COMMON.VAR'
1461 include 'COMMON.LANGEVIN'
1463 include 'COMMON.LANGEVIN.lang0'
1465 include 'COMMON.CHAIN'
1466 include 'COMMON.DERIV'
1467 include 'COMMON.GEO'
1468 include 'COMMON.LOCAL'
1469 include 'COMMON.INTERACT'
1470 include 'COMMON.IOUNITS'
1471 include 'COMMON.NAMES'
1472 include 'COMMON.REMD'
1473 real*8 energia_long(0:n_ene),
1474 & energia_short(0:n_ene),energia(0:n_ene),vcm(3),incr(3)
1475 double precision cm(3),L(3),xv,sigv,lowb,highb
1476 double precision varia(maxvar)
1484 c write(iout,*) "d_time", d_time
1485 c Compute the standard deviations of stochastic forces for Langevin dynamics
1486 c if the friction coefficients do not depend on surface area
1487 if (lang.gt.0 .and. .not.surfarea) then
1489 stdforcp(i)=stdfp*dsqrt(gamp)
1492 stdforcsc(i)=stdfsc(iabs(itype(i)))
1493 & *dsqrt(gamsc(iabs(itype(i))))
1496 c Open the pdb file for snapshotshots
1499 if (ilen(tmpdir).gt.0)
1500 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1501 & liczba(:ilen(liczba))//".pdb")
1502 #if defined(AIX) || defined(PGI)
1504 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1505 & //".pdb", position='append')
1508 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1509 & //".pdb", access='append')
1513 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1514 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1515 & liczba(:ilen(liczba))//".x")
1516 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1519 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1520 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1521 & liczba(:ilen(liczba))//".cx")
1522 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1528 if (ilen(tmpdir).gt.0)
1529 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1530 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1532 if (ilen(tmpdir).gt.0)
1533 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1534 cartname=prefix(:ilen(prefix))//"_MD.cx"
1538 write (qstr,'(256(1h ))')
1541 iq = qinfrag(i,iset)*10
1542 iw = wfrag(i,iset)/100
1544 if(me.eq.king.or..not.out1file)
1545 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1546 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1551 iq = qinpair(i,iset)*10
1552 iw = wpair(i,iset)/100
1554 if(me.eq.king.or..not.out1file)
1555 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1556 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1560 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1562 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1564 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1566 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1570 if (restart1file) then
1572 & inquire(file=mremd_rst_name,exist=file_exist)
1574 write (*,*) me," Before broadcast: file_exist",file_exist
1575 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1577 write (*,*) me," After broadcast: file_exist",file_exist
1578 c inquire(file=mremd_rst_name,exist=file_exist)
1580 if(me.eq.king.or..not.out1file)
1581 & write(iout,*) "Initial state read by master and distributed"
1583 if (ilen(tmpdir).gt.0)
1584 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1585 & //liczba(:ilen(liczba))//'.rst')
1586 inquire(file=rest2name,exist=file_exist)
1589 if(.not.restart1file) then
1590 if(me.eq.king.or..not.out1file)
1591 & write(iout,*) "Initial state will be read from file ",
1592 & rest2name(:ilen(rest2name))
1595 call rescale_weights(t_bath)
1597 if(me.eq.king.or..not.out1file)then
1598 if (restart1file) then
1599 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1602 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1605 write(iout,*) "Initial velocities randomly generated"
1612 c Generate initial velocities
1613 if(me.eq.king.or..not.out1file)
1614 & write(iout,*) "Initial velocities randomly generated"
1617 CtotTafm is the variable for AFM time which eclipsed during
1620 c rest2name = prefix(:ilen(prefix))//'.rst'
1621 if(me.eq.king.or..not.out1file)then
1622 write (iout,*) "Initial velocities"
1624 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1625 & (d_t(j,i+nres),j=1,3)
1627 c Zeroing the total angular momentum of the system
1628 write(iout,*) "Calling the zero-angular
1629 & momentum subroutine"
1632 c Getting the potential energy and forces and velocities and accelerations
1634 c write (iout,*) "velocity of the center of the mass:"
1635 c write (iout,*) (vcm(j),j=1,3)
1637 d_t(j,0)=d_t(j,0)-vcm(j)
1639 c Removing the velocity of the center of mass
1641 if(me.eq.king.or..not.out1file)then
1642 write (iout,*) "vcm right after adjustment:"
1643 write (iout,*) (vcm(j),j=1,3)
1646 if (start_from_model) then
1647 i_model=iran_num(1,constr_homology)
1648 write (iout,*) 'starting from model ',i_model
1651 c(j,i)=chomo(j,i,i_model)
1655 call int_from_cart(.true.,.false.)
1656 call sc_loc_geom(.false.)
1659 dc(j,i)=c(j,i+1)-c(j,i)
1660 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1665 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1666 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1669 call etotal(energia(0))
1671 c call chainbuild_cart
1672 write (iout,*) "Initial energies"
1673 call enerprint(energia(0))
1674 write (iout,*) "PREMINIM ",preminim
1675 if(iranconf.ne.0 .or. preminim) then
1677 write (iout,*) 'Calling OVERLAP_SC'
1678 call overlap_sc(fail)
1682 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1683 print *,'SC_move',nft_sc,etot
1684 if(me.eq.king.or..not.out1file)
1685 & write(iout,*) 'SC_move',nft_sc,etot
1689 print *, 'Calling MINIM_DC'
1690 call minim_dc(etot,iretcode,nfun)
1692 call geom_to_var(nvar,varia)
1693 print *,'Calling MINIMIZE.'
1694 call minimize(etot,varia,iretcode,nfun)
1695 call var_to_geom(nvar,varia)
1697 if(me.eq.king.or..not.out1file) then
1698 write(iout,*) "Minimized energy is",etot
1699 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1700 call etotal(potEcomp)
1701 call enerprint(potEcomp)
1704 if(iranconf.ne.0) then
1705 c 8/22/17 AL Loop to produce a low-energy random conformation
1708 if(me.eq.king.or..not.out1file)
1709 & write (iout,*) 'Calling OVERLAP_SC'
1710 call overlap_sc(fail)
1714 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1715 print *,'SC_move',nft_sc,etot
1716 if(me.eq.king.or..not.out1file)
1717 & write(iout,*) 'SC_move',nft_sc,etot
1721 print *, 'Calling MINIM_DC'
1722 call minim_dc(etot,iretcode,nfun)
1724 call geom_to_var(nvar,varia)
1725 print *,'Calling MINIMIZE.'
1726 call minimize(etot,varia,iretcode,nfun)
1727 call var_to_geom(nvar,varia)
1729 if(me.eq.king.or..not.out1file)
1730 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1732 if (isnan(etot) .or. etot.gt.1.0d4) then
1733 write (iout,*) "Energy too large",etot,
1734 & " trying another random conformation"
1737 call gen_rand_conf(itmp,*30)
1739 30 write (iout,*) 'Failed to generate random conformation',
1740 & ', itrial=',itrial
1741 write (*,*) 'Processor:',me,
1742 & ' Failed to generate random conformation',
1751 write (iout,'(a,i3,a)') 'Processor:',me,
1752 & ' error in generating random conformation.'
1753 write (*,'(a,i3,a)') 'Processor:',me,
1754 & ' error in generating random conformation.'
1757 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1767 write (iout,'(a,i3,a)') 'Processor:',me,
1768 & ' failed to generate a low-energy random conformation.'
1769 write (*,'(a,i3,a)') 'Processor:',me,
1770 & ' failed to generate a low-energy random conformation.'
1773 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1780 call chainbuild_cart
1785 kinetic_T=2.0d0/(dimen3*Rb)*EK
1786 if(me.eq.king.or..not.out1file)then
1796 call etotal(potEcomp)
1797 if (large) call enerprint(potEcomp)
1800 t_etotal=t_etotal+MPI_Wtime()-tt0
1802 t_etotal=t_etotal+tcpu()-tt0
1809 if (amax*d_time .gt. dvmax) then
1810 d_time=d_time*dvmax/amax
1811 if(me.eq.king.or..not.out1file) write (iout,*)
1812 & "Time step reduced to",d_time,
1813 & " because of too large initial acceleration."
1815 if(me.eq.king.or..not.out1file)then
1816 write(iout,*) "Potential energy and its components"
1817 call enerprint(potEcomp)
1818 c write(iout,*) (potEcomp(i),i=0,n_ene)
1820 potE=potEcomp(0)-potEcomp(20)
1823 if (ntwe.ne.0) call statout(itime)
1824 if(me.eq.king.or..not.out1file)
1825 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1826 & " Kinetic energy",EK," Potential energy",potE,
1827 & " Total energy",totE," Maximum acceleration ",
1830 write (iout,*) "Initial coordinates"
1832 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1833 & (c(j,i+nres),j=1,3)
1835 write (iout,*) "Initial dC"
1837 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1838 & (dc(j,i+nres),j=1,3)
1840 write (iout,*) "Initial velocities"
1841 write (iout,"(13x,' backbone ',23x,' side chain')")
1843 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1844 & (d_t(j,i+nres),j=1,3)
1846 write (iout,*) "Initial accelerations"
1848 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1849 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1850 & (d_a(j,i+nres),j=1,3)
1856 d_t_old(j,i)=d_t(j,i)
1857 d_a_old(j,i)=d_a(j,i)
1859 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1868 call etotal_short(energia_short)
1869 if (large) call enerprint(potEcomp)
1872 t_eshort=t_eshort+MPI_Wtime()-tt0
1874 t_eshort=t_eshort+tcpu()-tt0
1879 if(.not.out1file .and. large) then
1880 write (iout,*) "energia_long",energia_long(0),
1881 & " energia_short",energia_short(0),
1882 & " total",energia_long(0)+energia_short(0)
1883 write (iout,*) "Initial fast-force accelerations"
1885 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1886 & (d_a(j,i+nres),j=1,3)
1889 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1892 d_a_short(j,i)=d_a(j,i)
1901 call etotal_long(energia_long)
1902 if (large) call enerprint(potEcomp)
1905 t_elong=t_elong+MPI_Wtime()-tt0
1907 t_elong=t_elong+tcpu()-tt0
1912 if(.not.out1file .and. large) then
1913 write (iout,*) "energia_long",energia_long(0)
1914 write (iout,*) "Initial slow-force accelerations"
1916 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1917 & (d_a(j,i+nres),j=1,3)
1921 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1923 t_enegrad=t_enegrad+tcpu()-tt0
1928 c-----------------------------------------------------------
1929 subroutine random_vel
1930 implicit real*8 (a-h,o-z)
1931 include 'DIMENSIONS'
1932 include 'COMMON.CONTROL'
1933 include 'COMMON.VAR'
1936 include 'COMMON.LANGEVIN'
1938 include 'COMMON.LANGEVIN.lang0'
1940 include 'COMMON.CHAIN'
1941 include 'COMMON.DERIV'
1942 include 'COMMON.GEO'
1943 include 'COMMON.LOCAL'
1944 include 'COMMON.INTERACT'
1945 include 'COMMON.IOUNITS'
1946 include 'COMMON.NAMES'
1947 include 'COMMON.TIME1'
1948 double precision xv,sigv,lowb,highb,vec_afm(3)
1949 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1950 c First generate velocities in the eigenspace of the G matrix
1951 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1958 sigv=dsqrt((Rb*t_bath)/geigen(i))
1961 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1963 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1964 c & " d_t_work_new",d_t_work_new(ii)
1967 C if (SELFGUIDE.gt.0) then
1970 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
1971 C distance=distance+vec_afm(j)**2
1973 C distance=dsqrt(distance)
1975 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
1976 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
1977 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
1978 C & d_t_work_new(j+(afmend-1)*3)
1989 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1992 c write (iout,*) "Ek from eigenvectors",Ek1
1994 c Transform velocities to UNRES coordinate space
2000 d_t_work(ind)=d_t_work(ind)
2001 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2003 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2007 c Transfer to the d_t vector
2009 d_t(j,0)=d_t_work(j)
2015 d_t(j,i)=d_t_work(ind)
2019 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2022 d_t(j,i+nres)=d_t_work(ind)
2027 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2028 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2033 c-----------------------------------------------------------
2034 subroutine sd_verlet_p_setup
2035 c Sets up the parameters of stochastic Verlet algorithm
2036 implicit real*8 (a-h,o-z)
2037 include 'DIMENSIONS'
2041 include 'COMMON.CONTROL'
2042 include 'COMMON.VAR'
2045 include 'COMMON.LANGEVIN'
2047 include 'COMMON.LANGEVIN.lang0'
2049 include 'COMMON.CHAIN'
2050 include 'COMMON.DERIV'
2051 include 'COMMON.GEO'
2052 include 'COMMON.LOCAL'
2053 include 'COMMON.INTERACT'
2054 include 'COMMON.IOUNITS'
2055 include 'COMMON.NAMES'
2056 include 'COMMON.TIME1'
2057 double precision emgdt(MAXRES6),
2058 & pterm,vterm,rho,rhoc,vsig,
2059 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2060 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2061 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2062 logical lprn /.false./
2063 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2064 double precision ktm
2071 c AL 8/17/04 Code adapted from tinker
2073 c Get the frictional and random terms for stochastic dynamics in the
2074 c eigenspace of mass-scaled UNRES friction matrix
2077 gdt = fricgam(i) * d_time
2079 c Stochastic dynamics reduces to simple MD for zero friction
2081 if (gdt .le. zero) then
2082 pfric_vec(i) = 1.0d0
2083 vfric_vec(i) = d_time
2084 afric_vec(i) = 0.5d0 * d_time * d_time
2085 prand_vec(i) = 0.0d0
2086 vrand_vec1(i) = 0.0d0
2087 vrand_vec2(i) = 0.0d0
2089 c Analytical expressions when friction coefficient is large
2092 if (gdt .ge. gdt_radius) then
2095 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2096 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2097 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2098 vterm = 1.0d0 - egdt**2
2099 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2101 c Use series expansions when friction coefficient is small
2112 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2113 & - gdt5/120.0d0 + gdt6/720.0d0
2114 & - gdt7/5040.0d0 + gdt8/40320.0d0
2115 & - gdt9/362880.0d0) / fricgam(i)**2
2116 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2117 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2118 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2119 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2120 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2121 & + 127.0d0*gdt9/90720.0d0
2122 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2123 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2124 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2125 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2126 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2127 & - 17.0d0*gdt2/1280.0d0
2128 & + 17.0d0*gdt3/6144.0d0
2129 & + 40967.0d0*gdt4/34406400.0d0
2130 & - 57203.0d0*gdt5/275251200.0d0
2131 & - 1429487.0d0*gdt6/13212057600.0d0)
2134 c Compute the scaling factors of random terms for the nonzero friction case
2136 ktm = 0.5d0*d_time/fricgam(i)
2137 psig = dsqrt(ktm*pterm) / fricgam(i)
2138 vsig = dsqrt(ktm*vterm)
2139 rhoc = dsqrt(1.0d0 - rho*rho)
2141 vrand_vec1(i) = vsig * rho
2142 vrand_vec2(i) = vsig * rhoc
2147 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2150 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2151 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2155 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2158 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2159 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2160 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2161 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2162 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2163 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2166 t_sdsetup=t_sdsetup+MPI_Wtime()
2168 t_sdsetup=t_sdsetup+tcpu()-tt0
2172 c-------------------------------------------------------------
2173 subroutine eigtransf1(n,ndim,ab,d,c)
2176 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2182 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2188 c-------------------------------------------------------------
2189 subroutine eigtransf(n,ndim,a,b,d,c)
2192 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2198 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2204 c-------------------------------------------------------------
2205 subroutine sd_verlet1
2206 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2207 implicit real*8 (a-h,o-z)
2208 include 'DIMENSIONS'
2209 include 'COMMON.CONTROL'
2210 include 'COMMON.VAR'
2213 include 'COMMON.LANGEVIN'
2215 include 'COMMON.LANGEVIN.lang0'
2217 include 'COMMON.CHAIN'
2218 include 'COMMON.DERIV'
2219 include 'COMMON.GEO'
2220 include 'COMMON.LOCAL'
2221 include 'COMMON.INTERACT'
2222 include 'COMMON.IOUNITS'
2223 include 'COMMON.NAMES'
2224 double precision stochforcvec(MAXRES6)
2225 common /stochcalc/ stochforcvec
2226 logical lprn /.false./
2228 c write (iout,*) "dc_old"
2230 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2231 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2234 dc_work(j)=dc_old(j,0)
2235 d_t_work(j)=d_t_old(j,0)
2236 d_a_work(j)=d_a_old(j,0)
2241 dc_work(ind+j)=dc_old(j,i)
2242 d_t_work(ind+j)=d_t_old(j,i)
2243 d_a_work(ind+j)=d_a_old(j,i)
2248 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2250 dc_work(ind+j)=dc_old(j,i+nres)
2251 d_t_work(ind+j)=d_t_old(j,i+nres)
2252 d_a_work(ind+j)=d_a_old(j,i+nres)
2260 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2264 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2265 & vfric_mat(i,j),afric_mat(i,j),
2266 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2274 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2275 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2276 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2277 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2279 d_t_work_new(i)=ddt1+0.5d0*ddt2
2280 d_t_work(i)=ddt1+ddt2
2285 d_t(j,0)=d_t_work(j)
2290 dc(j,i)=dc_work(ind+j)
2291 d_t(j,i)=d_t_work(ind+j)
2296 if (itype(i).ne.10) then
2299 dc(j,inres)=dc_work(ind+j)
2300 d_t(j,inres)=d_t_work(ind+j)
2307 c--------------------------------------------------------------------------
2308 subroutine sd_verlet2
2309 c Calculating the adjusted velocities for accelerations
2310 implicit real*8 (a-h,o-z)
2311 include 'DIMENSIONS'
2312 include 'COMMON.CONTROL'
2313 include 'COMMON.VAR'
2316 include 'COMMON.LANGEVIN'
2318 include 'COMMON.LANGEVIN.lang0'
2320 include 'COMMON.CHAIN'
2321 include 'COMMON.DERIV'
2322 include 'COMMON.GEO'
2323 include 'COMMON.LOCAL'
2324 include 'COMMON.INTERACT'
2325 include 'COMMON.IOUNITS'
2326 include 'COMMON.NAMES'
2327 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2328 common /stochcalc/ stochforcvec
2330 c Compute the stochastic forces which contribute to velocity change
2332 call stochastic_force(stochforcvecV)
2339 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2340 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2341 & vrand_mat2(i,j)*stochforcvecV(j)
2343 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2347 d_t(j,0)=d_t_work(j)
2352 d_t(j,i)=d_t_work(ind+j)
2357 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2360 d_t(j,inres)=d_t_work(ind+j)
2367 c-----------------------------------------------------------
2368 subroutine sd_verlet_ciccotti_setup
2369 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2371 implicit real*8 (a-h,o-z)
2372 include 'DIMENSIONS'
2376 include 'COMMON.CONTROL'
2377 include 'COMMON.VAR'
2380 include 'COMMON.LANGEVIN'
2382 include 'COMMON.LANGEVIN.lang0'
2384 include 'COMMON.CHAIN'
2385 include 'COMMON.DERIV'
2386 include 'COMMON.GEO'
2387 include 'COMMON.LOCAL'
2388 include 'COMMON.INTERACT'
2389 include 'COMMON.IOUNITS'
2390 include 'COMMON.NAMES'
2391 include 'COMMON.TIME1'
2392 double precision emgdt(MAXRES6),
2393 & pterm,vterm,rho,rhoc,vsig,
2394 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2395 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2396 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2397 logical lprn /.false./
2398 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2399 double precision ktm
2406 c AL 8/17/04 Code adapted from tinker
2408 c Get the frictional and random terms for stochastic dynamics in the
2409 c eigenspace of mass-scaled UNRES friction matrix
2412 write (iout,*) "i",i," fricgam",fricgam(i)
2413 gdt = fricgam(i) * d_time
2415 c Stochastic dynamics reduces to simple MD for zero friction
2417 if (gdt .le. zero) then
2418 pfric_vec(i) = 1.0d0
2419 vfric_vec(i) = d_time
2420 afric_vec(i) = 0.5d0*d_time*d_time
2421 prand_vec(i) = afric_vec(i)
2422 vrand_vec2(i) = vfric_vec(i)
2424 c Analytical expressions when friction coefficient is large
2429 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2430 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2431 prand_vec(i) = afric_vec(i)
2432 vrand_vec2(i) = vfric_vec(i)
2434 c Compute the scaling factors of random terms for the nonzero friction case
2436 c ktm = 0.5d0*d_time/fricgam(i)
2437 c psig = dsqrt(ktm*pterm) / fricgam(i)
2438 c vsig = dsqrt(ktm*vterm)
2439 c prand_vec(i) = psig*afric_vec(i)
2440 c vrand_vec2(i) = vsig*vfric_vec(i)
2445 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2448 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2449 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2453 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2455 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2456 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2457 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2458 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2459 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2461 t_sdsetup=t_sdsetup+MPI_Wtime()
2463 t_sdsetup=t_sdsetup+tcpu()-tt0
2467 c-------------------------------------------------------------
2468 subroutine sd_verlet1_ciccotti
2469 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2470 implicit real*8 (a-h,o-z)
2471 include 'DIMENSIONS'
2475 include 'COMMON.CONTROL'
2476 include 'COMMON.VAR'
2479 include 'COMMON.LANGEVIN'
2481 include 'COMMON.LANGEVIN.lang0'
2483 include 'COMMON.CHAIN'
2484 include 'COMMON.DERIV'
2485 include 'COMMON.GEO'
2486 include 'COMMON.LOCAL'
2487 include 'COMMON.INTERACT'
2488 include 'COMMON.IOUNITS'
2489 include 'COMMON.NAMES'
2490 double precision stochforcvec(MAXRES6)
2491 common /stochcalc/ stochforcvec
2492 logical lprn /.false./
2494 c write (iout,*) "dc_old"
2496 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2497 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2500 dc_work(j)=dc_old(j,0)
2501 d_t_work(j)=d_t_old(j,0)
2502 d_a_work(j)=d_a_old(j,0)
2507 dc_work(ind+j)=dc_old(j,i)
2508 d_t_work(ind+j)=d_t_old(j,i)
2509 d_a_work(ind+j)=d_a_old(j,i)
2514 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2516 dc_work(ind+j)=dc_old(j,i+nres)
2517 d_t_work(ind+j)=d_t_old(j,i+nres)
2518 d_a_work(ind+j)=d_a_old(j,i+nres)
2527 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2531 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2532 & vfric_mat(i,j),afric_mat(i,j),
2533 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2541 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2542 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2543 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2544 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2546 d_t_work_new(i)=ddt1+0.5d0*ddt2
2547 d_t_work(i)=ddt1+ddt2
2552 d_t(j,0)=d_t_work(j)
2557 dc(j,i)=dc_work(ind+j)
2558 d_t(j,i)=d_t_work(ind+j)
2563 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2566 dc(j,inres)=dc_work(ind+j)
2567 d_t(j,inres)=d_t_work(ind+j)
2574 c--------------------------------------------------------------------------
2575 subroutine sd_verlet2_ciccotti
2576 c Calculating the adjusted velocities for accelerations
2577 implicit real*8 (a-h,o-z)
2578 include 'DIMENSIONS'
2579 include 'COMMON.CONTROL'
2580 include 'COMMON.VAR'
2583 include 'COMMON.LANGEVIN'
2585 include 'COMMON.LANGEVIN.lang0'
2587 include 'COMMON.CHAIN'
2588 include 'COMMON.DERIV'
2589 include 'COMMON.GEO'
2590 include 'COMMON.LOCAL'
2591 include 'COMMON.INTERACT'
2592 include 'COMMON.IOUNITS'
2593 include 'COMMON.NAMES'
2594 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2595 common /stochcalc/ stochforcvec
2597 c Compute the stochastic forces which contribute to velocity change
2599 call stochastic_force(stochforcvecV)
2606 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2607 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2608 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2610 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2614 d_t(j,0)=d_t_work(j)
2619 d_t(j,i)=d_t_work(ind+j)
2624 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2627 d_t(j,inres)=d_t_work(ind+j)