2 c------------------------------------------------
3 c The driver for molecular dynamics subroutines
4 c------------------------------------------------
11 include 'COMMON.SETUP'
12 include 'COMMON.CONTROL'
16 include 'COMMON.LAGRANGE.5diag'
18 include 'COMMON.LAGRANGE'
21 include 'COMMON.LANGEVIN'
24 include 'COMMON.LANGEVIN.lang0.5diag'
26 include 'COMMON.LANGEVIN.lang0'
29 include 'COMMON.CHAIN'
30 include 'COMMON.DERIV'
32 include 'COMMON.LOCAL'
33 include 'COMMON.INTERACT'
34 include 'COMMON.IOUNITS'
35 include 'COMMON.NAMES'
36 include 'COMMON.TIME1'
37 include 'COMMON.HAIRPIN'
38 double precision cm(3),L(3),vcm(3)
40 double precision v_work(maxres6),v_transf(maxres6)
48 integer i,j,icount_scale,itime_scal
49 integer nharp,iharp(4,maxres/3)
50 double precision scalfac
54 if (ilen(tmpdir).gt.0)
55 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
56 & //liczba(:ilen(liczba))//'.rst')
58 if (ilen(tmpdir).gt.0)
59 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
61 write (iout,*) "MD lang",lang
67 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
73 c Determine the inverse of the inertia matrix.
74 call setup_MD_matrices
78 t_MDsetup = MPI_Wtime()-tt0
80 t_MDsetup = tcpu()-tt0
83 c Entering the MD loop
89 if (lang.eq.2 .or. lang.eq.3) then
93 call sd_verlet_p_setup
95 call sd_verlet_ciccotti_setup
99 pfric0_mat(i,j,0)=pfric_mat(i,j)
100 afric0_mat(i,j,0)=afric_mat(i,j)
101 vfric0_mat(i,j,0)=vfric_mat(i,j)
102 prand0_mat(i,j,0)=prand_mat(i,j)
103 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
104 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
109 flag_stoch(i)=.false.
113 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
115 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
119 else if (lang.eq.1 .or. lang.eq.4) then
123 t_langsetup=MPI_Wtime()-tt0
126 t_langsetup=tcpu()-tt0
129 do itime=1,n_timestep
131 if (large.and. mod(itime,ntwe).eq.0)
132 & write (iout,*) "itime",itime
134 if (lang.gt.0 .and. surfarea .and.
135 & mod(itime,reset_fricmat).eq.0) then
136 if (lang.eq.2 .or. lang.eq.3) then
140 call sd_verlet_p_setup
142 call sd_verlet_ciccotti_setup
146 pfric0_mat(i,j,0)=pfric_mat(i,j)
147 afric0_mat(i,j,0)=afric_mat(i,j)
148 vfric0_mat(i,j,0)=vfric_mat(i,j)
149 prand0_mat(i,j,0)=prand_mat(i,j)
150 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
151 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
156 flag_stoch(i)=.false.
159 else if (lang.eq.1 .or. lang.eq.4) then
162 write (iout,'(a,i10)')
163 & "Friction matrix reset based on surface area, itime",itime
165 if (reset_vel .and. tbf .and. lang.eq.0
166 & .and. mod(itime,count_reset_vel).eq.0) then
168 write(iout,'(a,f20.2)')
169 & "Velocities reset to random values, time",totT
172 d_t_old(j,i)=d_t(j,i)
176 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
180 d_t(j,0)=d_t(j,0)-vcm(j)
183 kinetic_T=2.0d0/(dimen3*Rb)*EK
184 scalfac=dsqrt(T_bath/kinetic_T)
185 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
188 d_t_old(j,i)=scalfac*d_t(j,i)
194 c Time-reversible RESPA algorithm
195 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
196 call RESPA_step(itime)
198 c Variable time step algorithm.
199 call velverlet_step(itime)
203 call brown_step(itime)
205 print *,"Brown dynamics not here!"
207 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
214 if (mod(itime,ntwe).eq.0) then
216 C call enerprint(potEcomp)
217 C print *,itime,'AFM',Eafmforc,etot
231 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
234 v_work(ind)=d_t(j,i+nres)
239 write (66,'(80f10.5)')
240 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
244 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
246 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
248 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
251 if (mod(itime,ntwx).eq.0) then
252 c write(iout,*) 'time=',itime
253 C call check_ecartint
255 write (tytul,'("time",f8.2)') totT
257 call hairpin(.true.,nharp,iharp)
258 call secondary2(.true.)
259 call pdbout(potE,tytul,ipdb)
264 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
265 open(irest2,file=rest2name,status='unknown')
266 write(irest2,*) totT,EK,potE,totE,t_bath
268 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
271 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
282 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
284 & 'MD calculations setup:',t_MDsetup,
285 & 'Energy & gradient evaluation:',t_enegrad,
286 & 'Stochastic MD setup:',t_langsetup,
287 & 'Stochastic MD step setup:',t_sdsetup,
289 write (iout,'(/28(1h=),a25,27(1h=))')
290 & ' End of MD calculation '
292 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
294 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
295 & " time_fricmatmult",time_fricmatmult," time_fsample ",
300 c-------------------------------------------------------------------------------
301 subroutine velverlet_step(itime)
302 c-------------------------------------------------------------------------------
303 c Perform a single velocity Verlet step; the time step can be rescaled if
304 c increments in accelerations exceed the threshold
305 c-------------------------------------------------------------------------------
310 integer ierror,ierrcode,errcode
312 include 'COMMON.SETUP'
313 include 'COMMON.CONTROL'
317 include 'COMMON.LAGRANGE.5diag'
319 include 'COMMON.LAGRANGE'
322 include 'COMMON.LANGEVIN'
325 include 'COMMON.LANGEVIN.lang0.5diag'
327 include 'COMMON.LANGEVIN.lang0'
330 include 'COMMON.CHAIN'
331 include 'COMMON.DERIV'
333 include 'COMMON.LOCAL'
334 include 'COMMON.INTERACT'
335 include 'COMMON.IOUNITS'
336 include 'COMMON.NAMES'
337 include 'COMMON.TIME1'
338 include 'COMMON.MUCA'
339 double precision vcm(3),incr(3)
340 double precision cm(3),L(3)
341 integer ilen,count,rstcount
344 integer maxcount_scale /20/
346 double precision stochforcvec(MAXRES6)
347 common /stochcalc/ stochforcvec
348 integer itime,icount_scale,itime_scal,ifac_time,i,j,itt
350 double precision epdrift,fac_time
357 else if (lang.eq.2 .or. lang.eq.3) then
359 call stochastic_force(stochforcvec)
362 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
364 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
371 icount_scale=icount_scale+1
372 if (icount_scale.gt.maxcount_scale) then
374 & "ERROR: too many attempts at scaling down the time step. ",
375 & "amax=",amax,"epdrift=",epdrift,
376 & "damax=",damax,"edriftmax=",edriftmax,
380 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
384 c First step of the velocity Verlet algorithm
389 else if (lang.eq.3) then
391 call sd_verlet1_ciccotti
393 else if (lang.eq.1) then
398 c Build the chain from the newly calculated coordinates
400 if (rattle) call rattle1
402 if (large.and. mod(itime,ntwe).eq.0) then
403 write (iout,*) "Cartesian and internal coordinates: step 1"
408 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
409 & (dc(j,i+nres),j=1,3)
411 write (iout,*) "Accelerations"
413 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
414 & (d_a(j,i+nres),j=1,3)
416 write (iout,*) "Velocities, step 1"
418 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
419 & (d_t(j,i+nres),j=1,3)
428 c Calculate energy and forces
430 call etotal(potEcomp)
431 ! AL 4/17/17: Reduce the steps if NaNs occurred.
432 if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
437 if (large.and. mod(itime,ntwe).eq.0)
438 & call enerprint(potEcomp)
441 t_etotal=t_etotal+MPI_Wtime()-tt0
443 t_etotal=t_etotal+tcpu()-tt0
446 potE=potEcomp(0)-potEcomp(27)
448 c Get the new accelerations
451 t_enegrad=t_enegrad+MPI_Wtime()-tt0
453 t_enegrad=t_enegrad+tcpu()-tt0
455 c Determine maximum acceleration and scale down the timestep if needed
457 amax=amax/(itime_scal+1)**2
458 call predict_edrift(epdrift)
459 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
460 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
462 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
464 itime_scal=itime_scal+ifac_time
465 c fac_time=dmin1(damax/amax,0.5d0)
466 fac_time=0.5d0**ifac_time
467 d_time=d_time*fac_time
468 if (lang.eq.2 .or. lang.eq.3) then
470 c write (iout,*) "Calling sd_verlet_setup: 1"
471 c Rescale the stochastic forces and recalculate or restore
472 c the matrices of tinker integrator
473 if (itime_scal.gt.maxflag_stoch) then
474 if (large) write (iout,'(a,i5,a)')
475 & "Calculate matrices for stochastic step;",
476 & " itime_scal ",itime_scal
478 call sd_verlet_p_setup
480 call sd_verlet_ciccotti_setup
482 write (iout,'(2a,i3,a,i3,1h.)')
483 & "Warning: cannot store matrices for stochastic",
484 & " integration because the index",itime_scal,
485 & " is greater than",maxflag_stoch
486 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
487 & " integration Langevin algorithm for better efficiency."
488 else if (flag_stoch(itime_scal)) then
489 if (large) write (iout,'(a,i5,a,l1)')
490 & "Restore matrices for stochastic step; itime_scal ",
491 & itime_scal," flag ",flag_stoch(itime_scal)
494 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
495 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
496 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
497 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
498 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
499 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
503 if (large) write (iout,'(2a,i5,a,l1)')
504 & "Calculate & store matrices for stochastic step;",
505 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
507 call sd_verlet_p_setup
509 call sd_verlet_ciccotti_setup
511 flag_stoch(ifac_time)=.true.
514 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
515 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
516 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
517 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
518 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
519 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
523 fac_time=1.0d0/dsqrt(fac_time)
525 stochforcvec(i)=fac_time*stochforcvec(i)
528 else if (lang.eq.1) then
529 c Rescale the accelerations due to stochastic forces
530 fac_time=1.0d0/dsqrt(fac_time)
532 d_as_work(i)=d_as_work(i)*fac_time
535 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
536 & "itime",itime," Timestep scaled down to ",
537 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
539 c Second step of the velocity Verlet algorithm
544 else if (lang.eq.3) then
546 call sd_verlet2_ciccotti
548 else if (lang.eq.1) then
553 if (rattle) call rattle2
556 C print *,totTafm,"TU?"
557 if (d_time.ne.d_time0) then
560 if (lang.eq.2 .or. lang.eq.3) then
561 if (large) write (iout,'(a)')
562 & "Restore original matrices for stochastic step"
563 c write (iout,*) "Calling sd_verlet_setup: 2"
564 c Restore the matrices of tinker integrator if the time step has been restored
567 pfric_mat(i,j)=pfric0_mat(i,j,0)
568 afric_mat(i,j)=afric0_mat(i,j,0)
569 vfric_mat(i,j)=vfric0_mat(i,j,0)
570 prand_mat(i,j)=prand0_mat(i,j,0)
571 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
572 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
581 c Calculate the kinetic and the total energy and the kinetic temperature
586 c write (iout,*) "step",itime," EK",EK," EK1",EK1
588 c Couple the system to Berendsen bath if needed
589 if (tbf .and. lang.eq.0) then
592 kinetic_T=2.0d0/(dimen3*Rb)*EK
593 c Backup the coordinates, velocities, and accelerations
597 d_t_old(j,i)=d_t(j,i)
598 d_a_old(j,i)=d_a(j,i)
602 if (mod(itime,ntwe).eq.0 .and. large) then
603 write (iout,*) "Velocities, step 2"
605 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
606 & (d_t(j,i+nres),j=1,3)
612 c-------------------------------------------------------------------------------
613 subroutine RESPA_step(itime)
614 c-------------------------------------------------------------------------------
615 c Perform a single RESPA step.
616 c-------------------------------------------------------------------------------
621 integer IERROR,ERRCODE
623 include 'COMMON.SETUP'
624 include 'COMMON.CONTROL'
628 include 'COMMON.LAGRANGE.5diag'
630 include 'COMMON.LAGRANGE'
633 include 'COMMON.LANGEVIN'
636 include 'COMMON.LANGEVIN.lang0.5diag'
638 include 'COMMON.LANGEVIN.lang0'
641 include 'COMMON.CHAIN'
642 include 'COMMON.DERIV'
644 include 'COMMON.LOCAL'
645 include 'COMMON.INTERACT'
646 include 'COMMON.IOUNITS'
647 include 'COMMON.NAMES'
648 include 'COMMON.TIME1'
650 common /shortcheck/ lprint_short
651 double precision energia_short(0:n_ene),
652 & energia_long(0:n_ene)
653 double precision cm(3),L(3),vcm(3),incr(3)
654 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
655 & d_a_old0(3,0:maxres2)
657 double precision fac_time
658 logical PRINT_AMTS_MSG /.false./
659 integer ilen,count,rstcount
662 integer maxcount_scale /10/
664 double precision stochforcvec(MAXRES6)
665 common /stochcalc/ stochforcvec
669 common /cipiszcze/ itt
671 double precision epdrift,epdriftmax
675 if (large.and. mod(itime,ntwe).eq.0) then
676 write (iout,*) "***************** RESPA itime",itime
677 write (iout,*) "Cartesian and internal coordinates: step 0"
682 write (iout,*) "Accelerations from long-range forces"
684 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
685 & (d_a(j,i+nres),j=1,3)
687 write (iout,*) "Velocities, step 0"
689 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
690 & (d_t(j,i+nres),j=1,3)
695 c Perform the initial RESPA step (increment velocities)
696 c write (iout,*) "*********************** RESPA ini"
699 if (mod(itime,ntwe).eq.0 .and. large) then
700 write (iout,*) "Velocities, end"
702 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
703 & (d_t(j,i+nres),j=1,3)
707 c Compute the short-range forces
713 C 7/2/2009 commented out
715 c call etotal_short(energia_short)
718 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
721 d_a(j,i)=d_a_short(j,i)
725 if (large.and. mod(itime,ntwe).eq.0) then
726 write (iout,*) "energia_short",energia_short(0)
727 write (iout,*) "Accelerations from short-range forces"
729 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
730 & (d_a(j,i+nres),j=1,3)
735 t_enegrad=t_enegrad+MPI_Wtime()-tt0
737 t_enegrad=t_enegrad+tcpu()-tt0
742 d_t_old(j,i)=d_t(j,i)
743 d_a_old(j,i)=d_a(j,i)
746 c 6/30/08 A-MTS: attempt at increasing the split number
749 dc_old0(j,i)=dc_old(j,i)
750 d_t_old0(j,i)=d_t_old(j,i)
751 d_a_old0(j,i)=d_a_old(j,i)
754 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
755 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
762 c write (iout,*) "itime",itime," ntime_split",ntime_split
763 c Split the time step
764 d_time=d_time0/ntime_split
765 c Perform the short-range RESPA steps (velocity Verlet increments of
766 c positions and velocities using short-range forces)
767 c write (iout,*) "*********************** RESPA split"
768 do itsplit=1,ntime_split
771 else if (lang.eq.2 .or. lang.eq.3) then
773 call stochastic_force(stochforcvec)
776 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
778 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
783 c First step of the velocity Verlet algorithm
788 else if (lang.eq.3) then
790 call sd_verlet1_ciccotti
792 else if (lang.eq.1) then
797 c Build the chain from the newly calculated coordinates
799 if (rattle) call rattle1
801 if (large.and. mod(itime,ntwe).eq.0) then
802 write (iout,*) "***** ITSPLIT",itsplit
803 write (iout,*) "Cartesian and internal coordinates: step 1"
808 write (iout,*) "Velocities, step 1"
810 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
811 & (d_t(j,i+nres),j=1,3)
820 c Calculate energy and forces
821 c if (large.and. mod(itime,ntwe).eq.0) lprint_short=.true.
823 call etotal_short(energia_short)
824 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
825 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) )
827 if (PRINT_AMTS_MSG) write (iout,*)
828 & "Infinities/NaNs in energia_short",
829 & energia_short(0),"; increasing ntime_split to",ntime_split
830 ntime_split=ntime_split*2
831 if (ntime_split.gt.maxtime_split) then
833 write (iout,*) "Cannot rescue the run; aborting job.",
834 & " Retry with a smaller time step"
836 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
838 write (iout,*) "Cannot rescue the run; terminating.",
839 & " Retry with a smaller time step"
845 if (large.and. mod(itime,ntwe).eq.0)
846 & call enerprint(energia_short)
849 t_eshort=t_eshort+MPI_Wtime()-tt0
851 t_eshort=t_eshort+tcpu()-tt0
855 c Get the new accelerations
857 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
860 d_a_short(j,i)=d_a(j,i)
864 if (large.and. mod(itime,ntwe).eq.0) then
865 write (iout,*)"energia_short",energia_short(0)
866 write (iout,*) "Accelerations from short-range forces"
868 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
869 & (d_a(j,i+nres),j=1,3)
874 c Determine maximum acceleration and scale down the timestep if needed
876 amax=amax/ntime_split**2
877 call predict_edrift(epdrift)
878 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
879 & write (iout,*) "amax",amax," damax",damax,
880 & " epdrift",epdrift," epdriftmax",epdriftmax
881 c Exit loop and try with increased split number if the change of
882 c acceleration is too big
883 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
884 if (ntime_split.lt.maxtime_split) then
886 ntime_split=ntime_split*2
889 dc_old(j,i)=dc_old0(j,i)
890 d_t_old(j,i)=d_t_old0(j,i)
891 d_a_old(j,i)=d_a_old0(j,i)
894 if (PRINT_AMTS_MSG) then
895 write (iout,*) "acceleration/energy drift too large",amax,
896 & epdrift," split increased to ",ntime_split," itime",itime,
902 & "Uh-hu. Bumpy landscape. Maximum splitting number",
904 & " already reached!!! Trying to carry on!"
908 t_enegrad=t_enegrad+MPI_Wtime()-tt0
910 t_enegrad=t_enegrad+tcpu()-tt0
912 c Second step of the velocity Verlet algorithm
917 else if (lang.eq.3) then
919 call sd_verlet2_ciccotti
921 else if (lang.eq.1) then
926 if (rattle) call rattle2
927 c Backup the coordinates, velocities, and accelerations
931 d_t_old(j,i)=d_t(j,i)
932 d_a_old(j,i)=d_a(j,i)
939 c Restore the time step
941 c Compute long-range forces
948 call etotal_long(energia_long)
949 ! AL 4/17/2017 Handling NaNs
950 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
953 & "Infinitied/NaNs in energia_long, Aborting MPI job."
954 call enerprint(energia_long)
956 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
958 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
959 call enerprint(energia_long)
964 c lprint_short=.false.
965 if (large.and. mod(itime,ntwe).eq.0)
966 & call enerprint(energia_long)
969 t_elong=t_elong+MPI_Wtime()-tt0
971 t_elong=t_elong+tcpu()-tt0
977 t_enegrad=t_enegrad+MPI_Wtime()-tt0
979 t_enegrad=t_enegrad+tcpu()-tt0
981 c Compute accelerations from long-range forces
983 if (large.and. mod(itime,ntwe).eq.0) then
984 write (iout,*) "energia_long",energia_long(0)
985 write (iout,*) "Cartesian and internal coordinates: step 2"
990 write (iout,*) "Accelerations from long-range forces"
992 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
993 & (d_a(j,i+nres),j=1,3)
995 write (iout,*) "Velocities, step 2"
997 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
998 & (d_t(j,i+nres),j=1,3)
1002 c Compute the final RESPA step (increment velocities)
1003 c write (iout,*) "*********************** RESPA fin"
1005 c Compute the complete potential energy
1007 potEcomp(i)=energia_short(i)+energia_long(i)
1009 potE=potEcomp(0)-potEcomp(27)
1011 if (large.and. mod(itime,ntwe).eq.0) then
1012 call enerprint(potEcomp)
1013 write (iout,*) "potE",potE
1016 c potE=energia_short(0)+energia_long(0)
1019 c Calculate the kinetic and the total energy and the kinetic temperature
1022 c Couple the system to Berendsen bath if needed
1023 if (tbf .and. lang.eq.0) then
1026 kinetic_T=2.0d0/(dimen3*Rb)*EK
1027 c Backup the coordinates, velocities, and accelerations
1029 if (mod(itime,ntwe).eq.0 .and. large) then
1030 write (iout,*) "Velocities, end"
1032 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1033 & (d_t(j,i+nres),j=1,3)
1039 c---------------------------------------------------------------------
1040 subroutine RESPA_vel
1041 c First and last RESPA step (incrementing velocities using long-range
1044 include 'DIMENSIONS'
1045 include 'COMMON.CONTROL'
1046 include 'COMMON.VAR'
1049 include 'COMMON.LAGRANGE.5diag'
1051 include 'COMMON.LAGRANGE'
1053 include 'COMMON.CHAIN'
1054 include 'COMMON.DERIV'
1055 include 'COMMON.GEO'
1056 include 'COMMON.LOCAL'
1057 include 'COMMON.INTERACT'
1058 include 'COMMON.IOUNITS'
1059 include 'COMMON.NAMES'
1062 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1066 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1070 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1073 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1079 c-----------------------------------------------------------------
1081 c Applying velocity Verlet algorithm - step 1 to coordinates
1083 include 'DIMENSIONS'
1084 include 'COMMON.CONTROL'
1085 include 'COMMON.VAR'
1088 include 'COMMON.LAGRANGE.5diag'
1090 include 'COMMON.LAGRANGE'
1092 include 'COMMON.CHAIN'
1093 include 'COMMON.DERIV'
1094 include 'COMMON.GEO'
1095 include 'COMMON.LOCAL'
1096 include 'COMMON.INTERACT'
1097 include 'COMMON.IOUNITS'
1098 include 'COMMON.NAMES'
1099 double precision adt,adt2
1102 write (iout,*) "VELVERLET1 START: DC"
1104 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1105 & (dc(j,i+nres),j=1,3)
1109 adt=d_a_old(j,0)*d_time
1111 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1112 d_t_new(j,0)=d_t_old(j,0)+adt2
1113 d_t(j,0)=d_t_old(j,0)+adt
1119 write (iout,*) "i",i," d_a_old",(d_a_old(j,i),j=1,3)
1122 adt=d_a_old(j,i)*d_time
1124 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1125 d_t_new(j,i)=d_t_old(j,i)+adt2
1126 d_t(j,i)=d_t_old(j,i)+adt
1131 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1134 adt=d_a_old(j,inres)*d_time
1136 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1137 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1138 d_t(j,inres)=d_t_old(j,inres)+adt
1143 write (iout,*) "VELVERLET1 END: DC"
1145 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1146 & (dc(j,i+nres),j=1,3)
1151 c---------------------------------------------------------------------
1153 c Step 2 of the velocity Verlet algorithm: update velocities
1155 include 'DIMENSIONS'
1156 include 'COMMON.CONTROL'
1157 include 'COMMON.VAR'
1160 include 'COMMON.LAGRANGE.5diag'
1162 include 'COMMON.LAGRANGE'
1164 include 'COMMON.CHAIN'
1165 include 'COMMON.DERIV'
1166 include 'COMMON.GEO'
1167 include 'COMMON.LOCAL'
1168 include 'COMMON.INTERACT'
1169 include 'COMMON.IOUNITS'
1170 include 'COMMON.NAMES'
1173 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1177 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1181 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1184 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1190 c-----------------------------------------------------------------
1191 subroutine sddir_precalc
1192 c Applying velocity Verlet algorithm - step 1 to coordinates
1194 include 'DIMENSIONS'
1198 include 'COMMON.CONTROL'
1199 include 'COMMON.VAR'
1202 include 'COMMON.LAGRANGE.5diag'
1204 include 'COMMON.LAGRANGE'
1207 include 'COMMON.LANGEVIN'
1210 include 'COMMON.LANGEVIN.lang0.5diag'
1212 include 'COMMON.LANGEVIN.lang0'
1215 include 'COMMON.CHAIN'
1216 include 'COMMON.DERIV'
1217 include 'COMMON.GEO'
1218 include 'COMMON.LOCAL'
1219 include 'COMMON.INTERACT'
1220 include 'COMMON.IOUNITS'
1221 include 'COMMON.NAMES'
1222 include 'COMMON.TIME1'
1223 double precision time00
1224 double precision stochforcvec(MAXRES6)
1225 common /stochcalc/ stochforcvec
1228 c Compute friction and stochastic forces
1236 c write (iout,*) "After friction_force"
1238 time_fric=time_fric+MPI_Wtime()-time00
1241 time_fric=time_fric+tcpu()-time00
1244 call stochastic_force(stochforcvec)
1245 c write (iout,*) "After stochastic_force"
1247 time_stoch=time_stoch+MPI_Wtime()-time00
1249 time_stoch=time_stoch+tcpu()-time00
1252 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1253 c forces (d_as_work)
1256 c write (iout,*) "friction accelerations"
1257 call fivediaginv_mult(dimen,fric_work, d_af_work)
1258 c write (iout,*) "stochastic acceleratios"
1259 call fivediaginv_mult(dimen,stochforcvec, d_as_work)
1261 call ginv_mult(fric_work, d_af_work)
1262 call ginv_mult(stochforcvec, d_as_work)
1265 write (iout,*) "d_af_work"
1266 write (iout,'(3f10.5)') (d_af_work(i),i=1,dimen3)
1267 write (iout,*) "d_as_work"
1268 write (iout,'(3f10.5)') (d_as_work(i),i=1,dimen3)
1269 write (iout,*) "Leaving sddir_precalc"
1273 c---------------------------------------------------------------------
1274 subroutine sddir_verlet1
1275 c Applying velocity Verlet algorithm - step 1 to velocities
1277 include 'DIMENSIONS'
1278 include 'COMMON.CONTROL'
1279 include 'COMMON.VAR'
1282 include 'COMMON.LAGRANGE.5diag'
1284 include 'COMMON.LAGRANGE'
1287 include 'COMMON.LANGEVIN'
1290 include 'COMMON.LANGEVIN.lang0.5diag'
1292 include 'COMMON.LANGEVIN.lang0'
1295 include 'COMMON.CHAIN'
1296 include 'COMMON.DERIV'
1297 include 'COMMON.GEO'
1298 include 'COMMON.LOCAL'
1299 include 'COMMON.INTERACT'
1300 include 'COMMON.IOUNITS'
1301 include 'COMMON.NAMES'
1302 c Revised 3/31/05 AL: correlation between random contributions to
1303 c position and velocity increments included.
1304 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1305 double precision adt,adt2
1306 integer i,j,ind,inres
1308 c Add the contribution from BOTH friction and stochastic force to the
1309 c coordinates, but ONLY the contribution from the friction forces to velocities
1312 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1313 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1314 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1315 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1316 d_t(j,0)=d_t_old(j,0)+adt
1321 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1322 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1323 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1324 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1325 d_t(j,i)=d_t_old(j,i)+adt
1330 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1333 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1334 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1335 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1336 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1337 d_t(j,inres)=d_t_old(j,inres)+adt
1344 c---------------------------------------------------------------------
1345 subroutine sddir_verlet2
1346 c Calculating the adjusted velocities for accelerations
1348 include 'DIMENSIONS'
1349 include 'COMMON.CONTROL'
1350 include 'COMMON.VAR'
1353 include 'COMMON.LAGRANGE.5diag'
1355 include 'COMMON.LAGRANGE'
1358 include 'COMMON.LANGEVIN'
1361 include 'COMMON.LANGEVIN.lang0.5diag'
1363 include 'COMMON.LANGEVIN.lang0'
1366 include 'COMMON.CHAIN'
1367 include 'COMMON.DERIV'
1368 include 'COMMON.GEO'
1369 include 'COMMON.LOCAL'
1370 include 'COMMON.INTERACT'
1371 include 'COMMON.IOUNITS'
1372 include 'COMMON.NAMES'
1373 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1374 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1375 integer i,j,inres,ind
1376 c Revised 3/31/05 AL: correlation between random contributions to
1377 c position and velocity increments included.
1378 c The correlation coefficients are calculated at low-friction limit.
1379 c Also, friction forces are now not calculated with new velocities.
1381 c call friction_force
1382 call stochastic_force(stochforcvec)
1384 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1385 c forces (d_as_work)
1388 call fivediaginv_mult(maxres6,stochforcvec, d_as_work1)
1390 call ginv_mult(stochforcvec, d_as_work1)
1396 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1397 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1402 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1403 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1408 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1411 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1412 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1413 & +cos60*d_as_work1(ind+j))*d_time
1420 c---------------------------------------------------------------------
1421 subroutine max_accel
1423 c Find the maximum difference in the accelerations of the the sites
1424 c at the beginning and the end of the time step.
1427 include 'DIMENSIONS'
1428 include 'COMMON.CONTROL'
1429 include 'COMMON.VAR'
1432 include 'COMMON.LAGRANGE.5diag'
1434 include 'COMMON.LAGRANGE'
1436 include 'COMMON.CHAIN'
1437 include 'COMMON.DERIV'
1438 include 'COMMON.GEO'
1439 include 'COMMON.LOCAL'
1440 include 'COMMON.INTERACT'
1441 include 'COMMON.IOUNITS'
1443 double precision aux(3),accel(3),accel_old(3),dacc
1445 c aux(j)=d_a(j,0)-d_a_old(j,0)
1446 accel_old(j)=d_a_old(j,0)
1453 c 7/3/08 changed to asymmetric difference
1455 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1456 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1457 accel(j)=accel(j)+0.5d0*d_a(j,i)
1458 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1459 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1460 dacc=dabs(accel(j)-accel_old(j))
1461 c write (iout,*) i,dacc
1462 if (dacc.gt.amax) amax=dacc
1470 accel_old(j)=d_a_old(j,0)
1475 accel_old(j)=accel_old(j)+d_a_old(j,1)
1476 accel(j)=accel(j)+d_a(j,1)
1480 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1482 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1483 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1484 accel(j)=accel(j)+d_a(j,i+nres)
1488 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1489 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1490 dacc=dabs(accel(j)-accel_old(j))
1491 c write (iout,*) "side-chain",i,dacc
1492 if (dacc.gt.amax) amax=dacc
1496 accel_old(j)=accel_old(j)+d_a_old(j,i)
1497 accel(j)=accel(j)+d_a(j,i)
1498 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1503 c---------------------------------------------------------------------
1504 subroutine predict_edrift(epdrift)
1506 c Predict the drift of the potential energy
1509 include 'DIMENSIONS'
1510 include 'COMMON.CONTROL'
1511 include 'COMMON.VAR'
1514 include 'COMMON.LAGRANGE.5diag'
1516 include 'COMMON.LAGRANGE'
1518 include 'COMMON.CHAIN'
1519 include 'COMMON.DERIV'
1520 include 'COMMON.GEO'
1521 include 'COMMON.LOCAL'
1522 include 'COMMON.INTERACT'
1523 include 'COMMON.IOUNITS'
1524 include 'COMMON.MUCA'
1525 double precision epdrift,epdriftij
1527 c Drift of the potential energy
1533 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1534 if (lmuca) epdriftij=epdriftij*factor
1535 c write (iout,*) "back",i,j,epdriftij
1536 if (epdriftij.gt.epdrift) epdrift=epdriftij
1540 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1543 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1544 if (lmuca) epdriftij=epdriftij*factor
1545 c write (iout,*) "side",i,j,epdriftij
1546 if (epdriftij.gt.epdrift) epdrift=epdriftij
1550 epdrift=0.5d0*epdrift*d_time*d_time
1551 c write (iout,*) "epdrift",epdrift
1554 c-----------------------------------------------------------------------
1555 subroutine verlet_bath
1557 c Coupling to the thermostat by using the Berendsen algorithm
1560 include 'DIMENSIONS'
1561 include 'COMMON.CONTROL'
1562 include 'COMMON.VAR'
1565 include 'COMMON.LAGRANGE.5diag'
1567 include 'COMMON.LAGRANGE'
1571 include 'COMMON.LANGEVIN.lang0.5diag'
1573 include 'COMMON.LANGEVIN.lang0'
1576 include 'COMMON.LANGEVIN'
1578 include 'COMMON.CHAIN'
1579 include 'COMMON.DERIV'
1580 include 'COMMON.GEO'
1581 include 'COMMON.LOCAL'
1582 include 'COMMON.INTERACT'
1583 include 'COMMON.IOUNITS'
1584 include 'COMMON.NAMES'
1585 double precision T_half,fact
1587 T_half=2.0d0/(dimen3*Rb)*EK
1588 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1589 c write(iout,*) "T_half", T_half
1590 c write(iout,*) "EK", EK
1591 c write(iout,*) "fact", fact
1593 d_t(j,0)=fact*d_t(j,0)
1597 d_t(j,i)=fact*d_t(j,i)
1601 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1604 d_t(j,inres)=fact*d_t(j,inres)
1610 c---------------------------------------------------------
1612 c Set up the initial conditions of a MD simulation
1614 include 'DIMENSIONS'
1618 integer IERROR,ERRCODE,error_msg,ierr,ierrcode
1620 include 'COMMON.SETUP'
1621 include 'COMMON.CONTROL'
1622 include 'COMMON.VAR'
1625 include 'COMMON.LAGRANGE.5diag'
1627 include 'COMMON.LAGRANGE'
1629 include 'COMMON.QRESTR'
1631 include 'COMMON.LANGEVIN'
1634 include 'COMMON.LANGEVIN.lang0.5diag'
1637 include 'COMMON.LANGEVIN.lang0.5diag'
1639 include 'COMMON.LANGEVIN.lang0'
1643 include 'COMMON.CHAIN'
1644 include 'COMMON.DERIV'
1645 include 'COMMON.GEO'
1646 include 'COMMON.LOCAL'
1647 include 'COMMON.INTERACT'
1648 include 'COMMON.IOUNITS'
1649 include 'COMMON.NAMES'
1650 include 'COMMON.REMD'
1651 include 'COMMON.TIME1'
1655 common /lbfgstat/ status,niter,nfun
1657 integer n_model_try,list_model_try(max_template),k
1658 double precision tt0
1659 real*8 energia_long(0:n_ene),
1660 & energia_short(0:n_ene),vcm(3),incr(3)
1661 double precision cm(3),L(3),xv,sigv,lowb,highb
1662 double precision varia(maxvar),energia(0:n_ene)
1669 integer i,ipos,iq,iw,j,iranmin,nft_sc,iretcode,nfun,itrial,itmp,
1672 double precision etot
1674 integer i_start_models(0:nodes-1)
1675 write (iout,*) "init_MD INDPDB",indpdb
1677 c write(iout,*) "d_time", d_time
1678 c Compute the standard deviations of stochastic forces for Langevin dynamics
1679 c if the friction coefficients do not depend on surface area
1680 if (lang.gt.0 .and. .not.surfarea) then
1682 stdforcp(i)=stdfp*dsqrt(gamp)
1685 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1686 & *dsqrt(gamsc(iabs(itype(i))))
1689 c Open the pdb file for snapshotshots
1692 if (ilen(tmpdir).gt.0)
1693 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1694 & liczba(:ilen(liczba))//".pdb")
1696 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1700 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1701 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1702 & liczba(:ilen(liczba))//".x")
1703 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1706 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1707 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1708 & liczba(:ilen(liczba))//".cx")
1709 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1715 if (ilen(tmpdir).gt.0)
1716 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1717 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1719 if (ilen(tmpdir).gt.0)
1720 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1721 cartname=prefix(:ilen(prefix))//"_MD.cx"
1725 write (qstr,'(256(1h ))')
1728 iq = qinfrag(i,iset)*10
1729 iw = wfrag(i,iset)/100
1731 if(me.eq.king.or..not.out1file)
1732 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1733 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1738 iq = qinpair(i,iset)*10
1739 iw = wpair(i,iset)/100
1741 if(me.eq.king.or..not.out1file)
1742 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1743 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1747 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1749 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1751 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1753 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1756 write (iout,*) "REST ",rest
1758 if (restart1file) then
1760 & inquire(file=mremd_rst_name,exist=file_exist)
1762 c write (*,*) me," Before broadcast: file_exist",file_exist
1763 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1765 c write (*,*) me," After broadcast: file_exist",file_exist
1766 c inquire(file=mremd_rst_name,exist=file_exist)
1768 if(me.eq.king.or..not.out1file)
1769 & write(iout,*) "Initial state read by master and distributed"
1771 if (ilen(tmpdir).gt.0)
1772 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1773 & //liczba(:ilen(liczba))//'.rst')
1774 inquire(file=rest2name,exist=file_exist)
1777 if(.not.restart1file) then
1778 if(me.eq.king.or..not.out1file)
1779 & write(iout,*) "Initial state will be read from file ",
1780 & rest2name(:ilen(rest2name))
1783 call rescale_weights(t_bath)
1786 if(me.eq.king.or..not.out1file)then
1787 if (restart1file) then
1788 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1791 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1794 write(iout,*) "Initial velocities randomly generated"
1801 c Generate initial velocities
1802 if(me.eq.king.or..not.out1file)
1803 & write(iout,*) "Initial velocities randomly generated"
1806 CtotTafm is the variable for AFM time which eclipsed during
1809 c rest2name = prefix(:ilen(prefix))//'.rst'
1810 if(me.eq.king.or..not.out1file)then
1811 write (iout,*) "Initial velocities"
1813 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1814 & (d_t(j,i+nres),j=1,3)
1818 c Zeroing the total angular momentum of the system
1819 write(iout,*) "Calling the zero-angular
1820 & momentum subroutine"
1822 c Getting the potential energy and forces and velocities and accelerations
1824 write (iout,*) "velocity of the center of the mass:"
1825 write (iout,*) (vcm(j),j=1,3)
1828 d_t(j,0)=d_t(j,0)-vcm(j)
1830 c Removing the velocity of the center of mass
1832 if(me.eq.king.or..not.out1file)then
1833 write (iout,*) "vcm right after adjustment:"
1834 write (iout,*) (vcm(j),j=1,3)
1835 write (iout,*) "Initial velocities after adjustment"
1837 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1838 & (d_t(j,i+nres),j=1,3)
1843 c write (iout,*) "init_MD before initial structure REST ",rest
1844 if(start_from_model .and. (me.eq.king .or. .not. out1file))
1845 & write(iout,*) 'START_FROM_MODELS is ON'
1846 if(start_from_model .and. rest) then
1847 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1849 & 'START_FROM_MODELS is OFF because the run is restarted'
1850 write(iout,*) 'Remove restart keyword from input'
1853 c write (iout,*) "rest ",rest," start_from_model",start_from_model,
1854 c & " nmodel_start",nmodel_start," preminim",preminim
1857 if (iranconf.ne.0) then
1858 c 8/22/17 AL Loop to produce a low-energy random conformation
1862 print *, 'Calling OVERLAP_SC'
1863 call overlap_sc(fail)
1867 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1868 print *,'SC_move',nft_sc,etot
1869 if(me.eq.king.or..not.out1file)
1870 & write(iout,*) 'SC_move',nft_sc,etot
1874 if (me.eq.king.or..not.out1file) write(iout,*)
1875 & 'Minimizing random structure: Calling MINIM_DC'
1876 call minim_dc(etot,iretcode,nfun)
1878 call geom_to_var(nvar,varia)
1879 if (me.eq.king.or..not.out1file) write (iout,*)
1880 & 'Minimizing random structure: Calling MINIMIZE.'
1881 call minimize(etot,varia,iretcode,nfun)
1882 call var_to_geom(nvar,varia)
1884 if (me.eq.king.or..not.out1file)
1885 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1886 if (isnan(etot) .or. etot.gt.1.0d4) then
1887 write (iout,*) "Energy too large",etot,
1888 & " trying another random conformation"
1891 call gen_rand_conf(itmp,*30)
1893 30 write (iout,*) 'Failed to generate random conformation',
1894 & ', itrial=',itrial
1895 write (*,*) 'Processor:',me,
1896 & ' Failed to generate random conformation',
1906 write (iout,'(a,i3,a)') 'Processor:',me,
1907 & ' error in generating random conformation.'
1908 write (*,'(a,i3,a)') 'Processor:',me,
1909 & ' error in generating random conformation.'
1912 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1921 write (iout,'(a,i3,a)') 'Processor:',me,
1922 & ' failed to generate a low-energy random conformation.'
1923 write (*,'(a,i3,a)') 'Processor:',me,
1924 & ' failed to generate a low-energy random conformation.'
1927 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1932 else if (preminim) then
1933 if (start_from_model) then
1937 do while (fail .and. n_model_try.lt.nmodel_start)
1938 write (iout,*) "n_model_try",n_model_try
1940 i_model=iran_num(1,nmodel_start)
1942 if (i_model.eq.list_model_try(k)) exit
1944 if (k.gt.n_model_try) exit
1946 n_model_try=n_model_try+1
1947 list_model_try(n_model_try)=i_model
1948 if (me.eq.king .or. .not. out1file)
1949 & write (iout,*) 'Trying to start from model ',
1950 & pdbfiles_chomo(i_model)(:ilen(pdbfiles_chomo(i_model)))
1953 c(j,i)=chomo(j,i,i_model)
1956 call int_from_cart(.true.,.false.)
1957 call sc_loc_geom(.false.)
1961 dc(j,i)=c(j,i+1)-c(j,i)
1962 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1967 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1968 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1971 if (me.eq.king.or..not.out1file) then
1972 write (iout,*) "Energies before removing overlaps"
1973 call etotal(energia(0))
1974 call enerprint(energia(0))
1976 ! Remove SC overlaps if requested
1978 write (iout,*) 'Calling OVERLAP_SC'
1979 call overlap_sc(fail)
1982 & "Failed to remove overlap from model",i_model
1986 if (me.eq.king.or..not.out1file) then
1987 write (iout,*) "Energies after removing overlaps"
1988 call etotal(energia(0))
1989 call enerprint(energia(0))
1992 ! Search for better SC rotamers if requested
1994 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1995 print *,'SC_move',nft_sc,etot
1996 if (me.eq.king.or..not.out1file)
1997 & write(iout,*) 'SC_move',nft_sc,etot
1999 call etotal(energia(0))
2002 call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),
2003 & 1,MPI_INTEGER,king,CG_COMM,IERROR)
2004 if (n_model_try.gt.nmodel_start .and.
2005 & (me.eq.king .or. out1file)) then
2007 & "All models have irreparable overlaps. Trying randoms starts."
2009 i_model=nmodel_start+1
2013 ! Remove SC overlaps if requested
2015 write (iout,*) 'Calling OVERLAP_SC'
2016 call overlap_sc(fail)
2019 & "Failed to remove overlap"
2022 if (me.eq.king.or..not.out1file) then
2023 write (iout,*) "Energies after removing overlaps"
2024 call etotal(energia(0))
2025 call enerprint(energia(0))
2028 C 8/22/17 AL Minimize initial structure
2030 if (me.eq.king.or..not.out1file) write(iout,*)
2031 & 'Minimizing initial PDB structure: Calling MINIM_DC'
2032 call minim_dc(etot,iretcode,nfun)
2034 call geom_to_var(nvar,varia)
2035 if(me.eq.king.or..not.out1file) write (iout,*)
2036 & 'Minimizing initial PDB structure: Calling MINIMIZE.'
2037 call minimize(etot,varia,iretcode,nfun)
2038 call var_to_geom(nvar,varia)
2040 if (me.eq.king.or..not.out1file)
2041 & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2042 if(me.eq.king.or..not.out1file)
2043 & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2045 if (me.eq.king.or..not.out1file)
2046 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2047 if(me.eq.king.or..not.out1file)
2048 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2052 if (nmodel_start.gt.0 .and. me.eq.king) then
2053 write (iout,'(a)') "Task Starting model"
2055 if (i_start_models(i).gt.nmodel_start) then
2056 write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE"
2058 write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i))
2059 & (:ilen(pdbfiles_chomo(i_start_models(i))))
2064 call chainbuild_cart
2069 kinetic_T=2.0d0/(dimen3*Rb)*EK
2070 write (iout,*) "Initial kinetic energy",EK," kinetic T",kinetic_T
2071 if(me.eq.king.or..not.out1file)then
2081 call etotal(potEcomp)
2082 if (large) call enerprint(potEcomp)
2085 t_etotal=t_etotal+MPI_Wtime()-tt0
2087 t_etotal=t_etotal+tcpu()-tt0
2091 c write (iout,*) "PotE-homology",potE-potEcomp(27)
2095 if (amax*d_time .gt. dvmax) then
2096 d_time=d_time*dvmax/amax
2097 if(me.eq.king.or..not.out1file) write (iout,*)
2098 & "Time step reduced to",d_time,
2099 & " because of too large initial acceleration."
2101 if(me.eq.king.or..not.out1file)then
2102 write(iout,*) "Potential energy and its components"
2103 call enerprint(potEcomp)
2104 c write(iout,*) (potEcomp(i),i=0,n_ene)
2106 potE=potEcomp(0)-potEcomp(27)
2107 c write (iout,*) "PotE-homology",potE
2111 if (ntwe.ne.0) call statout(itime)
2112 if(me.eq.king.or..not.out1file)
2113 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
2114 & " Kinetic energy",EK," Potential energy",potE,
2115 & " Total energy",totE," Maximum acceleration ",
2118 write (iout,*) "Initial coordinates"
2120 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
2121 & (c(j,i+nres),j=1,3)
2123 write (iout,*) "Initial dC"
2125 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
2126 & (dc(j,i+nres),j=1,3)
2128 write (iout,*) "Initial velocities"
2129 write (iout,"(13x,' backbone ',23x,' side chain')")
2131 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
2132 & (d_t(j,i+nres),j=1,3)
2134 write (iout,*) "Initial accelerations"
2136 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2137 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
2138 & (d_a(j,i+nres),j=1,3)
2144 d_t_old(j,i)=d_t(j,i)
2145 d_a_old(j,i)=d_a(j,i)
2147 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2156 call etotal_short(energia_short)
2157 if (large) call enerprint(potEcomp)
2160 t_eshort=t_eshort+MPI_Wtime()-tt0
2162 t_eshort=t_eshort+tcpu()-tt0
2167 if(.not.out1file .and. large) then
2168 write (iout,*) "energia_long",energia_long(0),
2169 & " energia_short",energia_short(0),
2170 & " total",energia_long(0)+energia_short(0)
2171 write (iout,*) "Initial fast-force accelerations"
2173 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2174 & (d_a(j,i+nres),j=1,3)
2177 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2180 d_a_short(j,i)=d_a(j,i)
2189 call etotal_long(energia_long)
2190 if (large) call enerprint(potEcomp)
2193 t_elong=t_elong+MPI_Wtime()-tt0
2195 t_elong=t_elong+tcpu()-tt0
2200 if(.not.out1file .and. large) then
2201 write (iout,*) "energia_long",energia_long(0)
2202 write (iout,*) "Initial slow-force accelerations"
2204 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2205 & (d_a(j,i+nres),j=1,3)
2209 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2211 t_enegrad=t_enegrad+tcpu()-tt0
2216 c-----------------------------------------------------------
2217 subroutine random_vel
2219 include 'DIMENSIONS'
2220 include 'COMMON.CONTROL'
2221 include 'COMMON.VAR'
2224 include 'COMMON.LAGRANGE.5diag'
2226 include 'COMMON.LAGRANGE'
2229 include 'COMMON.LANGEVIN'
2232 include 'COMMON.LANGEVIN.lang0.5diag'
2234 include 'COMMON.LANGEVIN.lang0'
2237 include 'COMMON.CHAIN'
2238 include 'COMMON.DERIV'
2239 include 'COMMON.GEO'
2240 include 'COMMON.LOCAL'
2241 include 'COMMON.INTERACT'
2242 include 'COMMON.IOUNITS'
2243 include 'COMMON.NAMES'
2244 include 'COMMON.TIME1'
2245 double precision xv,sigv,lowb,highb,vec_afm(3),Ek1,Ek2,Ek3,aux
2246 integer i,ii,j,k,l,ind
2247 double precision anorm_distr
2248 logical lprn /.false./
2250 integer ichain,n,innt,inct,ibeg,ierr
2251 double precision work(8*maxres6)
2252 integer iwork(maxres6)
2253 double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),
2254 & Gvec(maxres2_chain,maxres2_chain)
2255 common /przechowalnia/Ghalf,Geigen,Gvec
2257 double precision inertia(maxres2_chain,maxres2_chain)
2259 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2260 c First generate velocities in the eigenspace of the G matrix
2261 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2264 write (iout,*) "Random_vel, fivediag"
2273 n=dimen_chain(ichain)
2274 innt=iposd_chain(ichain)
2277 write (iout,*) "Chain",ichain," n",n," start",innt
2279 if (i.lt.inct-1) then
2280 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),
2282 else if (i.eq.inct-1) then
2283 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
2285 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
2289 ghalf(ind+1)=dmorig(innt)
2290 ghalf(ind+2)=du1orig(innt)
2291 ghalf(ind+3)=dmorig(innt+1)
2295 c write (iout,*) "i",i," ind",ind," indu2",innt+i-2,
2296 c & " indu1",innt+i-1," indm",innt+i
2297 ghalf(ind+1)=du2orig(innt-1+i-2)
2298 ghalf(ind+2)=du1orig(innt-1+i-1)
2299 ghalf(ind+3)=dmorig(innt-1+i)
2300 c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
2301 c & "DU1",innt-1+i-1,"DM ",innt-1+i
2309 inertia(i,j)=ghalf(ind)
2310 inertia(j,i)=ghalf(ind)
2315 write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
2316 write (iout,*) "Five-diagonal inertia matrix, lower triangle"
2317 call matoutr(n,ghalf)
2319 call gldiag(maxres2_chain,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
2321 write (iout,'(//a,i3)')
2322 & "Eigenvectors and eigenvalues of the G matrix chain",ichain
2323 call eigout(n,n,maxres2_chain,maxres2_chain,Gvec,Geigen)
2326 c check diagonalization
2332 aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
2336 write (iout,*) i,j,aux,geigen(i)
2338 write (iout,*) i,j,aux
2348 sigv=dsqrt((Rb*t_bath)/geigen(i))
2351 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2352 EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
2353 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2354 c & " d_t_work_new",d_t_work_new(ii)
2362 d_t_work(ind)=d_t_work(ind)
2363 & +Gvec(i,j)*d_t_work_new((j-1)*3+k)
2365 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2374 aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
2380 c Transfer to the d_t vector
2381 innt=chain_border(1,ichain)
2382 inct=chain_border(2,ichain)
2384 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
2388 d_t(j,i)=d_t_work(ind)
2390 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2393 d_t(j,i+nres)=d_t_work(ind)
2400 write (iout,*) "Random velocities in the Calpha,SC space"
2402 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2403 & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2406 call kinetic_CASC(Ek1)
2408 ! Transform the velocities to virtual-bond space
2416 if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
2418 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2422 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2423 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2436 c write (iout,*) "Shifting accelerations"
2438 c write (iout,*) "ichain",chain_border1(1,ichain)-1,
2439 c & chain_border1(1,ichain)
2440 d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
2441 d_t(:,chain_border1(1,ichain))=0.0d0
2443 c write (iout,*) "Adding accelerations"
2445 c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2446 c & chain_border(2,ichain-1)
2447 d_t(:,chain_border1(1,ichain)-1)=
2448 & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2449 d_t(:,chain_border(2,ichain-1))=0.0d0
2452 write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2453 & chain_border(2,ichain-1)
2454 d_t(:,chain_border1(1,ichain)-1)=
2455 & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2456 d_t(:,chain_border(2,ichain-1))=0.0d0
2461 c d_t(j,0)=d_t(j,nnt)
2464 innt=chain_border(1,ichain)
2465 inct=chain_border(2,ichain)
2466 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
2467 c write (iout,*) "ibeg",ibeg
2469 d_t(j,ibeg)=d_t(j,innt)
2473 if (iabs(itype(i).eq.10)) then
2474 c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
2476 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2480 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2481 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2490 & "Random velocities in the virtual-bond-vector space"
2491 write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
2493 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2494 & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2497 write (iout,*) "Kinetic energy from inertia matrix eigenvalues",
2500 & "Kinetic temperatures from inertia matrix eigenvalues",
2503 write (iout,*) "Kinetic energy from inertia matrix",Ek3
2504 write (iout,*) "Kinetic temperatures from inertia",
2505 & 2*Ek3/(3*dimen*Rb)
2507 write (iout,*) "Kinetic energy from velocities in CA-SC space",
2510 & "Kinetic temperatures from velovities in CA-SC space",
2511 & 2*Ek1/(3*dimen*Rb)
2514 & "Kinetic energy from virtual-bond-vector velocities",Ek1
2516 & "Kinetic temperature from virtual-bond-vector velocities ",
2525 sigv=dsqrt((Rb*t_bath)/geigen(i))
2528 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2530 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2531 c & " d_t_work_new",d_t_work_new(ii)
2534 C if (SELFGUIDE.gt.0) then
2537 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
2538 C distance=distance+vec_afm(j)**2
2540 C distance=dsqrt(distance)
2542 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
2543 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
2544 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
2545 C & d_t_work_new(j+(afmend-1)*3)
2556 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2559 c write (iout,*) "Ek from eigenvectors",Ek1
2561 c Transform velocities to UNRES coordinate space
2567 d_t_work(ind)=d_t_work(ind)
2568 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2570 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2574 c Transfer to the d_t vector
2576 d_t(j,0)=d_t_work(j)
2582 d_t(j,i)=d_t_work(ind)
2586 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2589 d_t(j,i+nres)=d_t_work(ind)
2594 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2595 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2601 c-----------------------------------------------------------
2602 subroutine sd_verlet_p_setup
2603 c Sets up the parameters of stochastic Verlet algorithm
2605 include 'DIMENSIONS'
2609 include 'COMMON.CONTROL'
2610 include 'COMMON.VAR'
2613 include 'COMMON.LANGEVIN'
2616 include 'COMMON.LANGEVIN.lang0.5diag'
2618 include 'COMMON.LANGEVIN.lang0'
2621 include 'COMMON.CHAIN'
2622 include 'COMMON.DERIV'
2623 include 'COMMON.GEO'
2624 include 'COMMON.LOCAL'
2625 include 'COMMON.INTERACT'
2626 include 'COMMON.IOUNITS'
2627 include 'COMMON.NAMES'
2628 include 'COMMON.TIME1'
2629 double precision emgdt(MAXRES6),
2630 & pterm,vterm,rho,rhoc,vsig,
2631 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2632 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2633 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2634 logical lprn /.false./
2635 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2636 double precision ktm
2643 c AL 8/17/04 Code adapted from tinker
2645 c Get the frictional and random terms for stochastic dynamics in the
2646 c eigenspace of mass-scaled UNRES friction matrix
2649 gdt = fricgam(i) * d_time
2651 c Stochastic dynamics reduces to simple MD for zero friction
2653 if (gdt .le. zero) then
2654 pfric_vec(i) = 1.0d0
2655 vfric_vec(i) = d_time
2656 afric_vec(i) = 0.5d0 * d_time * d_time
2657 prand_vec(i) = 0.0d0
2658 vrand_vec1(i) = 0.0d0
2659 vrand_vec2(i) = 0.0d0
2661 c Analytical expressions when friction coefficient is large
2664 if (gdt .ge. gdt_radius) then
2667 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2668 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2669 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2670 vterm = 1.0d0 - egdt**2
2671 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2673 c Use series expansions when friction coefficient is small
2684 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2685 & - gdt5/120.0d0 + gdt6/720.0d0
2686 & - gdt7/5040.0d0 + gdt8/40320.0d0
2687 & - gdt9/362880.0d0) / fricgam(i)**2
2688 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2689 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2690 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2691 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2692 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2693 & + 127.0d0*gdt9/90720.0d0
2694 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2695 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2696 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2697 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2698 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2699 & - 17.0d0*gdt2/1280.0d0
2700 & + 17.0d0*gdt3/6144.0d0
2701 & + 40967.0d0*gdt4/34406400.0d0
2702 & - 57203.0d0*gdt5/275251200.0d0
2703 & - 1429487.0d0*gdt6/13212057600.0d0)
2706 c Compute the scaling factors of random terms for the nonzero friction case
2708 ktm = 0.5d0*d_time/fricgam(i)
2709 psig = dsqrt(ktm*pterm) / fricgam(i)
2710 vsig = dsqrt(ktm*vterm)
2711 rhoc = dsqrt(1.0d0 - rho*rho)
2713 vrand_vec1(i) = vsig * rho
2714 vrand_vec2(i) = vsig * rhoc
2719 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2722 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2723 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2727 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2729 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2730 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2731 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2732 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2733 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2734 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2736 t_sdsetup=t_sdsetup+MPI_Wtime()
2738 t_sdsetup=t_sdsetup+tcpu()-tt0
2742 c-------------------------------------------------------------
2743 subroutine eigtransf1(n,ndim,ab,d,c)
2746 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2752 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2758 c-------------------------------------------------------------
2759 subroutine eigtransf(n,ndim,a,b,d,c)
2762 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2768 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2774 c-------------------------------------------------------------
2775 subroutine sd_verlet1
2776 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2778 include 'DIMENSIONS'
2779 include 'COMMON.CONTROL'
2780 include 'COMMON.VAR'
2783 include 'COMMON.LAGRANGE.5diag'
2785 include 'COMMON.LAGRANGE'
2788 include 'COMMON.LANGEVIN'
2791 include 'COMMON.LANGEVIN.lang0.5diag'
2793 include 'COMMON.LANGEVIN.lang0'
2796 include 'COMMON.CHAIN'
2797 include 'COMMON.DERIV'
2798 include 'COMMON.GEO'
2799 include 'COMMON.LOCAL'
2800 include 'COMMON.INTERACT'
2801 include 'COMMON.IOUNITS'
2802 include 'COMMON.NAMES'
2803 double precision stochforcvec(MAXRES6)
2804 common /stochcalc/ stochforcvec
2805 logical lprn /.false./
2806 integer i,j,ind,inres
2808 c write (iout,*) "dc_old"
2810 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2811 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2814 dc_work(j)=dc_old(j,0)
2815 d_t_work(j)=d_t_old(j,0)
2816 d_a_work(j)=d_a_old(j,0)
2821 dc_work(ind+j)=dc_old(j,i)
2822 d_t_work(ind+j)=d_t_old(j,i)
2823 d_a_work(ind+j)=d_a_old(j,i)
2828 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2830 dc_work(ind+j)=dc_old(j,i+nres)
2831 d_t_work(ind+j)=d_t_old(j,i+nres)
2832 d_a_work(ind+j)=d_a_old(j,i+nres)
2840 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2844 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2845 & vfric_mat(i,j),afric_mat(i,j),
2846 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2854 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2855 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2856 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2857 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2859 d_t_work_new(i)=ddt1+0.5d0*ddt2
2860 d_t_work(i)=ddt1+ddt2
2865 d_t(j,0)=d_t_work(j)
2870 dc(j,i)=dc_work(ind+j)
2871 d_t(j,i)=d_t_work(ind+j)
2876 if (itype(i).ne.10) then
2879 dc(j,inres)=dc_work(ind+j)
2880 d_t(j,inres)=d_t_work(ind+j)
2887 c--------------------------------------------------------------------------
2888 subroutine sd_verlet2
2889 c Calculating the adjusted velocities for accelerations
2891 include 'DIMENSIONS'
2892 include 'COMMON.CONTROL'
2893 include 'COMMON.VAR'
2896 include 'COMMON.LAGRANGE.5diag'
2898 include 'COMMON.LAGRANGE'
2901 include 'COMMON.LANGEVIN'
2904 include 'COMMON.LANGEVIN.lang0.5diag'
2906 include 'COMMON.LANGEVIN.lang0'
2909 include 'COMMON.CHAIN'
2910 include 'COMMON.DERIV'
2911 include 'COMMON.GEO'
2912 include 'COMMON.LOCAL'
2913 include 'COMMON.INTERACT'
2914 include 'COMMON.IOUNITS'
2915 include 'COMMON.NAMES'
2916 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2917 common /stochcalc/ stochforcvec
2918 integer i,j,ind,inres
2920 c Compute the stochastic forces which contribute to velocity change
2922 call stochastic_force(stochforcvecV)
2929 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2930 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2931 & vrand_mat2(i,j)*stochforcvecV(j)
2933 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2937 d_t(j,0)=d_t_work(j)
2942 d_t(j,i)=d_t_work(ind+j)
2947 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2950 d_t(j,inres)=d_t_work(ind+j)
2957 c-----------------------------------------------------------
2958 subroutine sd_verlet_ciccotti_setup
2959 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2962 include 'DIMENSIONS'
2966 include 'COMMON.CONTROL'
2967 include 'COMMON.VAR'
2970 include 'COMMON.LAGRANGE.5diag'
2972 include 'COMMON.LAGRANGE'
2974 include 'COMMON.LANGEVIN'
2975 include 'COMMON.CHAIN'
2976 include 'COMMON.DERIV'
2977 include 'COMMON.GEO'
2978 include 'COMMON.LOCAL'
2979 include 'COMMON.INTERACT'
2980 include 'COMMON.IOUNITS'
2981 include 'COMMON.NAMES'
2982 include 'COMMON.TIME1'
2983 double precision emgdt(MAXRES6),
2984 & pterm,vterm,rho,rhoc,vsig,
2985 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2986 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2987 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2988 logical lprn /.false./
2989 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2990 double precision ktm
2998 c AL 8/17/04 Code adapted from tinker
3000 c Get the frictional and random terms for stochastic dynamics in the
3001 c eigenspace of mass-scaled UNRES friction matrix
3004 write (iout,*) "i",i," fricgam",fricgam(i)
3005 gdt = fricgam(i) * d_time
3007 c Stochastic dynamics reduces to simple MD for zero friction
3009 if (gdt .le. zero) then
3010 pfric_vec(i) = 1.0d0
3011 vfric_vec(i) = d_time
3012 afric_vec(i) = 0.5d0*d_time*d_time
3013 prand_vec(i) = afric_vec(i)
3014 vrand_vec2(i) = vfric_vec(i)
3016 c Analytical expressions when friction coefficient is large
3021 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3022 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3023 prand_vec(i) = afric_vec(i)
3024 vrand_vec2(i) = vfric_vec(i)
3026 c Compute the scaling factors of random terms for the nonzero friction case
3028 c ktm = 0.5d0*d_time/fricgam(i)
3029 c psig = dsqrt(ktm*pterm) / fricgam(i)
3030 c vsig = dsqrt(ktm*vterm)
3031 c prand_vec(i) = psig*afric_vec(i)
3032 c vrand_vec2(i) = vsig*vfric_vec(i)
3037 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
3040 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
3041 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3045 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3047 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3048 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3049 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3050 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3051 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3053 t_sdsetup=t_sdsetup+MPI_Wtime()
3055 t_sdsetup=t_sdsetup+tcpu()-tt0
3059 c-------------------------------------------------------------
3060 subroutine sd_verlet1_ciccotti
3061 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
3063 include 'DIMENSIONS'
3067 include 'COMMON.CONTROL'
3068 include 'COMMON.VAR'
3071 include 'COMMON.LAGRANGE.5diag'
3073 include 'COMMON.LAGRANGE'
3075 include 'COMMON.LANGEVIN'
3076 include 'COMMON.CHAIN'
3077 include 'COMMON.DERIV'
3078 include 'COMMON.GEO'
3079 include 'COMMON.LOCAL'
3080 include 'COMMON.INTERACT'
3081 include 'COMMON.IOUNITS'
3082 include 'COMMON.NAMES'
3083 double precision stochforcvec(MAXRES6)
3084 common /stochcalc/ stochforcvec
3085 logical lprn /.false./
3088 c write (iout,*) "dc_old"
3090 c write (iout,'(i5,3f10.5,5x,3f10.5)')
3091 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3094 dc_work(j)=dc_old(j,0)
3095 d_t_work(j)=d_t_old(j,0)
3096 d_a_work(j)=d_a_old(j,0)
3101 dc_work(ind+j)=dc_old(j,i)
3102 d_t_work(ind+j)=d_t_old(j,i)
3103 d_a_work(ind+j)=d_a_old(j,i)
3108 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3110 dc_work(ind+j)=dc_old(j,i+nres)
3111 d_t_work(ind+j)=d_t_old(j,i+nres)
3112 d_a_work(ind+j)=d_a_old(j,i+nres)
3121 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
3125 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
3126 & vfric_mat(i,j),afric_mat(i,j),
3127 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3135 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
3136 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3137 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3138 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3140 d_t_work_new(i)=ddt1+0.5d0*ddt2
3141 d_t_work(i)=ddt1+ddt2
3146 d_t(j,0)=d_t_work(j)
3151 dc(j,i)=dc_work(ind+j)
3152 d_t(j,i)=d_t_work(ind+j)
3157 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3160 dc(j,inres)=dc_work(ind+j)
3161 d_t(j,inres)=d_t_work(ind+j)
3168 c--------------------------------------------------------------------------
3169 subroutine sd_verlet2_ciccotti
3170 c Calculating the adjusted velocities for accelerations
3172 include 'DIMENSIONS'
3173 include 'COMMON.CONTROL'
3174 include 'COMMON.VAR'
3177 include 'COMMON.LAGRANGE.5diag'
3179 include 'COMMON.LAGRANGE'
3181 include 'COMMON.LANGEVIN'
3182 include 'COMMON.CHAIN'
3183 include 'COMMON.DERIV'
3184 include 'COMMON.GEO'
3185 include 'COMMON.LOCAL'
3186 include 'COMMON.INTERACT'
3187 include 'COMMON.IOUNITS'
3188 include 'COMMON.NAMES'
3189 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
3190 common /stochcalc/ stochforcvec
3193 c Compute the stochastic forces which contribute to velocity change
3195 call stochastic_force(stochforcvecV)
3202 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3203 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3204 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3206 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3210 d_t(j,0)=d_t_work(j)
3215 d_t(j,i)=d_t_work(ind+j)
3220 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3223 d_t(j,inres)=d_t_work(ind+j)