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 write (iout,*) "init_MD INDPDB",indpdb
1676 c write(iout,*) "d_time", d_time
1677 c Compute the standard deviations of stochastic forces for Langevin dynamics
1678 c if the friction coefficients do not depend on surface area
1679 if (lang.gt.0 .and. .not.surfarea) then
1681 stdforcp(i)=stdfp*dsqrt(gamp)
1684 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1685 & *dsqrt(gamsc(iabs(itype(i))))
1688 c Open the pdb file for snapshotshots
1691 if (ilen(tmpdir).gt.0)
1692 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1693 & liczba(:ilen(liczba))//".pdb")
1695 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1699 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1700 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1701 & liczba(:ilen(liczba))//".x")
1702 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1705 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1706 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1707 & liczba(:ilen(liczba))//".cx")
1708 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1714 if (ilen(tmpdir).gt.0)
1715 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1716 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1718 if (ilen(tmpdir).gt.0)
1719 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1720 cartname=prefix(:ilen(prefix))//"_MD.cx"
1724 write (qstr,'(256(1h ))')
1727 iq = qinfrag(i,iset)*10
1728 iw = wfrag(i,iset)/100
1730 if(me.eq.king.or..not.out1file)
1731 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1732 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1737 iq = qinpair(i,iset)*10
1738 iw = wpair(i,iset)/100
1740 if(me.eq.king.or..not.out1file)
1741 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1742 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1746 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1748 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1750 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1752 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1755 write (iout,*) "REST ",rest
1757 if (restart1file) then
1759 & inquire(file=mremd_rst_name,exist=file_exist)
1761 write (*,*) me," Before broadcast: file_exist",file_exist
1762 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1764 write (*,*) me," After broadcast: file_exist",file_exist
1765 c inquire(file=mremd_rst_name,exist=file_exist)
1767 if(me.eq.king.or..not.out1file)
1768 & write(iout,*) "Initial state read by master and distributed"
1770 if (ilen(tmpdir).gt.0)
1771 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1772 & //liczba(:ilen(liczba))//'.rst')
1773 inquire(file=rest2name,exist=file_exist)
1776 if(.not.restart1file) then
1777 if(me.eq.king.or..not.out1file)
1778 & write(iout,*) "Initial state will be read from file ",
1779 & rest2name(:ilen(rest2name))
1782 call rescale_weights(t_bath)
1785 if(me.eq.king.or..not.out1file)then
1786 if (restart1file) then
1787 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1790 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1793 write(iout,*) "Initial velocities randomly generated"
1800 c Generate initial velocities
1801 if(me.eq.king.or..not.out1file)
1802 & write(iout,*) "Initial velocities randomly generated"
1805 CtotTafm is the variable for AFM time which eclipsed during
1808 c rest2name = prefix(:ilen(prefix))//'.rst'
1809 if(me.eq.king.or..not.out1file)then
1810 write (iout,*) "Initial velocities"
1812 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1813 & (d_t(j,i+nres),j=1,3)
1817 c Zeroing the total angular momentum of the system
1818 write(iout,*) "Calling the zero-angular
1819 & momentum subroutine"
1821 c Getting the potential energy and forces and velocities and accelerations
1823 write (iout,*) "velocity of the center of the mass:"
1824 write (iout,*) (vcm(j),j=1,3)
1827 d_t(j,0)=d_t(j,0)-vcm(j)
1829 c Removing the velocity of the center of mass
1831 if(me.eq.king.or..not.out1file)then
1832 write (iout,*) "vcm right after adjustment:"
1833 write (iout,*) (vcm(j),j=1,3)
1834 write (iout,*) "Initial velocities after adjustment"
1836 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1837 & (d_t(j,i+nres),j=1,3)
1842 write (iout,*) "init_MD before initial structure REST ",rest
1845 if (iranconf.ne.0) then
1846 c 8/22/17 AL Loop to produce a low-energy random conformation
1850 print *, 'Calling OVERLAP_SC'
1851 call overlap_sc(fail)
1855 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1856 print *,'SC_move',nft_sc,etot
1857 if(me.eq.king.or..not.out1file)
1858 & write(iout,*) 'SC_move',nft_sc,etot
1862 if (me.eq.king.or..not.out1file) write(iout,*)
1863 & 'Minimizing random structure: Calling MINIM_DC'
1864 call minim_dc(etot,iretcode,nfun)
1866 call geom_to_var(nvar,varia)
1867 if (me.eq.king.or..not.out1file) write (iout,*)
1868 & 'Minimizing random structure: Calling MINIMIZE.'
1869 call minimize(etot,varia,iretcode,nfun)
1870 call var_to_geom(nvar,varia)
1872 if (me.eq.king.or..not.out1file)
1873 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1874 if (isnan(etot) .or. etot.gt.1.0d4) then
1875 write (iout,*) "Energy too large",etot,
1876 & " trying another random conformation"
1879 call gen_rand_conf(itmp,*30)
1881 30 write (iout,*) 'Failed to generate random conformation',
1882 & ', itrial=',itrial
1883 write (*,*) 'Processor:',me,
1884 & ' Failed to generate random conformation',
1894 write (iout,'(a,i3,a)') 'Processor:',me,
1895 & ' error in generating random conformation.'
1896 write (*,'(a,i3,a)') 'Processor:',me,
1897 & ' error in generating random conformation.'
1900 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1909 write (iout,'(a,i3,a)') 'Processor:',me,
1910 & ' failed to generate a low-energy random conformation.'
1911 write (*,'(a,i3,a)') 'Processor:',me,
1912 & ' failed to generate a low-energy random conformation.'
1915 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1920 else if (preminim) then
1921 if (start_from_model) then
1923 do while (fail .and. n_model_try.lt.constr_homology)
1925 i_model=iran_num(1,constr_homology)
1927 if (i_model.eq.list_model_try(k)) exit
1929 if (k.gt.n_model_try) exit
1931 n_model_try=n_model_try+1
1932 list_model_try(n_model_try)=i_model
1933 write (iout,*) 'starting from model ',i_model
1936 c(j,i)=chomo(j,i,i_model)
1939 call int_from_cart(.true.,.false.)
1940 call sc_loc_geom(.false.)
1943 dc(j,i)=c(j,i+1)-c(j,i)
1944 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1949 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1950 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1953 if (me.eq.king.or..not.out1file) then
1954 write (iout,*) "Energies before removing overlaps"
1955 call etotal(energia(0))
1956 call enerprint(energia(0))
1958 ! Remove SC overlaps if requested
1960 write (iout,*) 'Calling OVERLAP_SC'
1961 call overlap_sc(fail)
1964 & "Failed to remove overlap from model",i_model
1969 if (me.eq.king.or..not.out1file) then
1970 write (iout,*) "Energies after removing overlaps"
1971 call etotal(energia(0))
1972 call enerprint(energia(0))
1974 ! Search for better SC rotamers if requested
1976 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1977 print *,'SC_move',nft_sc,etot
1978 if (me.eq.king.or..not.out1file)
1979 & write(iout,*) 'SC_move',nft_sc,etot
1981 call etotal(energia(0))
1984 if (n_model_try.gt.constr_homology) then
1986 & "All models have irreparable overlaps. Trying randoms starts."
1991 C 8/22/17 AL Minimize initial structure
1993 if (me.eq.king.or..not.out1file) write(iout,*)
1994 & 'Minimizing initial PDB structure: Calling MINIM_DC'
1995 call minim_dc(etot,iretcode,nfun)
1997 call geom_to_var(nvar,varia)
1998 if(me.eq.king.or..not.out1file) write (iout,*)
1999 & 'Minimizing initial PDB structure: Calling MINIMIZE.'
2000 call minimize(etot,varia,iretcode,nfun)
2001 call var_to_geom(nvar,varia)
2003 if (me.eq.king.or..not.out1file)
2004 & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2005 if(me.eq.king.or..not.out1file)
2006 & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2008 if (me.eq.king.or..not.out1file)
2009 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2010 if(me.eq.king.or..not.out1file)
2011 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2016 call chainbuild_cart
2021 kinetic_T=2.0d0/(dimen3*Rb)*EK
2022 if(me.eq.king.or..not.out1file)then
2032 call etotal(potEcomp)
2033 if (large) call enerprint(potEcomp)
2036 t_etotal=t_etotal+MPI_Wtime()-tt0
2038 t_etotal=t_etotal+tcpu()-tt0
2042 c write (iout,*) "PotE-homology",potE-potEcomp(27)
2046 if (amax*d_time .gt. dvmax) then
2047 d_time=d_time*dvmax/amax
2048 if(me.eq.king.or..not.out1file) write (iout,*)
2049 & "Time step reduced to",d_time,
2050 & " because of too large initial acceleration."
2052 if(me.eq.king.or..not.out1file)then
2053 write(iout,*) "Potential energy and its components"
2054 call enerprint(potEcomp)
2055 c write(iout,*) (potEcomp(i),i=0,n_ene)
2057 potE=potEcomp(0)-potEcomp(27)
2058 c write (iout,*) "PotE-homology",potE
2062 if (ntwe.ne.0) call statout(itime)
2063 if(me.eq.king.or..not.out1file)
2064 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
2065 & " Kinetic energy",EK," Potential energy",potE,
2066 & " Total energy",totE," Maximum acceleration ",
2069 write (iout,*) "Initial coordinates"
2071 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
2072 & (c(j,i+nres),j=1,3)
2074 write (iout,*) "Initial dC"
2076 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
2077 & (dc(j,i+nres),j=1,3)
2079 write (iout,*) "Initial velocities"
2080 write (iout,"(13x,' backbone ',23x,' side chain')")
2082 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
2083 & (d_t(j,i+nres),j=1,3)
2085 write (iout,*) "Initial accelerations"
2087 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2088 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
2089 & (d_a(j,i+nres),j=1,3)
2095 d_t_old(j,i)=d_t(j,i)
2096 d_a_old(j,i)=d_a(j,i)
2098 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2107 call etotal_short(energia_short)
2108 if (large) call enerprint(potEcomp)
2111 t_eshort=t_eshort+MPI_Wtime()-tt0
2113 t_eshort=t_eshort+tcpu()-tt0
2118 if(.not.out1file .and. large) then
2119 write (iout,*) "energia_long",energia_long(0),
2120 & " energia_short",energia_short(0),
2121 & " total",energia_long(0)+energia_short(0)
2122 write (iout,*) "Initial fast-force accelerations"
2124 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2125 & (d_a(j,i+nres),j=1,3)
2128 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2131 d_a_short(j,i)=d_a(j,i)
2140 call etotal_long(energia_long)
2141 if (large) call enerprint(potEcomp)
2144 t_elong=t_elong+MPI_Wtime()-tt0
2146 t_elong=t_elong+tcpu()-tt0
2151 if(.not.out1file .and. large) then
2152 write (iout,*) "energia_long",energia_long(0)
2153 write (iout,*) "Initial slow-force accelerations"
2155 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2156 & (d_a(j,i+nres),j=1,3)
2160 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2162 t_enegrad=t_enegrad+tcpu()-tt0
2167 c-----------------------------------------------------------
2168 subroutine random_vel
2170 include 'DIMENSIONS'
2171 include 'COMMON.CONTROL'
2172 include 'COMMON.VAR'
2175 include 'COMMON.LAGRANGE.5diag'
2177 include 'COMMON.LAGRANGE'
2180 include 'COMMON.LANGEVIN'
2183 include 'COMMON.LANGEVIN.lang0.5diag'
2185 include 'COMMON.LANGEVIN.lang0'
2188 include 'COMMON.CHAIN'
2189 include 'COMMON.DERIV'
2190 include 'COMMON.GEO'
2191 include 'COMMON.LOCAL'
2192 include 'COMMON.INTERACT'
2193 include 'COMMON.IOUNITS'
2194 include 'COMMON.NAMES'
2195 include 'COMMON.TIME1'
2196 double precision xv,sigv,lowb,highb,vec_afm(3),Ek1,Ek2,Ek3,aux
2197 integer i,ii,j,k,l,ind
2198 double precision anorm_distr
2199 logical lprn /.false./
2201 integer ichain,n,innt,inct,ibeg,ierr
2202 double precision work(8*maxres6)
2203 integer iwork(maxres6)
2204 double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),
2205 & Gvec(maxres2_chain,maxres2_chain)
2206 common /przechowalnia/Ghalf,Geigen,Gvec
2208 double precision inertia(maxres2_chain,maxres2_chain)
2210 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2211 c First generate velocities in the eigenspace of the G matrix
2212 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2215 write (iout,*) "Random_vel, fivediag"
2224 n=dimen_chain(ichain)
2225 innt=iposd_chain(ichain)
2228 write (iout,*) "Chain",ichain," n",n," start",innt
2230 if (i.lt.inct-1) then
2231 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),
2233 else if (i.eq.inct-1) then
2234 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
2236 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
2240 ghalf(ind+1)=dmorig(innt)
2241 ghalf(ind+2)=du1orig(innt)
2242 ghalf(ind+3)=dmorig(innt+1)
2246 c write (iout,*) "i",i," ind",ind," indu2",innt+i-2,
2247 c & " indu1",innt+i-1," indm",innt+i
2248 ghalf(ind+1)=du2orig(innt-1+i-2)
2249 ghalf(ind+2)=du1orig(innt-1+i-1)
2250 ghalf(ind+3)=dmorig(innt-1+i)
2251 c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
2252 c & "DU1",innt-1+i-1,"DM ",innt-1+i
2260 inertia(i,j)=ghalf(ind)
2261 inertia(j,i)=ghalf(ind)
2266 write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
2267 write (iout,*) "Five-diagonal inertia matrix, lower triangle"
2268 call matoutr(n,ghalf)
2270 call gldiag(maxres2_chain,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
2272 write (iout,'(//a,i3)')
2273 & "Eigenvectors and eigenvalues of the G matrix chain",ichain
2274 call eigout(n,n,maxres2_chain,maxres2_chain,Gvec,Geigen)
2277 c check diagonalization
2283 aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
2287 write (iout,*) i,j,aux,geigen(i)
2289 write (iout,*) i,j,aux
2299 sigv=dsqrt((Rb*t_bath)/geigen(i))
2302 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2303 EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
2304 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2305 c & " d_t_work_new",d_t_work_new(ii)
2313 d_t_work(ind)=d_t_work(ind)
2314 & +Gvec(i,j)*d_t_work_new((j-1)*3+k)
2316 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2325 aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
2331 c Transfer to the d_t vector
2332 innt=chain_border(1,ichain)
2333 inct=chain_border(2,ichain)
2335 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
2339 d_t(j,i)=d_t_work(ind)
2341 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2344 d_t(j,i+nres)=d_t_work(ind)
2351 write (iout,*) "Random velocities in the Calpha,SC space"
2353 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2354 & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2357 call kinetic_CASC(Ek1)
2359 ! Transform the velocities to virtual-bond space
2367 if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
2369 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2373 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2374 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2387 c write (iout,*) "Shifting accelerations"
2389 c write (iout,*) "ichain",chain_border1(1,ichain)-1,
2390 c & chain_border1(1,ichain)
2391 d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
2392 d_t(:,chain_border1(1,ichain))=0.0d0
2394 c write (iout,*) "Adding accelerations"
2396 c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2397 c & chain_border(2,ichain-1)
2398 d_t(:,chain_border1(1,ichain)-1)=
2399 & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2400 d_t(:,chain_border(2,ichain-1))=0.0d0
2403 write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2404 & chain_border(2,ichain-1)
2405 d_t(:,chain_border1(1,ichain)-1)=
2406 & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2407 d_t(:,chain_border(2,ichain-1))=0.0d0
2412 c d_t(j,0)=d_t(j,nnt)
2415 innt=chain_border(1,ichain)
2416 inct=chain_border(2,ichain)
2417 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
2418 c write (iout,*) "ibeg",ibeg
2420 d_t(j,ibeg)=d_t(j,innt)
2424 if (iabs(itype(i).eq.10)) then
2425 c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
2427 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2431 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2432 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2441 & "Random velocities in the virtual-bond-vector space"
2442 write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
2444 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2445 & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2448 write (iout,*) "Kinetic energy from inertia matrix eigenvalues",
2451 & "Kinetic temperatures from inertia matrix eigenvalues",
2454 write (iout,*) "Kinetic energy from inertia matrix",Ek3
2455 write (iout,*) "Kinetic temperatures from inertia",
2456 & 2*Ek3/(3*dimen*Rb)
2458 write (iout,*) "Kinetic energy from velocities in CA-SC space",
2461 & "Kinetic temperatures from velovities in CA-SC space",
2462 & 2*Ek1/(3*dimen*Rb)
2465 & "Kinetic energy from virtual-bond-vector velocities",Ek1
2467 & "Kinetic temperature from virtual-bond-vector velocities ",
2476 sigv=dsqrt((Rb*t_bath)/geigen(i))
2479 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2481 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2482 c & " d_t_work_new",d_t_work_new(ii)
2485 C if (SELFGUIDE.gt.0) then
2488 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
2489 C distance=distance+vec_afm(j)**2
2491 C distance=dsqrt(distance)
2493 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
2494 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
2495 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
2496 C & d_t_work_new(j+(afmend-1)*3)
2507 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2510 c write (iout,*) "Ek from eigenvectors",Ek1
2512 c Transform velocities to UNRES coordinate space
2518 d_t_work(ind)=d_t_work(ind)
2519 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2521 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2525 c Transfer to the d_t vector
2527 d_t(j,0)=d_t_work(j)
2533 d_t(j,i)=d_t_work(ind)
2537 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2540 d_t(j,i+nres)=d_t_work(ind)
2545 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2546 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2552 c-----------------------------------------------------------
2553 subroutine sd_verlet_p_setup
2554 c Sets up the parameters of stochastic Verlet algorithm
2556 include 'DIMENSIONS'
2560 include 'COMMON.CONTROL'
2561 include 'COMMON.VAR'
2564 include 'COMMON.LANGEVIN'
2567 include 'COMMON.LANGEVIN.lang0.5diag'
2569 include 'COMMON.LANGEVIN.lang0'
2572 include 'COMMON.CHAIN'
2573 include 'COMMON.DERIV'
2574 include 'COMMON.GEO'
2575 include 'COMMON.LOCAL'
2576 include 'COMMON.INTERACT'
2577 include 'COMMON.IOUNITS'
2578 include 'COMMON.NAMES'
2579 include 'COMMON.TIME1'
2580 double precision emgdt(MAXRES6),
2581 & pterm,vterm,rho,rhoc,vsig,
2582 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2583 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2584 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2585 logical lprn /.false./
2586 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2587 double precision ktm
2594 c AL 8/17/04 Code adapted from tinker
2596 c Get the frictional and random terms for stochastic dynamics in the
2597 c eigenspace of mass-scaled UNRES friction matrix
2600 gdt = fricgam(i) * d_time
2602 c Stochastic dynamics reduces to simple MD for zero friction
2604 if (gdt .le. zero) then
2605 pfric_vec(i) = 1.0d0
2606 vfric_vec(i) = d_time
2607 afric_vec(i) = 0.5d0 * d_time * d_time
2608 prand_vec(i) = 0.0d0
2609 vrand_vec1(i) = 0.0d0
2610 vrand_vec2(i) = 0.0d0
2612 c Analytical expressions when friction coefficient is large
2615 if (gdt .ge. gdt_radius) then
2618 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2619 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2620 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2621 vterm = 1.0d0 - egdt**2
2622 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2624 c Use series expansions when friction coefficient is small
2635 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2636 & - gdt5/120.0d0 + gdt6/720.0d0
2637 & - gdt7/5040.0d0 + gdt8/40320.0d0
2638 & - gdt9/362880.0d0) / fricgam(i)**2
2639 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2640 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2641 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2642 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2643 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2644 & + 127.0d0*gdt9/90720.0d0
2645 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2646 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2647 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2648 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2649 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2650 & - 17.0d0*gdt2/1280.0d0
2651 & + 17.0d0*gdt3/6144.0d0
2652 & + 40967.0d0*gdt4/34406400.0d0
2653 & - 57203.0d0*gdt5/275251200.0d0
2654 & - 1429487.0d0*gdt6/13212057600.0d0)
2657 c Compute the scaling factors of random terms for the nonzero friction case
2659 ktm = 0.5d0*d_time/fricgam(i)
2660 psig = dsqrt(ktm*pterm) / fricgam(i)
2661 vsig = dsqrt(ktm*vterm)
2662 rhoc = dsqrt(1.0d0 - rho*rho)
2664 vrand_vec1(i) = vsig * rho
2665 vrand_vec2(i) = vsig * rhoc
2670 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2673 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2674 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2678 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2680 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2681 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2682 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2683 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2684 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2685 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2687 t_sdsetup=t_sdsetup+MPI_Wtime()
2689 t_sdsetup=t_sdsetup+tcpu()-tt0
2693 c-------------------------------------------------------------
2694 subroutine eigtransf1(n,ndim,ab,d,c)
2697 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2703 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2709 c-------------------------------------------------------------
2710 subroutine eigtransf(n,ndim,a,b,d,c)
2713 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2719 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2725 c-------------------------------------------------------------
2726 subroutine sd_verlet1
2727 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2729 include 'DIMENSIONS'
2730 include 'COMMON.CONTROL'
2731 include 'COMMON.VAR'
2734 include 'COMMON.LAGRANGE.5diag'
2736 include 'COMMON.LAGRANGE'
2739 include 'COMMON.LANGEVIN'
2742 include 'COMMON.LANGEVIN.lang0.5diag'
2744 include 'COMMON.LANGEVIN.lang0'
2747 include 'COMMON.CHAIN'
2748 include 'COMMON.DERIV'
2749 include 'COMMON.GEO'
2750 include 'COMMON.LOCAL'
2751 include 'COMMON.INTERACT'
2752 include 'COMMON.IOUNITS'
2753 include 'COMMON.NAMES'
2754 double precision stochforcvec(MAXRES6)
2755 common /stochcalc/ stochforcvec
2756 logical lprn /.false./
2757 integer i,j,ind,inres
2759 c write (iout,*) "dc_old"
2761 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2762 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2765 dc_work(j)=dc_old(j,0)
2766 d_t_work(j)=d_t_old(j,0)
2767 d_a_work(j)=d_a_old(j,0)
2772 dc_work(ind+j)=dc_old(j,i)
2773 d_t_work(ind+j)=d_t_old(j,i)
2774 d_a_work(ind+j)=d_a_old(j,i)
2779 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2781 dc_work(ind+j)=dc_old(j,i+nres)
2782 d_t_work(ind+j)=d_t_old(j,i+nres)
2783 d_a_work(ind+j)=d_a_old(j,i+nres)
2791 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2795 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2796 & vfric_mat(i,j),afric_mat(i,j),
2797 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2805 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2806 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2807 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2808 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2810 d_t_work_new(i)=ddt1+0.5d0*ddt2
2811 d_t_work(i)=ddt1+ddt2
2816 d_t(j,0)=d_t_work(j)
2821 dc(j,i)=dc_work(ind+j)
2822 d_t(j,i)=d_t_work(ind+j)
2827 if (itype(i).ne.10) then
2830 dc(j,inres)=dc_work(ind+j)
2831 d_t(j,inres)=d_t_work(ind+j)
2838 c--------------------------------------------------------------------------
2839 subroutine sd_verlet2
2840 c Calculating the adjusted velocities for accelerations
2842 include 'DIMENSIONS'
2843 include 'COMMON.CONTROL'
2844 include 'COMMON.VAR'
2847 include 'COMMON.LAGRANGE.5diag'
2849 include 'COMMON.LAGRANGE'
2852 include 'COMMON.LANGEVIN'
2855 include 'COMMON.LANGEVIN.lang0.5diag'
2857 include 'COMMON.LANGEVIN.lang0'
2860 include 'COMMON.CHAIN'
2861 include 'COMMON.DERIV'
2862 include 'COMMON.GEO'
2863 include 'COMMON.LOCAL'
2864 include 'COMMON.INTERACT'
2865 include 'COMMON.IOUNITS'
2866 include 'COMMON.NAMES'
2867 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2868 common /stochcalc/ stochforcvec
2869 integer i,j,ind,inres
2871 c Compute the stochastic forces which contribute to velocity change
2873 call stochastic_force(stochforcvecV)
2880 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2881 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2882 & vrand_mat2(i,j)*stochforcvecV(j)
2884 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2888 d_t(j,0)=d_t_work(j)
2893 d_t(j,i)=d_t_work(ind+j)
2898 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2901 d_t(j,inres)=d_t_work(ind+j)
2908 c-----------------------------------------------------------
2909 subroutine sd_verlet_ciccotti_setup
2910 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2913 include 'DIMENSIONS'
2917 include 'COMMON.CONTROL'
2918 include 'COMMON.VAR'
2921 include 'COMMON.LAGRANGE.5diag'
2923 include 'COMMON.LAGRANGE'
2925 include 'COMMON.LANGEVIN'
2926 include 'COMMON.CHAIN'
2927 include 'COMMON.DERIV'
2928 include 'COMMON.GEO'
2929 include 'COMMON.LOCAL'
2930 include 'COMMON.INTERACT'
2931 include 'COMMON.IOUNITS'
2932 include 'COMMON.NAMES'
2933 include 'COMMON.TIME1'
2934 double precision emgdt(MAXRES6),
2935 & pterm,vterm,rho,rhoc,vsig,
2936 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2937 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2938 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2939 logical lprn /.false./
2940 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2941 double precision ktm
2949 c AL 8/17/04 Code adapted from tinker
2951 c Get the frictional and random terms for stochastic dynamics in the
2952 c eigenspace of mass-scaled UNRES friction matrix
2955 write (iout,*) "i",i," fricgam",fricgam(i)
2956 gdt = fricgam(i) * d_time
2958 c Stochastic dynamics reduces to simple MD for zero friction
2960 if (gdt .le. zero) then
2961 pfric_vec(i) = 1.0d0
2962 vfric_vec(i) = d_time
2963 afric_vec(i) = 0.5d0*d_time*d_time
2964 prand_vec(i) = afric_vec(i)
2965 vrand_vec2(i) = vfric_vec(i)
2967 c Analytical expressions when friction coefficient is large
2972 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2973 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2974 prand_vec(i) = afric_vec(i)
2975 vrand_vec2(i) = vfric_vec(i)
2977 c Compute the scaling factors of random terms for the nonzero friction case
2979 c ktm = 0.5d0*d_time/fricgam(i)
2980 c psig = dsqrt(ktm*pterm) / fricgam(i)
2981 c vsig = dsqrt(ktm*vterm)
2982 c prand_vec(i) = psig*afric_vec(i)
2983 c vrand_vec2(i) = vsig*vfric_vec(i)
2988 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2991 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2992 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2996 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2998 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2999 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3000 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3001 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3002 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3004 t_sdsetup=t_sdsetup+MPI_Wtime()
3006 t_sdsetup=t_sdsetup+tcpu()-tt0
3010 c-------------------------------------------------------------
3011 subroutine sd_verlet1_ciccotti
3012 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
3014 include 'DIMENSIONS'
3018 include 'COMMON.CONTROL'
3019 include 'COMMON.VAR'
3022 include 'COMMON.LAGRANGE.5diag'
3024 include 'COMMON.LAGRANGE'
3026 include 'COMMON.LANGEVIN'
3027 include 'COMMON.CHAIN'
3028 include 'COMMON.DERIV'
3029 include 'COMMON.GEO'
3030 include 'COMMON.LOCAL'
3031 include 'COMMON.INTERACT'
3032 include 'COMMON.IOUNITS'
3033 include 'COMMON.NAMES'
3034 double precision stochforcvec(MAXRES6)
3035 common /stochcalc/ stochforcvec
3036 logical lprn /.false./
3039 c write (iout,*) "dc_old"
3041 c write (iout,'(i5,3f10.5,5x,3f10.5)')
3042 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3045 dc_work(j)=dc_old(j,0)
3046 d_t_work(j)=d_t_old(j,0)
3047 d_a_work(j)=d_a_old(j,0)
3052 dc_work(ind+j)=dc_old(j,i)
3053 d_t_work(ind+j)=d_t_old(j,i)
3054 d_a_work(ind+j)=d_a_old(j,i)
3059 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3061 dc_work(ind+j)=dc_old(j,i+nres)
3062 d_t_work(ind+j)=d_t_old(j,i+nres)
3063 d_a_work(ind+j)=d_a_old(j,i+nres)
3072 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
3076 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
3077 & vfric_mat(i,j),afric_mat(i,j),
3078 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3086 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
3087 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3088 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3089 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3091 d_t_work_new(i)=ddt1+0.5d0*ddt2
3092 d_t_work(i)=ddt1+ddt2
3097 d_t(j,0)=d_t_work(j)
3102 dc(j,i)=dc_work(ind+j)
3103 d_t(j,i)=d_t_work(ind+j)
3108 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3111 dc(j,inres)=dc_work(ind+j)
3112 d_t(j,inres)=d_t_work(ind+j)
3119 c--------------------------------------------------------------------------
3120 subroutine sd_verlet2_ciccotti
3121 c Calculating the adjusted velocities for accelerations
3123 include 'DIMENSIONS'
3124 include 'COMMON.CONTROL'
3125 include 'COMMON.VAR'
3128 include 'COMMON.LAGRANGE.5diag'
3130 include 'COMMON.LAGRANGE'
3132 include 'COMMON.LANGEVIN'
3133 include 'COMMON.CHAIN'
3134 include 'COMMON.DERIV'
3135 include 'COMMON.GEO'
3136 include 'COMMON.LOCAL'
3137 include 'COMMON.INTERACT'
3138 include 'COMMON.IOUNITS'
3139 include 'COMMON.NAMES'
3140 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
3141 common /stochcalc/ stochforcvec
3144 c Compute the stochastic forces which contribute to velocity change
3146 call stochastic_force(stochforcvecV)
3153 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3154 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3155 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3157 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3161 d_t(j,0)=d_t_work(j)
3166 d_t(j,i)=d_t_work(ind+j)
3171 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3174 d_t(j,inres)=d_t_work(ind+j)