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)
213 if (mod(itime,ntwe).eq.0) then
215 C call enerprint(potEcomp)
216 C print *,itime,'AFM',Eafmforc,etot
230 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
233 v_work(ind)=d_t(j,i+nres)
238 write (66,'(80f10.5)')
239 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
243 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
245 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
247 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
250 if (mod(itime,ntwx).eq.0) then
251 c write(iout,*) 'time=',itime
252 C call check_ecartint
254 write (tytul,'("time",f8.2)') totT
256 call hairpin(.true.,nharp,iharp)
257 call secondary2(.true.)
258 call pdbout(potE,tytul,ipdb)
263 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
264 open(irest2,file=rest2name,status='unknown')
265 write(irest2,*) totT,EK,potE,totE,t_bath
267 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
270 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
281 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
283 & 'MD calculations setup:',t_MDsetup,
284 & 'Energy & gradient evaluation:',t_enegrad,
285 & 'Stochastic MD setup:',t_langsetup,
286 & 'Stochastic MD step setup:',t_sdsetup,
288 write (iout,'(/28(1h=),a25,27(1h=))')
289 & ' End of MD calculation '
291 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
293 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
294 & " time_fricmatmult",time_fricmatmult," time_fsample ",
299 c-------------------------------------------------------------------------------
300 subroutine velverlet_step(itime)
301 c-------------------------------------------------------------------------------
302 c Perform a single velocity Verlet step; the time step can be rescaled if
303 c increments in accelerations exceed the threshold
304 c-------------------------------------------------------------------------------
309 integer ierror,ierrcode,errcode
311 include 'COMMON.SETUP'
312 include 'COMMON.CONTROL'
316 include 'COMMON.LAGRANGE.5diag'
318 include 'COMMON.LAGRANGE'
321 include 'COMMON.LANGEVIN'
324 include 'COMMON.LANGEVIN.lang0.5diag'
326 include 'COMMON.LANGEVIN.lang0'
329 include 'COMMON.CHAIN'
330 include 'COMMON.DERIV'
332 include 'COMMON.LOCAL'
333 include 'COMMON.INTERACT'
334 include 'COMMON.IOUNITS'
335 include 'COMMON.NAMES'
336 include 'COMMON.TIME1'
337 include 'COMMON.MUCA'
338 double precision vcm(3),incr(3)
339 double precision cm(3),L(3)
340 integer ilen,count,rstcount
343 integer maxcount_scale /20/
345 double precision stochforcvec(MAXRES6)
346 common /stochcalc/ stochforcvec
347 integer itime,icount_scale,itime_scal,ifac_time,i,j,itt
349 double precision epdrift,fac_time
356 else if (lang.eq.2 .or. lang.eq.3) then
358 call stochastic_force(stochforcvec)
361 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
363 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
370 icount_scale=icount_scale+1
371 if (icount_scale.gt.maxcount_scale) then
373 & "ERROR: too many attempts at scaling down the time step. ",
374 & "amax=",amax,"epdrift=",epdrift,
375 & "damax=",damax,"edriftmax=",edriftmax,
379 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
383 c First step of the velocity Verlet algorithm
388 else if (lang.eq.3) then
390 call sd_verlet1_ciccotti
392 else if (lang.eq.1) then
397 c Build the chain from the newly calculated coordinates
399 if (rattle) call rattle1
401 if (large.and. mod(itime,ntwe).eq.0) then
402 write (iout,*) "Cartesian and internal coordinates: step 1"
407 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
408 & (dc(j,i+nres),j=1,3)
410 write (iout,*) "Accelerations"
412 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
413 & (d_a(j,i+nres),j=1,3)
415 write (iout,*) "Velocities, step 1"
417 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
418 & (d_t(j,i+nres),j=1,3)
427 c Calculate energy and forces
429 call etotal(potEcomp)
430 ! AL 4/17/17: Reduce the steps if NaNs occurred.
431 if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
436 if (large.and. mod(itime,ntwe).eq.0)
437 & call enerprint(potEcomp)
440 t_etotal=t_etotal+MPI_Wtime()-tt0
442 t_etotal=t_etotal+tcpu()-tt0
445 potE=potEcomp(0)-potEcomp(27)
447 c Get the new accelerations
450 t_enegrad=t_enegrad+MPI_Wtime()-tt0
452 t_enegrad=t_enegrad+tcpu()-tt0
454 c Determine maximum acceleration and scale down the timestep if needed
456 amax=amax/(itime_scal+1)**2
457 call predict_edrift(epdrift)
458 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
459 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
461 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
463 itime_scal=itime_scal+ifac_time
464 c fac_time=dmin1(damax/amax,0.5d0)
465 fac_time=0.5d0**ifac_time
466 d_time=d_time*fac_time
467 if (lang.eq.2 .or. lang.eq.3) then
469 c write (iout,*) "Calling sd_verlet_setup: 1"
470 c Rescale the stochastic forces and recalculate or restore
471 c the matrices of tinker integrator
472 if (itime_scal.gt.maxflag_stoch) then
473 if (large) write (iout,'(a,i5,a)')
474 & "Calculate matrices for stochastic step;",
475 & " itime_scal ",itime_scal
477 call sd_verlet_p_setup
479 call sd_verlet_ciccotti_setup
481 write (iout,'(2a,i3,a,i3,1h.)')
482 & "Warning: cannot store matrices for stochastic",
483 & " integration because the index",itime_scal,
484 & " is greater than",maxflag_stoch
485 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
486 & " integration Langevin algorithm for better efficiency."
487 else if (flag_stoch(itime_scal)) then
488 if (large) write (iout,'(a,i5,a,l1)')
489 & "Restore matrices for stochastic step; itime_scal ",
490 & itime_scal," flag ",flag_stoch(itime_scal)
493 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
494 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
495 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
496 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
497 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
498 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
502 if (large) write (iout,'(2a,i5,a,l1)')
503 & "Calculate & store matrices for stochastic step;",
504 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
506 call sd_verlet_p_setup
508 call sd_verlet_ciccotti_setup
510 flag_stoch(ifac_time)=.true.
513 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
514 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
515 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
516 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
517 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
518 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
522 fac_time=1.0d0/dsqrt(fac_time)
524 stochforcvec(i)=fac_time*stochforcvec(i)
527 else if (lang.eq.1) then
528 c Rescale the accelerations due to stochastic forces
529 fac_time=1.0d0/dsqrt(fac_time)
531 d_as_work(i)=d_as_work(i)*fac_time
534 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
535 & "itime",itime," Timestep scaled down to ",
536 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
538 c Second step of the velocity Verlet algorithm
543 else if (lang.eq.3) then
545 call sd_verlet2_ciccotti
547 else if (lang.eq.1) then
552 if (rattle) call rattle2
555 C print *,totTafm,"TU?"
556 if (d_time.ne.d_time0) then
559 if (lang.eq.2 .or. lang.eq.3) then
560 if (large) write (iout,'(a)')
561 & "Restore original matrices for stochastic step"
562 c write (iout,*) "Calling sd_verlet_setup: 2"
563 c Restore the matrices of tinker integrator if the time step has been restored
566 pfric_mat(i,j)=pfric0_mat(i,j,0)
567 afric_mat(i,j)=afric0_mat(i,j,0)
568 vfric_mat(i,j)=vfric0_mat(i,j,0)
569 prand_mat(i,j)=prand0_mat(i,j,0)
570 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
571 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
580 c Calculate the kinetic and the total energy and the kinetic temperature
585 c write (iout,*) "step",itime," EK",EK," EK1",EK1
587 c Couple the system to Berendsen bath if needed
588 if (tbf .and. lang.eq.0) then
591 kinetic_T=2.0d0/(dimen3*Rb)*EK
592 c Backup the coordinates, velocities, and accelerations
596 d_t_old(j,i)=d_t(j,i)
597 d_a_old(j,i)=d_a(j,i)
601 if (mod(itime,ntwe).eq.0 .and. large) then
602 write (iout,*) "Velocities, step 2"
604 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
605 & (d_t(j,i+nres),j=1,3)
611 c-------------------------------------------------------------------------------
612 subroutine RESPA_step(itime)
613 c-------------------------------------------------------------------------------
614 c Perform a single RESPA step.
615 c-------------------------------------------------------------------------------
620 integer IERROR,ERRCODE
622 include 'COMMON.SETUP'
623 include 'COMMON.CONTROL'
627 include 'COMMON.LAGRANGE.5diag'
629 include 'COMMON.LAGRANGE'
632 include 'COMMON.LANGEVIN'
635 include 'COMMON.LANGEVIN.lang0.5diag'
637 include 'COMMON.LANGEVIN.lang0'
640 include 'COMMON.CHAIN'
641 include 'COMMON.DERIV'
643 include 'COMMON.LOCAL'
644 include 'COMMON.INTERACT'
645 include 'COMMON.IOUNITS'
646 include 'COMMON.NAMES'
647 include 'COMMON.TIME1'
649 common /shortcheck/ lprint_short
650 double precision energia_short(0:n_ene),
651 & energia_long(0:n_ene)
652 double precision cm(3),L(3),vcm(3),incr(3)
653 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
654 & d_a_old0(3,0:maxres2)
656 double precision fac_time
657 logical PRINT_AMTS_MSG /.false./
658 integer ilen,count,rstcount
661 integer maxcount_scale /10/
663 double precision stochforcvec(MAXRES6)
664 common /stochcalc/ stochforcvec
668 common /cipiszcze/ itt
670 double precision epdrift,epdriftmax
674 if (large.and. mod(itime,ntwe).eq.0) then
675 write (iout,*) "***************** RESPA itime",itime
676 write (iout,*) "Cartesian and internal coordinates: step 0"
681 write (iout,*) "Accelerations from long-range forces"
683 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
684 & (d_a(j,i+nres),j=1,3)
686 write (iout,*) "Velocities, step 0"
688 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
689 & (d_t(j,i+nres),j=1,3)
694 c Perform the initial RESPA step (increment velocities)
695 c write (iout,*) "*********************** RESPA ini"
698 if (mod(itime,ntwe).eq.0 .and. large) then
699 write (iout,*) "Velocities, end"
701 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
702 & (d_t(j,i+nres),j=1,3)
706 c Compute the short-range forces
712 C 7/2/2009 commented out
714 c call etotal_short(energia_short)
717 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
720 d_a(j,i)=d_a_short(j,i)
724 if (large.and. mod(itime,ntwe).eq.0) then
725 write (iout,*) "energia_short",energia_short(0)
726 write (iout,*) "Accelerations from short-range forces"
728 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
729 & (d_a(j,i+nres),j=1,3)
734 t_enegrad=t_enegrad+MPI_Wtime()-tt0
736 t_enegrad=t_enegrad+tcpu()-tt0
741 d_t_old(j,i)=d_t(j,i)
742 d_a_old(j,i)=d_a(j,i)
745 c 6/30/08 A-MTS: attempt at increasing the split number
748 dc_old0(j,i)=dc_old(j,i)
749 d_t_old0(j,i)=d_t_old(j,i)
750 d_a_old0(j,i)=d_a_old(j,i)
753 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
754 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
761 c write (iout,*) "itime",itime," ntime_split",ntime_split
762 c Split the time step
763 d_time=d_time0/ntime_split
764 c Perform the short-range RESPA steps (velocity Verlet increments of
765 c positions and velocities using short-range forces)
766 c write (iout,*) "*********************** RESPA split"
767 do itsplit=1,ntime_split
770 else if (lang.eq.2 .or. lang.eq.3) then
772 call stochastic_force(stochforcvec)
775 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
777 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
782 c First step of the velocity Verlet algorithm
787 else if (lang.eq.3) then
789 call sd_verlet1_ciccotti
791 else if (lang.eq.1) then
796 c Build the chain from the newly calculated coordinates
798 if (rattle) call rattle1
800 if (large.and. mod(itime,ntwe).eq.0) then
801 write (iout,*) "***** ITSPLIT",itsplit
802 write (iout,*) "Cartesian and internal coordinates: step 1"
807 write (iout,*) "Velocities, step 1"
809 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
810 & (d_t(j,i+nres),j=1,3)
819 c Calculate energy and forces
820 c if (large.and. mod(itime,ntwe).eq.0) lprint_short=.true.
822 call etotal_short(energia_short)
823 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
824 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) )
826 if (PRINT_AMTS_MSG) write (iout,*)
827 & "Infinities/NaNs in energia_short",
828 & energia_short(0),"; increasing ntime_split to",ntime_split
829 ntime_split=ntime_split*2
830 if (ntime_split.gt.maxtime_split) then
832 write (iout,*) "Cannot rescue the run; aborting job.",
833 & " Retry with a smaller time step"
835 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
837 write (iout,*) "Cannot rescue the run; terminating.",
838 & " Retry with a smaller time step"
844 if (large.and. mod(itime,ntwe).eq.0)
845 & call enerprint(energia_short)
848 t_eshort=t_eshort+MPI_Wtime()-tt0
850 t_eshort=t_eshort+tcpu()-tt0
854 c Get the new accelerations
856 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
859 d_a_short(j,i)=d_a(j,i)
863 if (large.and. mod(itime,ntwe).eq.0) then
864 write (iout,*)"energia_short",energia_short(0)
865 write (iout,*) "Accelerations from short-range forces"
867 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
868 & (d_a(j,i+nres),j=1,3)
873 c Determine maximum acceleration and scale down the timestep if needed
875 amax=amax/ntime_split**2
876 call predict_edrift(epdrift)
877 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
878 & write (iout,*) "amax",amax," damax",damax,
879 & " epdrift",epdrift," epdriftmax",epdriftmax
880 c Exit loop and try with increased split number if the change of
881 c acceleration is too big
882 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
883 if (ntime_split.lt.maxtime_split) then
885 ntime_split=ntime_split*2
888 dc_old(j,i)=dc_old0(j,i)
889 d_t_old(j,i)=d_t_old0(j,i)
890 d_a_old(j,i)=d_a_old0(j,i)
893 if (PRINT_AMTS_MSG) then
894 write (iout,*) "acceleration/energy drift too large",amax,
895 & epdrift," split increased to ",ntime_split," itime",itime,
901 & "Uh-hu. Bumpy landscape. Maximum splitting number",
903 & " already reached!!! Trying to carry on!"
907 t_enegrad=t_enegrad+MPI_Wtime()-tt0
909 t_enegrad=t_enegrad+tcpu()-tt0
911 c Second step of the velocity Verlet algorithm
916 else if (lang.eq.3) then
918 call sd_verlet2_ciccotti
920 else if (lang.eq.1) then
925 if (rattle) call rattle2
926 c Backup the coordinates, velocities, and accelerations
930 d_t_old(j,i)=d_t(j,i)
931 d_a_old(j,i)=d_a(j,i)
938 c Restore the time step
940 c Compute long-range forces
947 call etotal_long(energia_long)
948 ! AL 4/17/2017 Handling NaNs
949 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
952 & "Infinitied/NaNs in energia_long, Aborting MPI job."
953 call enerprint(energia_long)
955 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
957 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
958 call enerprint(energia_long)
963 c lprint_short=.false.
964 if (large.and. mod(itime,ntwe).eq.0)
965 & call enerprint(energia_long)
968 t_elong=t_elong+MPI_Wtime()-tt0
970 t_elong=t_elong+tcpu()-tt0
976 t_enegrad=t_enegrad+MPI_Wtime()-tt0
978 t_enegrad=t_enegrad+tcpu()-tt0
980 c Compute accelerations from long-range forces
982 if (large.and. mod(itime,ntwe).eq.0) then
983 write (iout,*) "energia_long",energia_long(0)
984 write (iout,*) "Cartesian and internal coordinates: step 2"
989 write (iout,*) "Accelerations from long-range forces"
991 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
992 & (d_a(j,i+nres),j=1,3)
994 write (iout,*) "Velocities, step 2"
996 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
997 & (d_t(j,i+nres),j=1,3)
1001 c Compute the final RESPA step (increment velocities)
1002 c write (iout,*) "*********************** RESPA fin"
1004 c Compute the complete potential energy
1006 potEcomp(i)=energia_short(i)+energia_long(i)
1008 potE=potEcomp(0)-potEcomp(27)
1010 if (large.and. mod(itime,ntwe).eq.0) then
1011 call enerprint(potEcomp)
1012 write (iout,*) "potE",potE
1015 c potE=energia_short(0)+energia_long(0)
1018 c Calculate the kinetic and the total energy and the kinetic temperature
1021 c Couple the system to Berendsen bath if needed
1022 if (tbf .and. lang.eq.0) then
1025 kinetic_T=2.0d0/(dimen3*Rb)*EK
1026 c Backup the coordinates, velocities, and accelerations
1028 if (mod(itime,ntwe).eq.0 .and. large) then
1029 write (iout,*) "Velocities, end"
1031 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1032 & (d_t(j,i+nres),j=1,3)
1038 c---------------------------------------------------------------------
1039 subroutine RESPA_vel
1040 c First and last RESPA step (incrementing velocities using long-range
1043 include 'DIMENSIONS'
1044 include 'COMMON.CONTROL'
1045 include 'COMMON.VAR'
1048 include 'COMMON.LAGRANGE.5diag'
1050 include 'COMMON.LAGRANGE'
1052 include 'COMMON.CHAIN'
1053 include 'COMMON.DERIV'
1054 include 'COMMON.GEO'
1055 include 'COMMON.LOCAL'
1056 include 'COMMON.INTERACT'
1057 include 'COMMON.IOUNITS'
1058 include 'COMMON.NAMES'
1061 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1065 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1069 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1072 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1078 c-----------------------------------------------------------------
1080 c Applying velocity Verlet algorithm - step 1 to coordinates
1082 include 'DIMENSIONS'
1083 include 'COMMON.CONTROL'
1084 include 'COMMON.VAR'
1087 include 'COMMON.LAGRANGE.5diag'
1089 include 'COMMON.LAGRANGE'
1091 include 'COMMON.CHAIN'
1092 include 'COMMON.DERIV'
1093 include 'COMMON.GEO'
1094 include 'COMMON.LOCAL'
1095 include 'COMMON.INTERACT'
1096 include 'COMMON.IOUNITS'
1097 include 'COMMON.NAMES'
1098 double precision adt,adt2
1101 write (iout,*) "VELVERLET1 START: DC"
1103 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1104 & (dc(j,i+nres),j=1,3)
1108 adt=d_a_old(j,0)*d_time
1110 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1111 d_t_new(j,0)=d_t_old(j,0)+adt2
1112 d_t(j,0)=d_t_old(j,0)+adt
1118 write (iout,*) "i",i," d_a_old",(d_a_old(j,i),j=1,3)
1121 adt=d_a_old(j,i)*d_time
1123 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1124 d_t_new(j,i)=d_t_old(j,i)+adt2
1125 d_t(j,i)=d_t_old(j,i)+adt
1130 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1133 adt=d_a_old(j,inres)*d_time
1135 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1136 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1137 d_t(j,inres)=d_t_old(j,inres)+adt
1142 write (iout,*) "VELVERLET1 END: DC"
1144 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1145 & (dc(j,i+nres),j=1,3)
1150 c---------------------------------------------------------------------
1152 c Step 2 of the velocity Verlet algorithm: update velocities
1154 include 'DIMENSIONS'
1155 include 'COMMON.CONTROL'
1156 include 'COMMON.VAR'
1159 include 'COMMON.LAGRANGE.5diag'
1161 include 'COMMON.LAGRANGE'
1163 include 'COMMON.CHAIN'
1164 include 'COMMON.DERIV'
1165 include 'COMMON.GEO'
1166 include 'COMMON.LOCAL'
1167 include 'COMMON.INTERACT'
1168 include 'COMMON.IOUNITS'
1169 include 'COMMON.NAMES'
1172 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1176 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1180 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1183 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1189 c-----------------------------------------------------------------
1190 subroutine sddir_precalc
1191 c Applying velocity Verlet algorithm - step 1 to coordinates
1193 include 'DIMENSIONS'
1197 include 'COMMON.CONTROL'
1198 include 'COMMON.VAR'
1201 include 'COMMON.LAGRANGE.5diag'
1203 include 'COMMON.LAGRANGE'
1206 include 'COMMON.LANGEVIN'
1209 include 'COMMON.LANGEVIN.lang0.5diag'
1211 include 'COMMON.LANGEVIN.lang0'
1214 include 'COMMON.CHAIN'
1215 include 'COMMON.DERIV'
1216 include 'COMMON.GEO'
1217 include 'COMMON.LOCAL'
1218 include 'COMMON.INTERACT'
1219 include 'COMMON.IOUNITS'
1220 include 'COMMON.NAMES'
1221 include 'COMMON.TIME1'
1222 double precision time00
1223 double precision stochforcvec(MAXRES6)
1224 common /stochcalc/ stochforcvec
1227 c Compute friction and stochastic forces
1235 c write (iout,*) "After friction_force"
1237 time_fric=time_fric+MPI_Wtime()-time00
1240 time_fric=time_fric+tcpu()-time00
1243 call stochastic_force(stochforcvec)
1244 c write (iout,*) "After stochastic_force"
1246 time_stoch=time_stoch+MPI_Wtime()-time00
1248 time_stoch=time_stoch+tcpu()-time00
1251 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1252 c forces (d_as_work)
1255 c write (iout,*) "friction accelerations"
1256 call fivediaginv_mult(dimen,fric_work, d_af_work)
1257 c write (iout,*) "stochastic acceleratios"
1258 call fivediaginv_mult(dimen,stochforcvec, d_as_work)
1259 c write (iout,*) "Leaving sddir_precalc"
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)
1272 c---------------------------------------------------------------------
1273 subroutine sddir_verlet1
1274 c Applying velocity Verlet algorithm - step 1 to velocities
1276 include 'DIMENSIONS'
1277 include 'COMMON.CONTROL'
1278 include 'COMMON.VAR'
1281 include 'COMMON.LAGRANGE.5diag'
1283 include 'COMMON.LAGRANGE'
1286 include 'COMMON.LANGEVIN'
1289 include 'COMMON.LANGEVIN.lang0.5diag'
1291 include 'COMMON.LANGEVIN.lang0'
1294 include 'COMMON.CHAIN'
1295 include 'COMMON.DERIV'
1296 include 'COMMON.GEO'
1297 include 'COMMON.LOCAL'
1298 include 'COMMON.INTERACT'
1299 include 'COMMON.IOUNITS'
1300 include 'COMMON.NAMES'
1301 c Revised 3/31/05 AL: correlation between random contributions to
1302 c position and velocity increments included.
1303 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1304 double precision adt,adt2
1305 integer i,j,ind,inres
1307 c Add the contribution from BOTH friction and stochastic force to the
1308 c coordinates, but ONLY the contribution from the friction forces to velocities
1311 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1312 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1313 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1314 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1315 d_t(j,0)=d_t_old(j,0)+adt
1320 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1321 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1322 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1323 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1324 d_t(j,i)=d_t_old(j,i)+adt
1329 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1332 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1333 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1334 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1335 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1336 d_t(j,inres)=d_t_old(j,inres)+adt
1343 c---------------------------------------------------------------------
1344 subroutine sddir_verlet2
1345 c Calculating the adjusted velocities for accelerations
1347 include 'DIMENSIONS'
1348 include 'COMMON.CONTROL'
1349 include 'COMMON.VAR'
1352 include 'COMMON.LAGRANGE.5diag'
1354 include 'COMMON.LAGRANGE'
1357 include 'COMMON.LANGEVIN'
1360 include 'COMMON.LANGEVIN.lang0.5diag'
1362 include 'COMMON.LANGEVIN.lang0'
1365 include 'COMMON.CHAIN'
1366 include 'COMMON.DERIV'
1367 include 'COMMON.GEO'
1368 include 'COMMON.LOCAL'
1369 include 'COMMON.INTERACT'
1370 include 'COMMON.IOUNITS'
1371 include 'COMMON.NAMES'
1372 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1373 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1374 integer i,j,inres,ind
1375 c Revised 3/31/05 AL: correlation between random contributions to
1376 c position and velocity increments included.
1377 c The correlation coefficients are calculated at low-friction limit.
1378 c Also, friction forces are now not calculated with new velocities.
1380 c call friction_force
1381 call stochastic_force(stochforcvec)
1383 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1384 c forces (d_as_work)
1387 call fivediaginv_mult(maxres6,stochforcvec, d_as_work1)
1389 call ginv_mult(stochforcvec, d_as_work1)
1395 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1396 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1401 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1402 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1407 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1410 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1411 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1412 & +cos60*d_as_work1(ind+j))*d_time
1419 c---------------------------------------------------------------------
1420 subroutine max_accel
1422 c Find the maximum difference in the accelerations of the the sites
1423 c at the beginning and the end of the time step.
1426 include 'DIMENSIONS'
1427 include 'COMMON.CONTROL'
1428 include 'COMMON.VAR'
1431 include 'COMMON.LAGRANGE.5diag'
1433 include 'COMMON.LAGRANGE'
1435 include 'COMMON.CHAIN'
1436 include 'COMMON.DERIV'
1437 include 'COMMON.GEO'
1438 include 'COMMON.LOCAL'
1439 include 'COMMON.INTERACT'
1440 include 'COMMON.IOUNITS'
1442 double precision aux(3),accel(3),accel_old(3),dacc
1444 c aux(j)=d_a(j,0)-d_a_old(j,0)
1445 accel_old(j)=d_a_old(j,0)
1452 c 7/3/08 changed to asymmetric difference
1454 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1455 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1456 accel(j)=accel(j)+0.5d0*d_a(j,i)
1457 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1458 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1459 dacc=dabs(accel(j)-accel_old(j))
1460 c write (iout,*) i,dacc
1461 if (dacc.gt.amax) amax=dacc
1469 accel_old(j)=d_a_old(j,0)
1474 accel_old(j)=accel_old(j)+d_a_old(j,1)
1475 accel(j)=accel(j)+d_a(j,1)
1479 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1481 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1482 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1483 accel(j)=accel(j)+d_a(j,i+nres)
1487 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1488 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1489 dacc=dabs(accel(j)-accel_old(j))
1490 c write (iout,*) "side-chain",i,dacc
1491 if (dacc.gt.amax) amax=dacc
1495 accel_old(j)=accel_old(j)+d_a_old(j,i)
1496 accel(j)=accel(j)+d_a(j,i)
1497 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1502 c---------------------------------------------------------------------
1503 subroutine predict_edrift(epdrift)
1505 c Predict the drift of the potential energy
1508 include 'DIMENSIONS'
1509 include 'COMMON.CONTROL'
1510 include 'COMMON.VAR'
1513 include 'COMMON.LAGRANGE.5diag'
1515 include 'COMMON.LAGRANGE'
1517 include 'COMMON.CHAIN'
1518 include 'COMMON.DERIV'
1519 include 'COMMON.GEO'
1520 include 'COMMON.LOCAL'
1521 include 'COMMON.INTERACT'
1522 include 'COMMON.IOUNITS'
1523 include 'COMMON.MUCA'
1524 double precision epdrift,epdriftij
1526 c Drift of the potential energy
1532 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1533 if (lmuca) epdriftij=epdriftij*factor
1534 c write (iout,*) "back",i,j,epdriftij
1535 if (epdriftij.gt.epdrift) epdrift=epdriftij
1539 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1542 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1543 if (lmuca) epdriftij=epdriftij*factor
1544 c write (iout,*) "side",i,j,epdriftij
1545 if (epdriftij.gt.epdrift) epdrift=epdriftij
1549 epdrift=0.5d0*epdrift*d_time*d_time
1550 c write (iout,*) "epdrift",epdrift
1553 c-----------------------------------------------------------------------
1554 subroutine verlet_bath
1556 c Coupling to the thermostat by using the Berendsen algorithm
1559 include 'DIMENSIONS'
1560 include 'COMMON.CONTROL'
1561 include 'COMMON.VAR'
1564 include 'COMMON.LAGRANGE.5diag'
1566 include 'COMMON.LAGRANGE'
1570 include 'COMMON.LANGEVIN.lang0.5diag'
1572 include 'COMMON.LANGEVIN.lang0'
1575 include 'COMMON.LANGEVIN'
1577 include 'COMMON.CHAIN'
1578 include 'COMMON.DERIV'
1579 include 'COMMON.GEO'
1580 include 'COMMON.LOCAL'
1581 include 'COMMON.INTERACT'
1582 include 'COMMON.IOUNITS'
1583 include 'COMMON.NAMES'
1584 double precision T_half,fact
1586 T_half=2.0d0/(dimen3*Rb)*EK
1587 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1588 c write(iout,*) "T_half", T_half
1589 c write(iout,*) "EK", EK
1590 c write(iout,*) "fact", fact
1592 d_t(j,0)=fact*d_t(j,0)
1596 d_t(j,i)=fact*d_t(j,i)
1600 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1603 d_t(j,inres)=fact*d_t(j,inres)
1609 c---------------------------------------------------------
1611 c Set up the initial conditions of a MD simulation
1613 include 'DIMENSIONS'
1617 integer IERROR,ERRCODE,error_msg,ierr,ierrcode
1619 include 'COMMON.SETUP'
1620 include 'COMMON.CONTROL'
1621 include 'COMMON.VAR'
1624 include 'COMMON.LAGRANGE.5diag'
1626 include 'COMMON.LAGRANGE'
1628 include 'COMMON.QRESTR'
1630 include 'COMMON.LANGEVIN'
1633 include 'COMMON.LANGEVIN.lang0.5diag'
1636 include 'COMMON.LANGEVIN.lang0.5diag'
1638 include 'COMMON.LANGEVIN.lang0'
1642 include 'COMMON.CHAIN'
1643 include 'COMMON.DERIV'
1644 include 'COMMON.GEO'
1645 include 'COMMON.LOCAL'
1646 include 'COMMON.INTERACT'
1647 include 'COMMON.IOUNITS'
1648 include 'COMMON.NAMES'
1649 include 'COMMON.REMD'
1650 include 'COMMON.TIME1'
1654 common /lbfgstat/ status,niter,nfun
1656 integer n_model_try,list_model_try(max_template),k
1657 double precision tt0
1658 real*8 energia_long(0:n_ene),
1659 & energia_short(0:n_ene),vcm(3),incr(3)
1660 double precision cm(3),L(3),xv,sigv,lowb,highb
1661 double precision varia(maxvar),energia(0:n_ene)
1668 integer i,ipos,iq,iw,j,iranmin,nft_sc,iretcode,nfun,itrial,itmp,
1671 double precision etot
1673 write (iout,*) "init_MD INDPDB",indpdb
1675 c write(iout,*) "d_time", d_time
1676 c Compute the standard deviations of stochastic forces for Langevin dynamics
1677 c if the friction coefficients do not depend on surface area
1678 if (lang.gt.0 .and. .not.surfarea) then
1680 stdforcp(i)=stdfp*dsqrt(gamp)
1683 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1684 & *dsqrt(gamsc(iabs(itype(i))))
1687 c Open the pdb file for snapshotshots
1690 if (ilen(tmpdir).gt.0)
1691 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1692 & liczba(:ilen(liczba))//".pdb")
1694 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1698 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1699 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1700 & liczba(:ilen(liczba))//".x")
1701 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1704 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1705 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1706 & liczba(:ilen(liczba))//".cx")
1707 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1713 if (ilen(tmpdir).gt.0)
1714 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1715 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1717 if (ilen(tmpdir).gt.0)
1718 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1719 cartname=prefix(:ilen(prefix))//"_MD.cx"
1723 write (qstr,'(256(1h ))')
1726 iq = qinfrag(i,iset)*10
1727 iw = wfrag(i,iset)/100
1729 if(me.eq.king.or..not.out1file)
1730 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1731 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1736 iq = qinpair(i,iset)*10
1737 iw = wpair(i,iset)/100
1739 if(me.eq.king.or..not.out1file)
1740 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1741 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1745 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1747 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1749 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1751 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1754 write (iout,*) "REST ",rest
1756 if (restart1file) then
1758 & inquire(file=mremd_rst_name,exist=file_exist)
1760 write (*,*) me," Before broadcast: file_exist",file_exist
1761 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1763 write (*,*) me," After broadcast: file_exist",file_exist
1764 c inquire(file=mremd_rst_name,exist=file_exist)
1766 if(me.eq.king.or..not.out1file)
1767 & write(iout,*) "Initial state read by master and distributed"
1769 if (ilen(tmpdir).gt.0)
1770 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1771 & //liczba(:ilen(liczba))//'.rst')
1772 inquire(file=rest2name,exist=file_exist)
1775 if(.not.restart1file) then
1776 if(me.eq.king.or..not.out1file)
1777 & write(iout,*) "Initial state will be read from file ",
1778 & rest2name(:ilen(rest2name))
1781 call rescale_weights(t_bath)
1784 if(me.eq.king.or..not.out1file)then
1785 if (restart1file) then
1786 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1789 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1792 write(iout,*) "Initial velocities randomly generated"
1799 c Generate initial velocities
1800 if(me.eq.king.or..not.out1file)
1801 & write(iout,*) "Initial velocities randomly generated"
1804 CtotTafm is the variable for AFM time which eclipsed during
1807 c rest2name = prefix(:ilen(prefix))//'.rst'
1808 if(me.eq.king.or..not.out1file)then
1809 write (iout,*) "Initial velocities"
1811 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1812 & (d_t(j,i+nres),j=1,3)
1814 c Zeroing the total angular momentum of the system
1815 write(iout,*) "Calling the zero-angular
1816 & momentum subroutine"
1819 c Getting the potential energy and forces and velocities and accelerations
1821 write (iout,*) "velocity of the center of the mass:"
1822 write (iout,*) (vcm(j),j=1,3)
1825 d_t(j,0)=d_t(j,0)-vcm(j)
1827 c Removing the velocity of the center of mass
1829 if(me.eq.king.or..not.out1file)then
1830 write (iout,*) "vcm right after adjustment:"
1831 write (iout,*) (vcm(j),j=1,3)
1834 write (iout,*) "init_MD before initial structure REST ",rest
1837 if (iranconf.ne.0) then
1838 c 8/22/17 AL Loop to produce a low-energy random conformation
1842 print *, 'Calling OVERLAP_SC'
1843 call overlap_sc(fail)
1847 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1848 print *,'SC_move',nft_sc,etot
1849 if(me.eq.king.or..not.out1file)
1850 & write(iout,*) 'SC_move',nft_sc,etot
1854 if (me.eq.king.or..not.out1file) write(iout,*)
1855 & 'Minimizing random structure: Calling MINIM_DC'
1856 call minim_dc(etot,iretcode,nfun)
1858 call geom_to_var(nvar,varia)
1859 if (me.eq.king.or..not.out1file) write (iout,*)
1860 & 'Minimizing random structure: Calling MINIMIZE.'
1861 call minimize(etot,varia,iretcode,nfun)
1862 call var_to_geom(nvar,varia)
1864 if (me.eq.king.or..not.out1file)
1865 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1866 if (isnan(etot) .or. etot.gt.1.0d4) then
1867 write (iout,*) "Energy too large",etot,
1868 & " trying another random conformation"
1871 call gen_rand_conf(itmp,*30)
1873 30 write (iout,*) 'Failed to generate random conformation',
1874 & ', itrial=',itrial
1875 write (*,*) 'Processor:',me,
1876 & ' Failed to generate random conformation',
1886 write (iout,'(a,i3,a)') 'Processor:',me,
1887 & ' error in generating random conformation.'
1888 write (*,'(a,i3,a)') 'Processor:',me,
1889 & ' error in generating random conformation.'
1892 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1901 write (iout,'(a,i3,a)') 'Processor:',me,
1902 & ' failed to generate a low-energy random conformation.'
1903 write (*,'(a,i3,a)') 'Processor:',me,
1904 & ' failed to generate a low-energy random conformation.'
1907 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1912 else if (preminim) then
1913 if (start_from_model) then
1915 do while (fail .and. n_model_try.lt.constr_homology)
1917 i_model=iran_num(1,constr_homology)
1919 if (i_model.eq.list_model_try(k)) exit
1921 if (k.gt.n_model_try) exit
1923 n_model_try=n_model_try+1
1924 list_model_try(n_model_try)=i_model
1925 write (iout,*) 'starting from model ',i_model
1928 c(j,i)=chomo(j,i,i_model)
1931 call int_from_cart(.true.,.false.)
1932 call sc_loc_geom(.false.)
1935 dc(j,i)=c(j,i+1)-c(j,i)
1936 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1941 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1942 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1945 if (me.eq.king.or..not.out1file) then
1946 write (iout,*) "Energies before removing overlaps"
1947 call etotal(energia(0))
1948 call enerprint(energia(0))
1950 ! Remove SC overlaps if requested
1952 write (iout,*) 'Calling OVERLAP_SC'
1953 call overlap_sc(fail)
1956 & "Failed to remove overlap from model",i_model
1961 if (me.eq.king.or..not.out1file) then
1962 write (iout,*) "Energies after removing overlaps"
1963 call etotal(energia(0))
1964 call enerprint(energia(0))
1966 ! Search for better SC rotamers if requested
1968 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1969 print *,'SC_move',nft_sc,etot
1970 if (me.eq.king.or..not.out1file)
1971 & write(iout,*) 'SC_move',nft_sc,etot
1973 call etotal(energia(0))
1976 if (n_model_try.gt.constr_homology) then
1978 & "All models have irreparable overlaps. Trying randoms starts."
1983 C 8/22/17 AL Minimize initial structure
1985 if (me.eq.king.or..not.out1file) write(iout,*)
1986 & 'Minimizing initial PDB structure: Calling MINIM_DC'
1987 call minim_dc(etot,iretcode,nfun)
1989 call geom_to_var(nvar,varia)
1990 if(me.eq.king.or..not.out1file) write (iout,*)
1991 & 'Minimizing initial PDB structure: Calling MINIMIZE.'
1992 call minimize(etot,varia,iretcode,nfun)
1993 call var_to_geom(nvar,varia)
1995 if (me.eq.king.or..not.out1file)
1996 & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
1997 if(me.eq.king.or..not.out1file)
1998 & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2000 if (me.eq.king.or..not.out1file)
2001 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2002 if(me.eq.king.or..not.out1file)
2003 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2008 call chainbuild_cart
2013 kinetic_T=2.0d0/(dimen3*Rb)*EK
2014 if(me.eq.king.or..not.out1file)then
2024 call etotal(potEcomp)
2025 if (large) call enerprint(potEcomp)
2028 t_etotal=t_etotal+MPI_Wtime()-tt0
2030 t_etotal=t_etotal+tcpu()-tt0
2034 c write (iout,*) "PotE-homology",potE-potEcomp(27)
2038 if (amax*d_time .gt. dvmax) then
2039 d_time=d_time*dvmax/amax
2040 if(me.eq.king.or..not.out1file) write (iout,*)
2041 & "Time step reduced to",d_time,
2042 & " because of too large initial acceleration."
2044 if(me.eq.king.or..not.out1file)then
2045 write(iout,*) "Potential energy and its components"
2046 call enerprint(potEcomp)
2047 c write(iout,*) (potEcomp(i),i=0,n_ene)
2049 potE=potEcomp(0)-potEcomp(27)
2050 c write (iout,*) "PotE-homology",potE
2053 if (ntwe.ne.0) call statout(itime)
2054 if(me.eq.king.or..not.out1file)
2055 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
2056 & " Kinetic energy",EK," Potential energy",potE,
2057 & " Total energy",totE," Maximum acceleration ",
2060 write (iout,*) "Initial coordinates"
2062 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
2063 & (c(j,i+nres),j=1,3)
2065 write (iout,*) "Initial dC"
2067 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
2068 & (dc(j,i+nres),j=1,3)
2070 write (iout,*) "Initial velocities"
2071 write (iout,"(13x,' backbone ',23x,' side chain')")
2073 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
2074 & (d_t(j,i+nres),j=1,3)
2076 write (iout,*) "Initial accelerations"
2078 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2079 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
2080 & (d_a(j,i+nres),j=1,3)
2086 d_t_old(j,i)=d_t(j,i)
2087 d_a_old(j,i)=d_a(j,i)
2089 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2098 call etotal_short(energia_short)
2099 if (large) call enerprint(potEcomp)
2102 t_eshort=t_eshort+MPI_Wtime()-tt0
2104 t_eshort=t_eshort+tcpu()-tt0
2109 if(.not.out1file .and. large) then
2110 write (iout,*) "energia_long",energia_long(0),
2111 & " energia_short",energia_short(0),
2112 & " total",energia_long(0)+energia_short(0)
2113 write (iout,*) "Initial fast-force accelerations"
2115 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2116 & (d_a(j,i+nres),j=1,3)
2119 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2122 d_a_short(j,i)=d_a(j,i)
2131 call etotal_long(energia_long)
2132 if (large) call enerprint(potEcomp)
2135 t_elong=t_elong+MPI_Wtime()-tt0
2137 t_elong=t_elong+tcpu()-tt0
2142 if(.not.out1file .and. large) then
2143 write (iout,*) "energia_long",energia_long(0)
2144 write (iout,*) "Initial slow-force accelerations"
2146 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2147 & (d_a(j,i+nres),j=1,3)
2151 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2153 t_enegrad=t_enegrad+tcpu()-tt0
2158 c-----------------------------------------------------------
2159 subroutine random_vel
2161 include 'DIMENSIONS'
2162 include 'COMMON.CONTROL'
2163 include 'COMMON.VAR'
2166 include 'COMMON.LAGRANGE.5diag'
2168 include 'COMMON.LAGRANGE'
2171 include 'COMMON.LANGEVIN'
2174 include 'COMMON.LANGEVIN.lang0.5diag'
2176 include 'COMMON.LANGEVIN.lang0'
2179 include 'COMMON.CHAIN'
2180 include 'COMMON.DERIV'
2181 include 'COMMON.GEO'
2182 include 'COMMON.LOCAL'
2183 include 'COMMON.INTERACT'
2184 include 'COMMON.IOUNITS'
2185 include 'COMMON.NAMES'
2186 include 'COMMON.TIME1'
2187 double precision xv,sigv,lowb,highb,vec_afm(3),Ek1,Ek2,Ek3,aux
2188 integer i,ii,j,k,l,ind
2189 double precision anorm_distr
2190 logical lprn /.true./
2192 integer ichain,n,innt,inct,ibeg,ierr
2193 double precision work(8*maxres6)
2194 integer iwork(maxres6)
2195 double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),
2196 & Gvec(maxres2_chain,maxres2_chain)
2197 common /przechowalnia/Ghalf,Geigen,Gvec
2199 double precision inertia(maxres2_chain,maxres2_chain)
2201 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2202 c First generate velocities in the eigenspace of the G matrix
2203 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2206 write (iout,*) "Random_vel, fivediag"
2215 n=dimen_chain(ichain)
2216 innt=iposd_chain(ichain)
2219 write (iout,*) "Chain",ichain," n",n," start",innt
2221 if (i.lt.inct-1) then
2222 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),
2224 else if (i.eq.inct-1) then
2225 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
2227 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
2231 ghalf(ind+1)=dmorig(innt)
2232 ghalf(ind+2)=du1orig(innt)
2233 ghalf(ind+3)=dmorig(innt+1)
2237 c write (iout,*) "i",i," ind",ind," indu2",innt+i-2,
2238 c & " indu1",innt+i-1," indm",innt+i
2239 ghalf(ind+1)=du2orig(innt-1+i-2)
2240 ghalf(ind+2)=du1orig(innt-1+i-1)
2241 ghalf(ind+3)=dmorig(innt-1+i)
2242 c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
2243 c & "DU1",innt-1+i-1,"DM ",innt-1+i
2251 inertia(i,j)=ghalf(ind)
2252 inertia(j,i)=ghalf(ind)
2257 write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
2258 write (iout,*) "Five-diagonal inertia matrix, lower triangle"
2259 call matoutr(n,ghalf)
2261 call gldiag(maxres2_chain,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
2263 write (iout,'(//a,i3)')
2264 & "Eigenvectors and eigenvalues of the G matrix chain",ichain
2265 call eigout(n,n,maxres2_chain,maxres2_chain,Gvec,Geigen)
2268 c check diagonalization
2274 aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
2278 write (iout,*) i,j,aux,geigen(i)
2280 write (iout,*) i,j,aux
2290 sigv=dsqrt((Rb*t_bath)/geigen(i))
2293 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2294 EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
2295 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2296 c & " d_t_work_new",d_t_work_new(ii)
2304 d_t_work(ind)=d_t_work(ind)
2305 & +Gvec(i,j)*d_t_work_new((j-1)*3+k)
2307 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2316 aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
2322 c Transfer to the d_t vector
2323 innt=chain_border(1,ichain)
2324 inct=chain_border(2,ichain)
2326 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
2330 d_t(j,i)=d_t_work(ind)
2332 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2335 d_t(j,i+nres)=d_t_work(ind)
2342 write (iout,*) "Random velocities in the Calpha,SC space"
2344 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2345 & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2348 call kinetic_CASC(Ek1)
2350 ! Transform the velocities to virtual-bond space
2355 if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
2357 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2361 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2362 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2369 c write (iout,*) "Shifting accelerations"
2371 c write (iout,*) "ichain",chain_border1(1,ichain)-1,
2372 c & chain_border1(1,ichain)
2373 d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
2374 d_t(:,chain_border1(1,ichain))=0.0d0
2376 c write (iout,*) "Adding accelerations"
2378 c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2379 c & chain_border(2,ichain-1)
2380 d_t(:,chain_border1(1,ichain)-1)=
2381 & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2382 d_t(:,chain_border(2,ichain-1))=0.0d0
2385 write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2386 & chain_border(2,ichain-1)
2387 d_t(:,chain_border1(1,ichain)-1)=
2388 & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2389 d_t(:,chain_border(2,ichain-1))=0.0d0
2394 c d_t(j,0)=d_t(j,nnt)
2397 innt=chain_border(1,ichain)
2398 inct=chain_border(2,ichain)
2399 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
2400 c write (iout,*) "ibeg",ibeg
2402 d_t(j,ibeg)=d_t(j,innt)
2406 if (iabs(itype(i).eq.10)) then
2407 c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
2409 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2413 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2414 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2423 & "Random velocities in the virtual-bond-vector space"
2424 write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
2426 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2427 & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2430 write (iout,*) "Kinetic energy from inertia matrix eigenvalues",
2433 & "Kinetic temperatures from inertia matrix eigenvalues",
2436 write (iout,*) "Kinetic energy from inertia matrix",Ek3
2437 write (iout,*) "Kinetic temperatures from inertia",
2438 & 2*Ek3/(3*dimen*Rb)
2440 write (iout,*) "Kinetic energy from velocities in CA-SC space",
2443 & "Kinetic temperatures from velovities in CA-SC space",
2444 & 2*Ek1/(3*dimen*Rb)
2447 & "Kinetic energy from virtual-bond-vector velocities",Ek1
2449 & "Kinetic temperature from virtual-bond-vector velocities ",
2458 sigv=dsqrt((Rb*t_bath)/geigen(i))
2461 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2463 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2464 c & " d_t_work_new",d_t_work_new(ii)
2467 C if (SELFGUIDE.gt.0) then
2470 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
2471 C distance=distance+vec_afm(j)**2
2473 C distance=dsqrt(distance)
2475 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
2476 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
2477 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
2478 C & d_t_work_new(j+(afmend-1)*3)
2489 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2492 c write (iout,*) "Ek from eigenvectors",Ek1
2494 c Transform velocities to UNRES coordinate space
2500 d_t_work(ind)=d_t_work(ind)
2501 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2503 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2507 c Transfer to the d_t vector
2509 d_t(j,0)=d_t_work(j)
2515 d_t(j,i)=d_t_work(ind)
2519 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2522 d_t(j,i+nres)=d_t_work(ind)
2527 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2528 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2534 c-----------------------------------------------------------
2535 subroutine sd_verlet_p_setup
2536 c Sets up the parameters of stochastic Verlet algorithm
2538 include 'DIMENSIONS'
2542 include 'COMMON.CONTROL'
2543 include 'COMMON.VAR'
2546 include 'COMMON.LANGEVIN'
2549 include 'COMMON.LANGEVIN.lang0.5diag'
2551 include 'COMMON.LANGEVIN.lang0'
2554 include 'COMMON.CHAIN'
2555 include 'COMMON.DERIV'
2556 include 'COMMON.GEO'
2557 include 'COMMON.LOCAL'
2558 include 'COMMON.INTERACT'
2559 include 'COMMON.IOUNITS'
2560 include 'COMMON.NAMES'
2561 include 'COMMON.TIME1'
2562 double precision emgdt(MAXRES6),
2563 & pterm,vterm,rho,rhoc,vsig,
2564 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2565 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2566 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2567 logical lprn /.false./
2568 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2569 double precision ktm
2576 c AL 8/17/04 Code adapted from tinker
2578 c Get the frictional and random terms for stochastic dynamics in the
2579 c eigenspace of mass-scaled UNRES friction matrix
2582 gdt = fricgam(i) * d_time
2584 c Stochastic dynamics reduces to simple MD for zero friction
2586 if (gdt .le. zero) then
2587 pfric_vec(i) = 1.0d0
2588 vfric_vec(i) = d_time
2589 afric_vec(i) = 0.5d0 * d_time * d_time
2590 prand_vec(i) = 0.0d0
2591 vrand_vec1(i) = 0.0d0
2592 vrand_vec2(i) = 0.0d0
2594 c Analytical expressions when friction coefficient is large
2597 if (gdt .ge. gdt_radius) then
2600 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2601 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2602 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2603 vterm = 1.0d0 - egdt**2
2604 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2606 c Use series expansions when friction coefficient is small
2617 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2618 & - gdt5/120.0d0 + gdt6/720.0d0
2619 & - gdt7/5040.0d0 + gdt8/40320.0d0
2620 & - gdt9/362880.0d0) / fricgam(i)**2
2621 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2622 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2623 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2624 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2625 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2626 & + 127.0d0*gdt9/90720.0d0
2627 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2628 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2629 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2630 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2631 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2632 & - 17.0d0*gdt2/1280.0d0
2633 & + 17.0d0*gdt3/6144.0d0
2634 & + 40967.0d0*gdt4/34406400.0d0
2635 & - 57203.0d0*gdt5/275251200.0d0
2636 & - 1429487.0d0*gdt6/13212057600.0d0)
2639 c Compute the scaling factors of random terms for the nonzero friction case
2641 ktm = 0.5d0*d_time/fricgam(i)
2642 psig = dsqrt(ktm*pterm) / fricgam(i)
2643 vsig = dsqrt(ktm*vterm)
2644 rhoc = dsqrt(1.0d0 - rho*rho)
2646 vrand_vec1(i) = vsig * rho
2647 vrand_vec2(i) = vsig * rhoc
2652 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2655 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2656 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2660 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2662 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2663 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2664 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2665 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2666 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2667 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2669 t_sdsetup=t_sdsetup+MPI_Wtime()
2671 t_sdsetup=t_sdsetup+tcpu()-tt0
2675 c-------------------------------------------------------------
2676 subroutine eigtransf1(n,ndim,ab,d,c)
2679 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2685 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2691 c-------------------------------------------------------------
2692 subroutine eigtransf(n,ndim,a,b,d,c)
2695 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2701 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2707 c-------------------------------------------------------------
2708 subroutine sd_verlet1
2709 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2711 include 'DIMENSIONS'
2712 include 'COMMON.CONTROL'
2713 include 'COMMON.VAR'
2716 include 'COMMON.LAGRANGE.5diag'
2718 include 'COMMON.LAGRANGE'
2721 include 'COMMON.LANGEVIN'
2724 include 'COMMON.LANGEVIN.lang0.5diag'
2726 include 'COMMON.LANGEVIN.lang0'
2729 include 'COMMON.CHAIN'
2730 include 'COMMON.DERIV'
2731 include 'COMMON.GEO'
2732 include 'COMMON.LOCAL'
2733 include 'COMMON.INTERACT'
2734 include 'COMMON.IOUNITS'
2735 include 'COMMON.NAMES'
2736 double precision stochforcvec(MAXRES6)
2737 common /stochcalc/ stochforcvec
2738 logical lprn /.false./
2739 integer i,j,ind,inres
2741 c write (iout,*) "dc_old"
2743 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2744 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2747 dc_work(j)=dc_old(j,0)
2748 d_t_work(j)=d_t_old(j,0)
2749 d_a_work(j)=d_a_old(j,0)
2754 dc_work(ind+j)=dc_old(j,i)
2755 d_t_work(ind+j)=d_t_old(j,i)
2756 d_a_work(ind+j)=d_a_old(j,i)
2761 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2763 dc_work(ind+j)=dc_old(j,i+nres)
2764 d_t_work(ind+j)=d_t_old(j,i+nres)
2765 d_a_work(ind+j)=d_a_old(j,i+nres)
2773 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2777 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2778 & vfric_mat(i,j),afric_mat(i,j),
2779 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2787 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2788 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2789 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2790 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2792 d_t_work_new(i)=ddt1+0.5d0*ddt2
2793 d_t_work(i)=ddt1+ddt2
2798 d_t(j,0)=d_t_work(j)
2803 dc(j,i)=dc_work(ind+j)
2804 d_t(j,i)=d_t_work(ind+j)
2809 if (itype(i).ne.10) then
2812 dc(j,inres)=dc_work(ind+j)
2813 d_t(j,inres)=d_t_work(ind+j)
2820 c--------------------------------------------------------------------------
2821 subroutine sd_verlet2
2822 c Calculating the adjusted velocities for accelerations
2824 include 'DIMENSIONS'
2825 include 'COMMON.CONTROL'
2826 include 'COMMON.VAR'
2829 include 'COMMON.LAGRANGE.5diag'
2831 include 'COMMON.LAGRANGE'
2834 include 'COMMON.LANGEVIN'
2837 include 'COMMON.LANGEVIN.lang0.5diag'
2839 include 'COMMON.LANGEVIN.lang0'
2842 include 'COMMON.CHAIN'
2843 include 'COMMON.DERIV'
2844 include 'COMMON.GEO'
2845 include 'COMMON.LOCAL'
2846 include 'COMMON.INTERACT'
2847 include 'COMMON.IOUNITS'
2848 include 'COMMON.NAMES'
2849 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2850 common /stochcalc/ stochforcvec
2851 integer i,j,ind,inres
2853 c Compute the stochastic forces which contribute to velocity change
2855 call stochastic_force(stochforcvecV)
2862 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2863 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2864 & vrand_mat2(i,j)*stochforcvecV(j)
2866 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2870 d_t(j,0)=d_t_work(j)
2875 d_t(j,i)=d_t_work(ind+j)
2880 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2883 d_t(j,inres)=d_t_work(ind+j)
2890 c-----------------------------------------------------------
2891 subroutine sd_verlet_ciccotti_setup
2892 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2895 include 'DIMENSIONS'
2899 include 'COMMON.CONTROL'
2900 include 'COMMON.VAR'
2903 include 'COMMON.LAGRANGE.5diag'
2905 include 'COMMON.LAGRANGE'
2907 include 'COMMON.LANGEVIN'
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 include 'COMMON.TIME1'
2916 double precision emgdt(MAXRES6),
2917 & pterm,vterm,rho,rhoc,vsig,
2918 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2919 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2920 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2921 logical lprn /.false./
2922 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2923 double precision ktm
2931 c AL 8/17/04 Code adapted from tinker
2933 c Get the frictional and random terms for stochastic dynamics in the
2934 c eigenspace of mass-scaled UNRES friction matrix
2937 write (iout,*) "i",i," fricgam",fricgam(i)
2938 gdt = fricgam(i) * d_time
2940 c Stochastic dynamics reduces to simple MD for zero friction
2942 if (gdt .le. zero) then
2943 pfric_vec(i) = 1.0d0
2944 vfric_vec(i) = d_time
2945 afric_vec(i) = 0.5d0*d_time*d_time
2946 prand_vec(i) = afric_vec(i)
2947 vrand_vec2(i) = vfric_vec(i)
2949 c Analytical expressions when friction coefficient is large
2954 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2955 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2956 prand_vec(i) = afric_vec(i)
2957 vrand_vec2(i) = vfric_vec(i)
2959 c Compute the scaling factors of random terms for the nonzero friction case
2961 c ktm = 0.5d0*d_time/fricgam(i)
2962 c psig = dsqrt(ktm*pterm) / fricgam(i)
2963 c vsig = dsqrt(ktm*vterm)
2964 c prand_vec(i) = psig*afric_vec(i)
2965 c vrand_vec2(i) = vsig*vfric_vec(i)
2970 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2973 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2974 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2978 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2980 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2981 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2982 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2983 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2984 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2986 t_sdsetup=t_sdsetup+MPI_Wtime()
2988 t_sdsetup=t_sdsetup+tcpu()-tt0
2992 c-------------------------------------------------------------
2993 subroutine sd_verlet1_ciccotti
2994 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2996 include 'DIMENSIONS'
3000 include 'COMMON.CONTROL'
3001 include 'COMMON.VAR'
3004 include 'COMMON.LAGRANGE.5diag'
3006 include 'COMMON.LAGRANGE'
3008 include 'COMMON.LANGEVIN'
3009 include 'COMMON.CHAIN'
3010 include 'COMMON.DERIV'
3011 include 'COMMON.GEO'
3012 include 'COMMON.LOCAL'
3013 include 'COMMON.INTERACT'
3014 include 'COMMON.IOUNITS'
3015 include 'COMMON.NAMES'
3016 double precision stochforcvec(MAXRES6)
3017 common /stochcalc/ stochforcvec
3018 logical lprn /.false./
3021 c write (iout,*) "dc_old"
3023 c write (iout,'(i5,3f10.5,5x,3f10.5)')
3024 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3027 dc_work(j)=dc_old(j,0)
3028 d_t_work(j)=d_t_old(j,0)
3029 d_a_work(j)=d_a_old(j,0)
3034 dc_work(ind+j)=dc_old(j,i)
3035 d_t_work(ind+j)=d_t_old(j,i)
3036 d_a_work(ind+j)=d_a_old(j,i)
3041 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3043 dc_work(ind+j)=dc_old(j,i+nres)
3044 d_t_work(ind+j)=d_t_old(j,i+nres)
3045 d_a_work(ind+j)=d_a_old(j,i+nres)
3054 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
3058 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
3059 & vfric_mat(i,j),afric_mat(i,j),
3060 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3068 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
3069 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3070 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3071 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3073 d_t_work_new(i)=ddt1+0.5d0*ddt2
3074 d_t_work(i)=ddt1+ddt2
3079 d_t(j,0)=d_t_work(j)
3084 dc(j,i)=dc_work(ind+j)
3085 d_t(j,i)=d_t_work(ind+j)
3090 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3093 dc(j,inres)=dc_work(ind+j)
3094 d_t(j,inres)=d_t_work(ind+j)
3101 c--------------------------------------------------------------------------
3102 subroutine sd_verlet2_ciccotti
3103 c Calculating the adjusted velocities for accelerations
3105 include 'DIMENSIONS'
3106 include 'COMMON.CONTROL'
3107 include 'COMMON.VAR'
3110 include 'COMMON.LAGRANGE.5diag'
3112 include 'COMMON.LAGRANGE'
3114 include 'COMMON.LANGEVIN'
3115 include 'COMMON.CHAIN'
3116 include 'COMMON.DERIV'
3117 include 'COMMON.GEO'
3118 include 'COMMON.LOCAL'
3119 include 'COMMON.INTERACT'
3120 include 'COMMON.IOUNITS'
3121 include 'COMMON.NAMES'
3122 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
3123 common /stochcalc/ stochforcvec
3126 c Compute the stochastic forces which contribute to velocity change
3128 call stochastic_force(stochforcvecV)
3135 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3136 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3137 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3139 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3143 d_t(j,0)=d_t_work(j)
3148 d_t(j,i)=d_t_work(ind+j)
3153 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3156 d_t(j,inres)=d_t_work(ind+j)