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.)
1960 dc(j,i)=c(j,i+1)-c(j,i)
1961 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1966 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1967 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1970 if (me.eq.king.or..not.out1file) then
1971 write (iout,*) "Energies before removing overlaps"
1972 call etotal(energia(0))
1973 call enerprint(energia(0))
1975 ! Remove SC overlaps if requested
1977 write (iout,*) 'Calling OVERLAP_SC'
1978 call overlap_sc(fail)
1981 & "Failed to remove overlap from model",i_model
1985 if (me.eq.king.or..not.out1file) then
1986 write (iout,*) "Energies after removing overlaps"
1987 call etotal(energia(0))
1988 call enerprint(energia(0))
1991 ! Search for better SC rotamers if requested
1993 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1994 print *,'SC_move',nft_sc,etot
1995 if (me.eq.king.or..not.out1file)
1996 & write(iout,*) 'SC_move',nft_sc,etot
1998 call etotal(energia(0))
2001 call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),
2002 & 1,MPI_INTEGER,king,CG_COMM,IERROR)
2003 if (n_model_try.gt.nmodel_start .and.
2004 & (me.eq.king .or. out1file)) then
2006 & "All models have irreparable overlaps. Trying randoms starts."
2008 i_model=nmodel_start+1
2012 ! Remove SC overlaps if requested
2014 write (iout,*) 'Calling OVERLAP_SC'
2015 call overlap_sc(fail)
2018 & "Failed to remove overlap"
2021 if (me.eq.king.or..not.out1file) then
2022 write (iout,*) "Energies after removing overlaps"
2023 call etotal(energia(0))
2024 call enerprint(energia(0))
2027 C 8/22/17 AL Minimize initial structure
2029 if (me.eq.king.or..not.out1file) write(iout,*)
2030 & 'Minimizing initial PDB structure: Calling MINIM_DC'
2031 call minim_dc(etot,iretcode,nfun)
2033 call geom_to_var(nvar,varia)
2034 if(me.eq.king.or..not.out1file) write (iout,*)
2035 & 'Minimizing initial PDB structure: Calling MINIMIZE.'
2036 call minimize(etot,varia,iretcode,nfun)
2037 call var_to_geom(nvar,varia)
2039 if (me.eq.king.or..not.out1file)
2040 & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2041 if(me.eq.king.or..not.out1file)
2042 & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2044 if (me.eq.king.or..not.out1file)
2045 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2046 if(me.eq.king.or..not.out1file)
2047 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2051 if (nmodel_start.gt.0 .and. me.eq.king) then
2052 write (iout,'(a)') "Task Starting model"
2054 if (i_start_models(i).gt.nmodel_start) then
2055 write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE"
2057 write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i))
2058 & (:ilen(pdbfiles_chomo(i_start_models(i))))
2063 call chainbuild_cart
2068 kinetic_T=2.0d0/(dimen3*Rb)*EK
2069 write (iout,*) "Initial kinetic energy",EK," kinetic T",kinetic_T
2070 if(me.eq.king.or..not.out1file)then
2080 call etotal(potEcomp)
2081 if (large) call enerprint(potEcomp)
2084 t_etotal=t_etotal+MPI_Wtime()-tt0
2086 t_etotal=t_etotal+tcpu()-tt0
2090 c write (iout,*) "PotE-homology",potE-potEcomp(27)
2094 if (amax*d_time .gt. dvmax) then
2095 d_time=d_time*dvmax/amax
2096 if(me.eq.king.or..not.out1file) write (iout,*)
2097 & "Time step reduced to",d_time,
2098 & " because of too large initial acceleration."
2100 if(me.eq.king.or..not.out1file)then
2101 write(iout,*) "Potential energy and its components"
2102 call enerprint(potEcomp)
2103 c write(iout,*) (potEcomp(i),i=0,n_ene)
2105 potE=potEcomp(0)-potEcomp(27)
2106 c write (iout,*) "PotE-homology",potE
2110 if (ntwe.ne.0) call statout(itime)
2111 if(me.eq.king.or..not.out1file)
2112 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
2113 & " Kinetic energy",EK," Potential energy",potE,
2114 & " Total energy",totE," Maximum acceleration ",
2117 write (iout,*) "Initial coordinates"
2119 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
2120 & (c(j,i+nres),j=1,3)
2122 write (iout,*) "Initial dC"
2124 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
2125 & (dc(j,i+nres),j=1,3)
2127 write (iout,*) "Initial velocities"
2128 write (iout,"(13x,' backbone ',23x,' side chain')")
2130 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
2131 & (d_t(j,i+nres),j=1,3)
2133 write (iout,*) "Initial accelerations"
2135 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2136 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
2137 & (d_a(j,i+nres),j=1,3)
2143 d_t_old(j,i)=d_t(j,i)
2144 d_a_old(j,i)=d_a(j,i)
2146 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2155 call etotal_short(energia_short)
2156 if (large) call enerprint(potEcomp)
2159 t_eshort=t_eshort+MPI_Wtime()-tt0
2161 t_eshort=t_eshort+tcpu()-tt0
2166 if(.not.out1file .and. large) then
2167 write (iout,*) "energia_long",energia_long(0),
2168 & " energia_short",energia_short(0),
2169 & " total",energia_long(0)+energia_short(0)
2170 write (iout,*) "Initial fast-force accelerations"
2172 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2173 & (d_a(j,i+nres),j=1,3)
2176 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2179 d_a_short(j,i)=d_a(j,i)
2188 call etotal_long(energia_long)
2189 if (large) call enerprint(potEcomp)
2192 t_elong=t_elong+MPI_Wtime()-tt0
2194 t_elong=t_elong+tcpu()-tt0
2199 if(.not.out1file .and. large) then
2200 write (iout,*) "energia_long",energia_long(0)
2201 write (iout,*) "Initial slow-force accelerations"
2203 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2204 & (d_a(j,i+nres),j=1,3)
2208 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2210 t_enegrad=t_enegrad+tcpu()-tt0
2215 c-----------------------------------------------------------
2216 subroutine random_vel
2218 include 'DIMENSIONS'
2219 include 'COMMON.CONTROL'
2220 include 'COMMON.VAR'
2223 include 'COMMON.LAGRANGE.5diag'
2225 include 'COMMON.LAGRANGE'
2228 include 'COMMON.LANGEVIN'
2231 include 'COMMON.LANGEVIN.lang0.5diag'
2233 include 'COMMON.LANGEVIN.lang0'
2236 include 'COMMON.CHAIN'
2237 include 'COMMON.DERIV'
2238 include 'COMMON.GEO'
2239 include 'COMMON.LOCAL'
2240 include 'COMMON.INTERACT'
2241 include 'COMMON.IOUNITS'
2242 include 'COMMON.NAMES'
2243 include 'COMMON.TIME1'
2244 double precision xv,sigv,lowb,highb,vec_afm(3),Ek1,Ek2,Ek3,aux
2245 integer i,ii,j,k,l,ind
2246 double precision anorm_distr
2247 logical lprn /.false./
2249 integer ichain,n,innt,inct,ibeg,ierr
2250 double precision work(8*maxres6)
2251 integer iwork(maxres6)
2252 double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),
2253 & Gvec(maxres2_chain,maxres2_chain)
2254 common /przechowalnia/Ghalf,Geigen,Gvec
2256 double precision inertia(maxres2_chain,maxres2_chain)
2258 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2259 c First generate velocities in the eigenspace of the G matrix
2260 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2263 write (iout,*) "Random_vel, fivediag"
2272 n=dimen_chain(ichain)
2273 innt=iposd_chain(ichain)
2276 write (iout,*) "Chain",ichain," n",n," start",innt
2278 if (i.lt.inct-1) then
2279 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),
2281 else if (i.eq.inct-1) then
2282 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
2284 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
2288 ghalf(ind+1)=dmorig(innt)
2289 ghalf(ind+2)=du1orig(innt)
2290 ghalf(ind+3)=dmorig(innt+1)
2294 c write (iout,*) "i",i," ind",ind," indu2",innt+i-2,
2295 c & " indu1",innt+i-1," indm",innt+i
2296 ghalf(ind+1)=du2orig(innt-1+i-2)
2297 ghalf(ind+2)=du1orig(innt-1+i-1)
2298 ghalf(ind+3)=dmorig(innt-1+i)
2299 c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
2300 c & "DU1",innt-1+i-1,"DM ",innt-1+i
2308 inertia(i,j)=ghalf(ind)
2309 inertia(j,i)=ghalf(ind)
2314 write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
2315 write (iout,*) "Five-diagonal inertia matrix, lower triangle"
2316 call matoutr(n,ghalf)
2318 call gldiag(maxres2_chain,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
2320 write (iout,'(//a,i3)')
2321 & "Eigenvectors and eigenvalues of the G matrix chain",ichain
2322 call eigout(n,n,maxres2_chain,maxres2_chain,Gvec,Geigen)
2325 c check diagonalization
2331 aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
2335 write (iout,*) i,j,aux,geigen(i)
2337 write (iout,*) i,j,aux
2347 sigv=dsqrt((Rb*t_bath)/geigen(i))
2350 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2351 EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
2352 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2353 c & " d_t_work_new",d_t_work_new(ii)
2361 d_t_work(ind)=d_t_work(ind)
2362 & +Gvec(i,j)*d_t_work_new((j-1)*3+k)
2364 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2373 aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
2379 c Transfer to the d_t vector
2380 innt=chain_border(1,ichain)
2381 inct=chain_border(2,ichain)
2383 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
2387 d_t(j,i)=d_t_work(ind)
2389 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2392 d_t(j,i+nres)=d_t_work(ind)
2399 write (iout,*) "Random velocities in the Calpha,SC space"
2401 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2402 & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2405 call kinetic_CASC(Ek1)
2407 ! Transform the velocities to virtual-bond space
2415 if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
2417 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2421 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2422 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2435 c write (iout,*) "Shifting accelerations"
2437 c write (iout,*) "ichain",chain_border1(1,ichain)-1,
2438 c & chain_border1(1,ichain)
2439 d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
2440 d_t(:,chain_border1(1,ichain))=0.0d0
2442 c write (iout,*) "Adding accelerations"
2444 c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2445 c & chain_border(2,ichain-1)
2446 d_t(:,chain_border1(1,ichain)-1)=
2447 & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2448 d_t(:,chain_border(2,ichain-1))=0.0d0
2451 write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2452 & chain_border(2,ichain-1)
2453 d_t(:,chain_border1(1,ichain)-1)=
2454 & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2455 d_t(:,chain_border(2,ichain-1))=0.0d0
2460 c d_t(j,0)=d_t(j,nnt)
2463 innt=chain_border(1,ichain)
2464 inct=chain_border(2,ichain)
2465 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
2466 c write (iout,*) "ibeg",ibeg
2468 d_t(j,ibeg)=d_t(j,innt)
2472 if (iabs(itype(i).eq.10)) then
2473 c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
2475 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2479 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2480 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2489 & "Random velocities in the virtual-bond-vector space"
2490 write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
2492 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2493 & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2496 write (iout,*) "Kinetic energy from inertia matrix eigenvalues",
2499 & "Kinetic temperatures from inertia matrix eigenvalues",
2502 write (iout,*) "Kinetic energy from inertia matrix",Ek3
2503 write (iout,*) "Kinetic temperatures from inertia",
2504 & 2*Ek3/(3*dimen*Rb)
2506 write (iout,*) "Kinetic energy from velocities in CA-SC space",
2509 & "Kinetic temperatures from velovities in CA-SC space",
2510 & 2*Ek1/(3*dimen*Rb)
2513 & "Kinetic energy from virtual-bond-vector velocities",Ek1
2515 & "Kinetic temperature from virtual-bond-vector velocities ",
2524 sigv=dsqrt((Rb*t_bath)/geigen(i))
2527 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2529 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2530 c & " d_t_work_new",d_t_work_new(ii)
2533 C if (SELFGUIDE.gt.0) then
2536 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
2537 C distance=distance+vec_afm(j)**2
2539 C distance=dsqrt(distance)
2541 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
2542 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
2543 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
2544 C & d_t_work_new(j+(afmend-1)*3)
2555 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2558 c write (iout,*) "Ek from eigenvectors",Ek1
2560 c Transform velocities to UNRES coordinate space
2566 d_t_work(ind)=d_t_work(ind)
2567 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2569 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2573 c Transfer to the d_t vector
2575 d_t(j,0)=d_t_work(j)
2581 d_t(j,i)=d_t_work(ind)
2585 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2588 d_t(j,i+nres)=d_t_work(ind)
2593 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2594 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2600 c-----------------------------------------------------------
2601 subroutine sd_verlet_p_setup
2602 c Sets up the parameters of stochastic Verlet algorithm
2604 include 'DIMENSIONS'
2608 include 'COMMON.CONTROL'
2609 include 'COMMON.VAR'
2612 include 'COMMON.LANGEVIN'
2615 include 'COMMON.LANGEVIN.lang0.5diag'
2617 include 'COMMON.LANGEVIN.lang0'
2620 include 'COMMON.CHAIN'
2621 include 'COMMON.DERIV'
2622 include 'COMMON.GEO'
2623 include 'COMMON.LOCAL'
2624 include 'COMMON.INTERACT'
2625 include 'COMMON.IOUNITS'
2626 include 'COMMON.NAMES'
2627 include 'COMMON.TIME1'
2628 double precision emgdt(MAXRES6),
2629 & pterm,vterm,rho,rhoc,vsig,
2630 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2631 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2632 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2633 logical lprn /.false./
2634 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2635 double precision ktm
2642 c AL 8/17/04 Code adapted from tinker
2644 c Get the frictional and random terms for stochastic dynamics in the
2645 c eigenspace of mass-scaled UNRES friction matrix
2648 gdt = fricgam(i) * d_time
2650 c Stochastic dynamics reduces to simple MD for zero friction
2652 if (gdt .le. zero) then
2653 pfric_vec(i) = 1.0d0
2654 vfric_vec(i) = d_time
2655 afric_vec(i) = 0.5d0 * d_time * d_time
2656 prand_vec(i) = 0.0d0
2657 vrand_vec1(i) = 0.0d0
2658 vrand_vec2(i) = 0.0d0
2660 c Analytical expressions when friction coefficient is large
2663 if (gdt .ge. gdt_radius) then
2666 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2667 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2668 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2669 vterm = 1.0d0 - egdt**2
2670 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2672 c Use series expansions when friction coefficient is small
2683 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2684 & - gdt5/120.0d0 + gdt6/720.0d0
2685 & - gdt7/5040.0d0 + gdt8/40320.0d0
2686 & - gdt9/362880.0d0) / fricgam(i)**2
2687 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2688 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2689 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2690 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2691 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2692 & + 127.0d0*gdt9/90720.0d0
2693 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2694 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2695 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2696 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2697 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2698 & - 17.0d0*gdt2/1280.0d0
2699 & + 17.0d0*gdt3/6144.0d0
2700 & + 40967.0d0*gdt4/34406400.0d0
2701 & - 57203.0d0*gdt5/275251200.0d0
2702 & - 1429487.0d0*gdt6/13212057600.0d0)
2705 c Compute the scaling factors of random terms for the nonzero friction case
2707 ktm = 0.5d0*d_time/fricgam(i)
2708 psig = dsqrt(ktm*pterm) / fricgam(i)
2709 vsig = dsqrt(ktm*vterm)
2710 rhoc = dsqrt(1.0d0 - rho*rho)
2712 vrand_vec1(i) = vsig * rho
2713 vrand_vec2(i) = vsig * rhoc
2718 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2721 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2722 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2726 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2728 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2729 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2730 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2731 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2732 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2733 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2735 t_sdsetup=t_sdsetup+MPI_Wtime()
2737 t_sdsetup=t_sdsetup+tcpu()-tt0
2741 c-------------------------------------------------------------
2742 subroutine eigtransf1(n,ndim,ab,d,c)
2745 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2751 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2757 c-------------------------------------------------------------
2758 subroutine eigtransf(n,ndim,a,b,d,c)
2761 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2767 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2773 c-------------------------------------------------------------
2774 subroutine sd_verlet1
2775 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2777 include 'DIMENSIONS'
2778 include 'COMMON.CONTROL'
2779 include 'COMMON.VAR'
2782 include 'COMMON.LAGRANGE.5diag'
2784 include 'COMMON.LAGRANGE'
2787 include 'COMMON.LANGEVIN'
2790 include 'COMMON.LANGEVIN.lang0.5diag'
2792 include 'COMMON.LANGEVIN.lang0'
2795 include 'COMMON.CHAIN'
2796 include 'COMMON.DERIV'
2797 include 'COMMON.GEO'
2798 include 'COMMON.LOCAL'
2799 include 'COMMON.INTERACT'
2800 include 'COMMON.IOUNITS'
2801 include 'COMMON.NAMES'
2802 double precision stochforcvec(MAXRES6)
2803 common /stochcalc/ stochforcvec
2804 logical lprn /.false./
2805 integer i,j,ind,inres
2807 c write (iout,*) "dc_old"
2809 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2810 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2813 dc_work(j)=dc_old(j,0)
2814 d_t_work(j)=d_t_old(j,0)
2815 d_a_work(j)=d_a_old(j,0)
2820 dc_work(ind+j)=dc_old(j,i)
2821 d_t_work(ind+j)=d_t_old(j,i)
2822 d_a_work(ind+j)=d_a_old(j,i)
2827 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2829 dc_work(ind+j)=dc_old(j,i+nres)
2830 d_t_work(ind+j)=d_t_old(j,i+nres)
2831 d_a_work(ind+j)=d_a_old(j,i+nres)
2839 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2843 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2844 & vfric_mat(i,j),afric_mat(i,j),
2845 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2853 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2854 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2855 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2856 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2858 d_t_work_new(i)=ddt1+0.5d0*ddt2
2859 d_t_work(i)=ddt1+ddt2
2864 d_t(j,0)=d_t_work(j)
2869 dc(j,i)=dc_work(ind+j)
2870 d_t(j,i)=d_t_work(ind+j)
2875 if (itype(i).ne.10) then
2878 dc(j,inres)=dc_work(ind+j)
2879 d_t(j,inres)=d_t_work(ind+j)
2886 c--------------------------------------------------------------------------
2887 subroutine sd_verlet2
2888 c Calculating the adjusted velocities for accelerations
2890 include 'DIMENSIONS'
2891 include 'COMMON.CONTROL'
2892 include 'COMMON.VAR'
2895 include 'COMMON.LAGRANGE.5diag'
2897 include 'COMMON.LAGRANGE'
2900 include 'COMMON.LANGEVIN'
2903 include 'COMMON.LANGEVIN.lang0.5diag'
2905 include 'COMMON.LANGEVIN.lang0'
2908 include 'COMMON.CHAIN'
2909 include 'COMMON.DERIV'
2910 include 'COMMON.GEO'
2911 include 'COMMON.LOCAL'
2912 include 'COMMON.INTERACT'
2913 include 'COMMON.IOUNITS'
2914 include 'COMMON.NAMES'
2915 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2916 common /stochcalc/ stochforcvec
2917 integer i,j,ind,inres
2919 c Compute the stochastic forces which contribute to velocity change
2921 call stochastic_force(stochforcvecV)
2928 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2929 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2930 & vrand_mat2(i,j)*stochforcvecV(j)
2932 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2936 d_t(j,0)=d_t_work(j)
2941 d_t(j,i)=d_t_work(ind+j)
2946 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2949 d_t(j,inres)=d_t_work(ind+j)
2956 c-----------------------------------------------------------
2957 subroutine sd_verlet_ciccotti_setup
2958 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2961 include 'DIMENSIONS'
2965 include 'COMMON.CONTROL'
2966 include 'COMMON.VAR'
2969 include 'COMMON.LAGRANGE.5diag'
2971 include 'COMMON.LAGRANGE'
2973 include 'COMMON.LANGEVIN'
2974 include 'COMMON.CHAIN'
2975 include 'COMMON.DERIV'
2976 include 'COMMON.GEO'
2977 include 'COMMON.LOCAL'
2978 include 'COMMON.INTERACT'
2979 include 'COMMON.IOUNITS'
2980 include 'COMMON.NAMES'
2981 include 'COMMON.TIME1'
2982 double precision emgdt(MAXRES6),
2983 & pterm,vterm,rho,rhoc,vsig,
2984 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2985 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2986 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2987 logical lprn /.false./
2988 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2989 double precision ktm
2997 c AL 8/17/04 Code adapted from tinker
2999 c Get the frictional and random terms for stochastic dynamics in the
3000 c eigenspace of mass-scaled UNRES friction matrix
3003 write (iout,*) "i",i," fricgam",fricgam(i)
3004 gdt = fricgam(i) * d_time
3006 c Stochastic dynamics reduces to simple MD for zero friction
3008 if (gdt .le. zero) then
3009 pfric_vec(i) = 1.0d0
3010 vfric_vec(i) = d_time
3011 afric_vec(i) = 0.5d0*d_time*d_time
3012 prand_vec(i) = afric_vec(i)
3013 vrand_vec2(i) = vfric_vec(i)
3015 c Analytical expressions when friction coefficient is large
3020 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3021 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3022 prand_vec(i) = afric_vec(i)
3023 vrand_vec2(i) = vfric_vec(i)
3025 c Compute the scaling factors of random terms for the nonzero friction case
3027 c ktm = 0.5d0*d_time/fricgam(i)
3028 c psig = dsqrt(ktm*pterm) / fricgam(i)
3029 c vsig = dsqrt(ktm*vterm)
3030 c prand_vec(i) = psig*afric_vec(i)
3031 c vrand_vec2(i) = vsig*vfric_vec(i)
3036 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
3039 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
3040 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3044 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3046 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3047 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3048 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3049 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3050 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3052 t_sdsetup=t_sdsetup+MPI_Wtime()
3054 t_sdsetup=t_sdsetup+tcpu()-tt0
3058 c-------------------------------------------------------------
3059 subroutine sd_verlet1_ciccotti
3060 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
3062 include 'DIMENSIONS'
3066 include 'COMMON.CONTROL'
3067 include 'COMMON.VAR'
3070 include 'COMMON.LAGRANGE.5diag'
3072 include 'COMMON.LAGRANGE'
3074 include 'COMMON.LANGEVIN'
3075 include 'COMMON.CHAIN'
3076 include 'COMMON.DERIV'
3077 include 'COMMON.GEO'
3078 include 'COMMON.LOCAL'
3079 include 'COMMON.INTERACT'
3080 include 'COMMON.IOUNITS'
3081 include 'COMMON.NAMES'
3082 double precision stochforcvec(MAXRES6)
3083 common /stochcalc/ stochforcvec
3084 logical lprn /.false./
3087 c write (iout,*) "dc_old"
3089 c write (iout,'(i5,3f10.5,5x,3f10.5)')
3090 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3093 dc_work(j)=dc_old(j,0)
3094 d_t_work(j)=d_t_old(j,0)
3095 d_a_work(j)=d_a_old(j,0)
3100 dc_work(ind+j)=dc_old(j,i)
3101 d_t_work(ind+j)=d_t_old(j,i)
3102 d_a_work(ind+j)=d_a_old(j,i)
3107 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3109 dc_work(ind+j)=dc_old(j,i+nres)
3110 d_t_work(ind+j)=d_t_old(j,i+nres)
3111 d_a_work(ind+j)=d_a_old(j,i+nres)
3120 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
3124 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
3125 & vfric_mat(i,j),afric_mat(i,j),
3126 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3134 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
3135 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3136 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3137 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3139 d_t_work_new(i)=ddt1+0.5d0*ddt2
3140 d_t_work(i)=ddt1+ddt2
3145 d_t(j,0)=d_t_work(j)
3150 dc(j,i)=dc_work(ind+j)
3151 d_t(j,i)=d_t_work(ind+j)
3156 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3159 dc(j,inres)=dc_work(ind+j)
3160 d_t(j,inres)=d_t_work(ind+j)
3167 c--------------------------------------------------------------------------
3168 subroutine sd_verlet2_ciccotti
3169 c Calculating the adjusted velocities for accelerations
3171 include 'DIMENSIONS'
3172 include 'COMMON.CONTROL'
3173 include 'COMMON.VAR'
3176 include 'COMMON.LAGRANGE.5diag'
3178 include 'COMMON.LAGRANGE'
3180 include 'COMMON.LANGEVIN'
3181 include 'COMMON.CHAIN'
3182 include 'COMMON.DERIV'
3183 include 'COMMON.GEO'
3184 include 'COMMON.LOCAL'
3185 include 'COMMON.INTERACT'
3186 include 'COMMON.IOUNITS'
3187 include 'COMMON.NAMES'
3188 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
3189 common /stochcalc/ stochforcvec
3192 c Compute the stochastic forces which contribute to velocity change
3194 call stochastic_force(stochforcvecV)
3201 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3202 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3203 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3205 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3209 d_t(j,0)=d_t_work(j)
3214 d_t(j,i)=d_t_work(ind+j)
3219 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3222 d_t(j,inres)=d_t_work(ind+j)