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 adt=d_a_old(j,i)*d_time
1120 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1121 d_t_new(j,i)=d_t_old(j,i)+adt2
1122 d_t(j,i)=d_t_old(j,i)+adt
1127 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1130 adt=d_a_old(j,inres)*d_time
1132 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1133 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1134 d_t(j,inres)=d_t_old(j,inres)+adt
1139 write (iout,*) "VELVERLET1 END: DC"
1141 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1142 & (dc(j,i+nres),j=1,3)
1147 c---------------------------------------------------------------------
1149 c Step 2 of the velocity Verlet algorithm: update velocities
1151 include 'DIMENSIONS'
1152 include 'COMMON.CONTROL'
1153 include 'COMMON.VAR'
1156 include 'COMMON.LAGRANGE.5diag'
1158 include 'COMMON.LAGRANGE'
1160 include 'COMMON.CHAIN'
1161 include 'COMMON.DERIV'
1162 include 'COMMON.GEO'
1163 include 'COMMON.LOCAL'
1164 include 'COMMON.INTERACT'
1165 include 'COMMON.IOUNITS'
1166 include 'COMMON.NAMES'
1169 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1173 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1177 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1180 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1186 c-----------------------------------------------------------------
1187 subroutine sddir_precalc
1188 c Applying velocity Verlet algorithm - step 1 to coordinates
1190 include 'DIMENSIONS'
1194 include 'COMMON.CONTROL'
1195 include 'COMMON.VAR'
1198 include 'COMMON.LAGRANGE.5diag'
1200 include 'COMMON.LAGRANGE'
1203 include 'COMMON.LANGEVIN'
1206 include 'COMMON.LANGEVIN.lang0.5diag'
1208 include 'COMMON.LANGEVIN.lang0'
1211 include 'COMMON.CHAIN'
1212 include 'COMMON.DERIV'
1213 include 'COMMON.GEO'
1214 include 'COMMON.LOCAL'
1215 include 'COMMON.INTERACT'
1216 include 'COMMON.IOUNITS'
1217 include 'COMMON.NAMES'
1218 include 'COMMON.TIME1'
1219 double precision time00
1220 double precision stochforcvec(MAXRES6)
1221 common /stochcalc/ stochforcvec
1224 c Compute friction and stochastic forces
1232 c write (iout,*) "After friction_force"
1234 time_fric=time_fric+MPI_Wtime()-time00
1237 time_fric=time_fric+tcpu()-time00
1240 call stochastic_force(stochforcvec)
1241 c write (iout,*) "After stochastic_force"
1243 time_stoch=time_stoch+MPI_Wtime()-time00
1245 time_stoch=time_stoch+tcpu()-time00
1248 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1249 c forces (d_as_work)
1252 c write (iout,*) "friction accelerations"
1253 call fivediaginv_mult(dimen,fric_work, d_af_work)
1254 c write (iout,*) "stochastic acceleratios"
1255 call fivediaginv_mult(dimen,stochforcvec, d_as_work)
1256 c write (iout,*) "Leaving sddir_precalc"
1258 call ginv_mult(fric_work, d_af_work)
1259 call ginv_mult(stochforcvec, d_as_work)
1262 write (iout,*) "d_af_work"
1263 write (iout,'(3f10.5)') (d_af_work(i),i=1,dimen3)
1264 write (iout,*) "d_as_work"
1265 write (iout,'(3f10.5)') (d_as_work(i),i=1,dimen3)
1269 c---------------------------------------------------------------------
1270 subroutine sddir_verlet1
1271 c Applying velocity Verlet algorithm - step 1 to velocities
1273 include 'DIMENSIONS'
1274 include 'COMMON.CONTROL'
1275 include 'COMMON.VAR'
1278 include 'COMMON.LAGRANGE.5diag'
1280 include 'COMMON.LAGRANGE'
1283 include 'COMMON.LANGEVIN'
1286 include 'COMMON.LANGEVIN.lang0.5diag'
1288 include 'COMMON.LANGEVIN.lang0'
1291 include 'COMMON.CHAIN'
1292 include 'COMMON.DERIV'
1293 include 'COMMON.GEO'
1294 include 'COMMON.LOCAL'
1295 include 'COMMON.INTERACT'
1296 include 'COMMON.IOUNITS'
1297 include 'COMMON.NAMES'
1298 c Revised 3/31/05 AL: correlation between random contributions to
1299 c position and velocity increments included.
1300 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1301 double precision adt,adt2
1302 integer i,j,ind,inres
1304 c Add the contribution from BOTH friction and stochastic force to the
1305 c coordinates, but ONLY the contribution from the friction forces to velocities
1308 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1309 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1310 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1311 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1312 d_t(j,0)=d_t_old(j,0)+adt
1317 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1318 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1319 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1320 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1321 d_t(j,i)=d_t_old(j,i)+adt
1326 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1329 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1330 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1331 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1332 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1333 d_t(j,inres)=d_t_old(j,inres)+adt
1340 c---------------------------------------------------------------------
1341 subroutine sddir_verlet2
1342 c Calculating the adjusted velocities for accelerations
1344 include 'DIMENSIONS'
1345 include 'COMMON.CONTROL'
1346 include 'COMMON.VAR'
1349 include 'COMMON.LAGRANGE.5diag'
1351 include 'COMMON.LAGRANGE'
1354 include 'COMMON.LANGEVIN'
1357 include 'COMMON.LANGEVIN.lang0.5diag'
1359 include 'COMMON.LANGEVIN.lang0'
1362 include 'COMMON.CHAIN'
1363 include 'COMMON.DERIV'
1364 include 'COMMON.GEO'
1365 include 'COMMON.LOCAL'
1366 include 'COMMON.INTERACT'
1367 include 'COMMON.IOUNITS'
1368 include 'COMMON.NAMES'
1369 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1370 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1371 integer i,j,inres,ind
1372 c Revised 3/31/05 AL: correlation between random contributions to
1373 c position and velocity increments included.
1374 c The correlation coefficients are calculated at low-friction limit.
1375 c Also, friction forces are now not calculated with new velocities.
1377 c call friction_force
1378 call stochastic_force(stochforcvec)
1380 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1381 c forces (d_as_work)
1384 call fivediaginv_mult(maxres6,stochforcvec, d_as_work1)
1386 call ginv_mult(stochforcvec, d_as_work1)
1392 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1393 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1398 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1399 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1404 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1407 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1408 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1409 & +cos60*d_as_work1(ind+j))*d_time
1416 c---------------------------------------------------------------------
1417 subroutine max_accel
1419 c Find the maximum difference in the accelerations of the the sites
1420 c at the beginning and the end of the time step.
1423 include 'DIMENSIONS'
1424 include 'COMMON.CONTROL'
1425 include 'COMMON.VAR'
1428 include 'COMMON.LAGRANGE.5diag'
1430 include 'COMMON.LAGRANGE'
1432 include 'COMMON.CHAIN'
1433 include 'COMMON.DERIV'
1434 include 'COMMON.GEO'
1435 include 'COMMON.LOCAL'
1436 include 'COMMON.INTERACT'
1437 include 'COMMON.IOUNITS'
1439 double precision aux(3),accel(3),accel_old(3),dacc
1441 c aux(j)=d_a(j,0)-d_a_old(j,0)
1442 accel_old(j)=d_a_old(j,0)
1449 c 7/3/08 changed to asymmetric difference
1451 c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1452 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1453 accel(j)=accel(j)+0.5d0*d_a(j,i)
1454 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1455 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1456 dacc=dabs(accel(j)-accel_old(j))
1457 c write (iout,*) i,dacc
1458 if (dacc.gt.amax) amax=dacc
1466 accel_old(j)=d_a_old(j,0)
1471 accel_old(j)=accel_old(j)+d_a_old(j,1)
1472 accel(j)=accel(j)+d_a(j,1)
1476 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1478 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1479 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1480 accel(j)=accel(j)+d_a(j,i+nres)
1484 c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1485 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1486 dacc=dabs(accel(j)-accel_old(j))
1487 c write (iout,*) "side-chain",i,dacc
1488 if (dacc.gt.amax) amax=dacc
1492 accel_old(j)=accel_old(j)+d_a_old(j,i)
1493 accel(j)=accel(j)+d_a(j,i)
1494 c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1499 c---------------------------------------------------------------------
1500 subroutine predict_edrift(epdrift)
1502 c Predict the drift of the potential energy
1505 include 'DIMENSIONS'
1506 include 'COMMON.CONTROL'
1507 include 'COMMON.VAR'
1510 include 'COMMON.LAGRANGE.5diag'
1512 include 'COMMON.LAGRANGE'
1514 include 'COMMON.CHAIN'
1515 include 'COMMON.DERIV'
1516 include 'COMMON.GEO'
1517 include 'COMMON.LOCAL'
1518 include 'COMMON.INTERACT'
1519 include 'COMMON.IOUNITS'
1520 include 'COMMON.MUCA'
1521 double precision epdrift,epdriftij
1523 c Drift of the potential energy
1529 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1530 if (lmuca) epdriftij=epdriftij*factor
1531 c write (iout,*) "back",i,j,epdriftij
1532 if (epdriftij.gt.epdrift) epdrift=epdriftij
1536 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1539 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1540 if (lmuca) epdriftij=epdriftij*factor
1541 c write (iout,*) "side",i,j,epdriftij
1542 if (epdriftij.gt.epdrift) epdrift=epdriftij
1546 epdrift=0.5d0*epdrift*d_time*d_time
1547 c write (iout,*) "epdrift",epdrift
1550 c-----------------------------------------------------------------------
1551 subroutine verlet_bath
1553 c Coupling to the thermostat by using the Berendsen algorithm
1556 include 'DIMENSIONS'
1557 include 'COMMON.CONTROL'
1558 include 'COMMON.VAR'
1561 include 'COMMON.LAGRANGE.5diag'
1563 include 'COMMON.LAGRANGE'
1567 include 'COMMON.LANGEVIN.lang0.5diag'
1569 include 'COMMON.LANGEVIN.lang0'
1572 include 'COMMON.LANGEVIN'
1574 include 'COMMON.CHAIN'
1575 include 'COMMON.DERIV'
1576 include 'COMMON.GEO'
1577 include 'COMMON.LOCAL'
1578 include 'COMMON.INTERACT'
1579 include 'COMMON.IOUNITS'
1580 include 'COMMON.NAMES'
1581 double precision T_half,fact
1583 T_half=2.0d0/(dimen3*Rb)*EK
1584 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1585 c write(iout,*) "T_half", T_half
1586 c write(iout,*) "EK", EK
1587 c write(iout,*) "fact", fact
1589 d_t(j,0)=fact*d_t(j,0)
1593 d_t(j,i)=fact*d_t(j,i)
1597 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1600 d_t(j,inres)=fact*d_t(j,inres)
1606 c---------------------------------------------------------
1608 c Set up the initial conditions of a MD simulation
1610 include 'DIMENSIONS'
1614 integer IERROR,ERRCODE,error_msg,ierr,ierrcode
1616 include 'COMMON.SETUP'
1617 include 'COMMON.CONTROL'
1618 include 'COMMON.VAR'
1621 include 'COMMON.LAGRANGE.5diag'
1623 include 'COMMON.LAGRANGE'
1625 include 'COMMON.QRESTR'
1627 include 'COMMON.LANGEVIN'
1630 include 'COMMON.LANGEVIN.lang0.5diag'
1633 include 'COMMON.LANGEVIN.lang0.5diag'
1635 include 'COMMON.LANGEVIN.lang0'
1639 include 'COMMON.CHAIN'
1640 include 'COMMON.DERIV'
1641 include 'COMMON.GEO'
1642 include 'COMMON.LOCAL'
1643 include 'COMMON.INTERACT'
1644 include 'COMMON.IOUNITS'
1645 include 'COMMON.NAMES'
1646 include 'COMMON.REMD'
1647 include 'COMMON.TIME1'
1648 double precision tt0
1649 real*8 energia_long(0:n_ene),
1650 & energia_short(0:n_ene),vcm(3),incr(3)
1651 double precision cm(3),L(3),xv,sigv,lowb,highb
1652 double precision varia(maxvar),energia(0:n_ene)
1659 integer i,ipos,iq,iw,j,iranmin,nft_sc,iretcode,nfun,itrial,itmp,
1662 double precision etot
1664 write (iout,*) "init_MD INDPDB",indpdb
1666 c write(iout,*) "d_time", d_time
1667 c Compute the standard deviations of stochastic forces for Langevin dynamics
1668 c if the friction coefficients do not depend on surface area
1669 if (lang.gt.0 .and. .not.surfarea) then
1671 stdforcp(i)=stdfp*dsqrt(gamp)
1674 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1675 & *dsqrt(gamsc(iabs(itype(i))))
1678 c Open the pdb file for snapshotshots
1681 if (ilen(tmpdir).gt.0)
1682 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1683 & liczba(:ilen(liczba))//".pdb")
1685 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1689 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1690 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1691 & liczba(:ilen(liczba))//".x")
1692 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1695 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1696 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1697 & liczba(:ilen(liczba))//".cx")
1698 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1704 if (ilen(tmpdir).gt.0)
1705 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1706 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1708 if (ilen(tmpdir).gt.0)
1709 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1710 cartname=prefix(:ilen(prefix))//"_MD.cx"
1714 write (qstr,'(256(1h ))')
1717 iq = qinfrag(i,iset)*10
1718 iw = wfrag(i,iset)/100
1720 if(me.eq.king.or..not.out1file)
1721 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1722 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1727 iq = qinpair(i,iset)*10
1728 iw = wpair(i,iset)/100
1730 if(me.eq.king.or..not.out1file)
1731 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1732 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1736 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1738 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1740 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1742 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1745 write (iout,*) "REST ",rest
1747 if (restart1file) then
1749 & inquire(file=mremd_rst_name,exist=file_exist)
1751 write (*,*) me," Before broadcast: file_exist",file_exist
1752 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1754 write (*,*) me," After broadcast: file_exist",file_exist
1755 c inquire(file=mremd_rst_name,exist=file_exist)
1757 if(me.eq.king.or..not.out1file)
1758 & write(iout,*) "Initial state read by master and distributed"
1760 if (ilen(tmpdir).gt.0)
1761 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1762 & //liczba(:ilen(liczba))//'.rst')
1763 inquire(file=rest2name,exist=file_exist)
1766 if(.not.restart1file) then
1767 if(me.eq.king.or..not.out1file)
1768 & write(iout,*) "Initial state will be read from file ",
1769 & rest2name(:ilen(rest2name))
1772 call rescale_weights(t_bath)
1775 if(me.eq.king.or..not.out1file)then
1776 if (restart1file) then
1777 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1780 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1783 write(iout,*) "Initial velocities randomly generated"
1790 c Generate initial velocities
1791 if(me.eq.king.or..not.out1file)
1792 & write(iout,*) "Initial velocities randomly generated"
1795 CtotTafm is the variable for AFM time which eclipsed during
1798 c rest2name = prefix(:ilen(prefix))//'.rst'
1799 if(me.eq.king.or..not.out1file)then
1800 write (iout,*) "Initial velocities"
1802 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1803 & (d_t(j,i+nres),j=1,3)
1805 c Zeroing the total angular momentum of the system
1806 write(iout,*) "Calling the zero-angular
1807 & momentum subroutine"
1810 c Getting the potential energy and forces and velocities and accelerations
1812 write (iout,*) "velocity of the center of the mass:"
1813 write (iout,*) (vcm(j),j=1,3)
1816 d_t(j,0)=d_t(j,0)-vcm(j)
1818 c Removing the velocity of the center of mass
1820 if(me.eq.king.or..not.out1file)then
1821 write (iout,*) "vcm right after adjustment:"
1822 write (iout,*) (vcm(j),j=1,3)
1825 write (iout,*) "init_MD before initial structure REST ",rest
1827 if (iranconf.ne.0) then
1828 c 8/22/17 AL Loop to produce a low-energy random conformation
1832 print *, 'Calling OVERLAP_SC'
1833 call overlap_sc(fail)
1837 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1838 print *,'SC_move',nft_sc,etot
1839 if(me.eq.king.or..not.out1file)
1840 & write(iout,*) 'SC_move',nft_sc,etot
1844 if (me.eq.king.or..not.out1file) write(iout,*)
1845 & 'Minimizing random structure: Calling MINIM_DC'
1846 call minim_dc(etot,iretcode,nfun)
1848 call geom_to_var(nvar,varia)
1849 if (me.eq.king.or..not.out1file) write (iout,*)
1850 & 'Minimizing random structure: Calling MINIMIZE.'
1851 call minimize(etot,varia,iretcode,nfun)
1852 call var_to_geom(nvar,varia)
1854 if (me.eq.king.or..not.out1file)
1855 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1856 if (isnan(etot) .or. etot.gt.1.0d4) then
1857 write (iout,*) "Energy too large",etot,
1858 & " trying another random conformation"
1861 call gen_rand_conf(itmp,*30)
1863 30 write (iout,*) 'Failed to generate random conformation',
1864 & ', itrial=',itrial
1865 write (*,*) 'Processor:',me,
1866 & ' Failed to generate random conformation',
1876 write (iout,'(a,i3,a)') 'Processor:',me,
1877 & ' error in generating random conformation.'
1878 write (*,'(a,i3,a)') 'Processor:',me,
1879 & ' error in generating random conformation.'
1882 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1891 write (iout,'(a,i3,a)') 'Processor:',me,
1892 & ' failed to generate a low-energy random conformation.'
1893 write (*,'(a,i3,a)') 'Processor:',me,
1894 & ' failed to generate a low-energy random conformation.'
1897 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1902 else if (preminim) then
1903 if (start_from_model) then
1904 i_model=iran_num(1,constr_homology)
1905 write (iout,*) 'starting from model ',i_model
1908 c(j,i)=chomo(j,i,i_model)
1911 call int_from_cart(.true.,.false.)
1912 call sc_loc_geom(.false.)
1915 dc(j,i)=c(j,i+1)-c(j,i)
1916 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1921 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1922 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1926 ! Remove SC overlaps if requested
1928 write (iout,*) 'Calling OVERLAP_SC'
1929 call overlap_sc(fail)
1931 ! Search for better SC rotamers if requested
1933 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1934 print *,'SC_move',nft_sc,etot
1935 if (me.eq.king.or..not.out1file)
1936 & write(iout,*) 'SC_move',nft_sc,etot
1938 call etotal(energia(0))
1939 C 8/22/17 AL Minimize initial structure
1941 if (me.eq.king.or..not.out1file) write(iout,*)
1942 & 'Minimizing initial PDB structure: Calling MINIM_DC'
1943 call minim_dc(etot,iretcode,nfun)
1945 call geom_to_var(nvar,varia)
1946 if(me.eq.king.or..not.out1file) write (iout,*)
1947 & 'Minimizing initial PDB structure: Calling MINIMIZE.'
1948 call minimize(etot,varia,iretcode,nfun)
1949 call var_to_geom(nvar,varia)
1951 if (me.eq.king.or..not.out1file)
1952 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1953 if(me.eq.king.or..not.out1file)
1954 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1957 call chainbuild_cart
1962 kinetic_T=2.0d0/(dimen3*Rb)*EK
1963 if(me.eq.king.or..not.out1file)then
1973 call etotal(potEcomp)
1974 if (large) call enerprint(potEcomp)
1977 t_etotal=t_etotal+MPI_Wtime()-tt0
1979 t_etotal=t_etotal+tcpu()-tt0
1983 c write (iout,*) "PotE-homology",potE-potEcomp(27)
1987 if (amax*d_time .gt. dvmax) then
1988 d_time=d_time*dvmax/amax
1989 if(me.eq.king.or..not.out1file) write (iout,*)
1990 & "Time step reduced to",d_time,
1991 & " because of too large initial acceleration."
1993 if(me.eq.king.or..not.out1file)then
1994 write(iout,*) "Potential energy and its components"
1995 call enerprint(potEcomp)
1996 c write(iout,*) (potEcomp(i),i=0,n_ene)
1998 potE=potEcomp(0)-potEcomp(27)
1999 c write (iout,*) "PotE-homology",potE
2002 if (ntwe.ne.0) call statout(itime)
2003 if(me.eq.king.or..not.out1file)
2004 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
2005 & " Kinetic energy",EK," Potential energy",potE,
2006 & " Total energy",totE," Maximum acceleration ",
2009 write (iout,*) "Initial coordinates"
2011 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
2012 & (c(j,i+nres),j=1,3)
2014 write (iout,*) "Initial dC"
2016 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
2017 & (dc(j,i+nres),j=1,3)
2019 write (iout,*) "Initial velocities"
2020 write (iout,"(13x,' backbone ',23x,' side chain')")
2022 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
2023 & (d_t(j,i+nres),j=1,3)
2025 write (iout,*) "Initial accelerations"
2027 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2028 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
2029 & (d_a(j,i+nres),j=1,3)
2035 d_t_old(j,i)=d_t(j,i)
2036 d_a_old(j,i)=d_a(j,i)
2038 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2047 call etotal_short(energia_short)
2048 if (large) call enerprint(potEcomp)
2051 t_eshort=t_eshort+MPI_Wtime()-tt0
2053 t_eshort=t_eshort+tcpu()-tt0
2058 if(.not.out1file .and. large) then
2059 write (iout,*) "energia_long",energia_long(0),
2060 & " energia_short",energia_short(0),
2061 & " total",energia_long(0)+energia_short(0)
2062 write (iout,*) "Initial fast-force accelerations"
2064 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2065 & (d_a(j,i+nres),j=1,3)
2068 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2071 d_a_short(j,i)=d_a(j,i)
2080 call etotal_long(energia_long)
2081 if (large) call enerprint(potEcomp)
2084 t_elong=t_elong+MPI_Wtime()-tt0
2086 t_elong=t_elong+tcpu()-tt0
2091 if(.not.out1file .and. large) then
2092 write (iout,*) "energia_long",energia_long(0)
2093 write (iout,*) "Initial slow-force accelerations"
2095 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2096 & (d_a(j,i+nres),j=1,3)
2100 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2102 t_enegrad=t_enegrad+tcpu()-tt0
2107 c-----------------------------------------------------------
2108 subroutine random_vel
2110 include 'DIMENSIONS'
2111 include 'COMMON.CONTROL'
2112 include 'COMMON.VAR'
2115 include 'COMMON.LAGRANGE.5diag'
2117 include 'COMMON.LAGRANGE'
2120 include 'COMMON.LANGEVIN'
2123 include 'COMMON.LANGEVIN.lang0.5diag'
2125 include 'COMMON.LANGEVIN.lang0'
2128 include 'COMMON.CHAIN'
2129 include 'COMMON.DERIV'
2130 include 'COMMON.GEO'
2131 include 'COMMON.LOCAL'
2132 include 'COMMON.INTERACT'
2133 include 'COMMON.IOUNITS'
2134 include 'COMMON.NAMES'
2135 include 'COMMON.TIME1'
2136 double precision xv,sigv,lowb,highb,vec_afm(3),Ek1,Ek2,Ek3,aux
2137 integer i,ii,j,k,l,ind
2138 double precision anorm_distr
2139 logical lprn /.true./
2141 integer ichain,n,innt,inct,ibeg,ierr
2142 double precision work(8*maxres6)
2143 integer iwork(maxres6)
2144 double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),
2145 & Gvec(maxres2_chain,maxres2_chain)
2146 common /przechowalnia/Ghalf,Geigen,Gvec
2148 double precision inertia(maxres2_chain,maxres2_chain)
2150 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2151 c First generate velocities in the eigenspace of the G matrix
2152 c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2155 write (iout,*) "Random_vel, fivediag"
2164 n=dimen_chain(ichain)
2165 innt=iposd_chain(ichain)
2168 write (iout,*) "Chain",ichain," n",n," start",innt
2170 if (i.lt.inct-1) then
2171 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),
2173 else if (i.eq.inct-1) then
2174 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
2176 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
2180 ghalf(ind+1)=dmorig(innt)
2181 ghalf(ind+2)=du1orig(innt)
2182 ghalf(ind+3)=dmorig(innt+1)
2186 c write (iout,*) "i",i," ind",ind," indu2",innt+i-2,
2187 c & " indu1",innt+i-1," indm",innt+i
2188 ghalf(ind+1)=du2orig(innt-1+i-2)
2189 ghalf(ind+2)=du1orig(innt-1+i-1)
2190 ghalf(ind+3)=dmorig(innt-1+i)
2191 c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
2192 c & "DU1",innt-1+i-1,"DM ",innt-1+i
2200 inertia(i,j)=ghalf(ind)
2201 inertia(j,i)=ghalf(ind)
2206 write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
2207 write (iout,*) "Five-diagonal inertia matrix, lower triangle"
2208 call matoutr(n,ghalf)
2210 call gldiag(maxres2_chain,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
2212 write (iout,'(//a,i3)')
2213 & "Eigenvectors and eigenvalues of the G matrix chain",ichain
2214 call eigout(n,n,maxres2_chain,maxres2_chain,Gvec,Geigen)
2217 c check diagonalization
2223 aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
2227 write (iout,*) i,j,aux,geigen(i)
2229 write (iout,*) i,j,aux
2239 sigv=dsqrt((Rb*t_bath)/geigen(i))
2242 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2243 EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
2244 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2245 c & " d_t_work_new",d_t_work_new(ii)
2253 d_t_work(ind)=d_t_work(ind)
2254 & +Gvec(i,j)*d_t_work_new((j-1)*3+k)
2256 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2265 aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
2271 c Transfer to the d_t vector
2272 innt=chain_border(1,ichain)
2273 inct=chain_border(2,ichain)
2275 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
2279 d_t(j,i)=d_t_work(ind)
2281 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2284 d_t(j,i+nres)=d_t_work(ind)
2291 write (iout,*) "Random velocities in the Calpha,SC space"
2293 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2294 & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2297 call kinetic_CASC(Ek1)
2299 ! Transform the velocities to virtual-bond space
2304 if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
2306 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2310 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2311 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2317 c write (iout,*) "Shifting accelerations"
2319 c write (iout,*) "ichain",chain_border1(1,ichain)-1,
2320 c & chain_border1(1,ichain)
2321 d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
2322 d_t(:,chain_border1(1,ichain))=0.0d0
2324 c write (iout,*) "Adding accelerations"
2326 c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2327 c & chain_border(2,ichain-1)
2328 d_t(:,chain_border1(1,ichain)-1)=
2329 & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2330 d_t(:,chain_border(2,ichain-1))=0.0d0
2335 c d_t(j,0)=d_t(j,nnt)
2338 innt=chain_border(1,ichain)
2339 inct=chain_border(2,ichain)
2340 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
2341 c write (iout,*) "ibeg",ibeg
2343 d_t(j,ibeg)=d_t(j,innt)
2347 if (iabs(itype(i).eq.10)) then
2348 c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
2350 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2354 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2355 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2364 & "Random velocities in the virtual-bond-vector space"
2365 write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
2367 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2368 & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2371 write (iout,*) "Kinetic energy from inertia matrix eigenvalues",
2374 & "Kinetic temperatures from inertia matrix eigenvalues",
2377 write (iout,*) "Kinetic energy from inertia matrix",Ek3
2378 write (iout,*) "Kinetic temperatures from inertia",
2379 & 2*Ek3/(3*dimen*Rb)
2381 write (iout,*) "Kinetic energy from velocities in CA-SC space",
2384 & "Kinetic temperatures from velovities in CA-SC space",
2385 & 2*Ek1/(3*dimen*Rb)
2388 & "Kinetic energy from virtual-bond-vector velocities",Ek1
2390 & "Kinetic temperature from virtual-bond-vector velocities ",
2399 sigv=dsqrt((Rb*t_bath)/geigen(i))
2402 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2404 c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2405 c & " d_t_work_new",d_t_work_new(ii)
2408 C if (SELFGUIDE.gt.0) then
2411 C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
2412 C distance=distance+vec_afm(j)**2
2414 C distance=dsqrt(distance)
2416 C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
2417 C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
2418 C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
2419 C & d_t_work_new(j+(afmend-1)*3)
2430 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2433 c write (iout,*) "Ek from eigenvectors",Ek1
2435 c Transform velocities to UNRES coordinate space
2441 d_t_work(ind)=d_t_work(ind)
2442 & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2444 c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2448 c Transfer to the d_t vector
2450 d_t(j,0)=d_t_work(j)
2456 d_t(j,i)=d_t_work(ind)
2460 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2463 d_t(j,i+nres)=d_t_work(ind)
2468 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2469 c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2475 c-----------------------------------------------------------
2476 subroutine sd_verlet_p_setup
2477 c Sets up the parameters of stochastic Verlet algorithm
2479 include 'DIMENSIONS'
2483 include 'COMMON.CONTROL'
2484 include 'COMMON.VAR'
2487 include 'COMMON.LANGEVIN'
2490 include 'COMMON.LANGEVIN.lang0.5diag'
2492 include 'COMMON.LANGEVIN.lang0'
2495 include 'COMMON.CHAIN'
2496 include 'COMMON.DERIV'
2497 include 'COMMON.GEO'
2498 include 'COMMON.LOCAL'
2499 include 'COMMON.INTERACT'
2500 include 'COMMON.IOUNITS'
2501 include 'COMMON.NAMES'
2502 include 'COMMON.TIME1'
2503 double precision emgdt(MAXRES6),
2504 & pterm,vterm,rho,rhoc,vsig,
2505 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2506 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2507 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2508 logical lprn /.false./
2509 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2510 double precision ktm
2517 c AL 8/17/04 Code adapted from tinker
2519 c Get the frictional and random terms for stochastic dynamics in the
2520 c eigenspace of mass-scaled UNRES friction matrix
2523 gdt = fricgam(i) * d_time
2525 c Stochastic dynamics reduces to simple MD for zero friction
2527 if (gdt .le. zero) then
2528 pfric_vec(i) = 1.0d0
2529 vfric_vec(i) = d_time
2530 afric_vec(i) = 0.5d0 * d_time * d_time
2531 prand_vec(i) = 0.0d0
2532 vrand_vec1(i) = 0.0d0
2533 vrand_vec2(i) = 0.0d0
2535 c Analytical expressions when friction coefficient is large
2538 if (gdt .ge. gdt_radius) then
2541 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2542 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2543 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2544 vterm = 1.0d0 - egdt**2
2545 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2547 c Use series expansions when friction coefficient is small
2558 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2559 & - gdt5/120.0d0 + gdt6/720.0d0
2560 & - gdt7/5040.0d0 + gdt8/40320.0d0
2561 & - gdt9/362880.0d0) / fricgam(i)**2
2562 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2563 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2564 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2565 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2566 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2567 & + 127.0d0*gdt9/90720.0d0
2568 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2569 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2570 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2571 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2572 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2573 & - 17.0d0*gdt2/1280.0d0
2574 & + 17.0d0*gdt3/6144.0d0
2575 & + 40967.0d0*gdt4/34406400.0d0
2576 & - 57203.0d0*gdt5/275251200.0d0
2577 & - 1429487.0d0*gdt6/13212057600.0d0)
2580 c Compute the scaling factors of random terms for the nonzero friction case
2582 ktm = 0.5d0*d_time/fricgam(i)
2583 psig = dsqrt(ktm*pterm) / fricgam(i)
2584 vsig = dsqrt(ktm*vterm)
2585 rhoc = dsqrt(1.0d0 - rho*rho)
2587 vrand_vec1(i) = vsig * rho
2588 vrand_vec2(i) = vsig * rhoc
2593 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2596 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2597 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2601 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2603 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2604 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2605 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2606 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2607 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2608 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2610 t_sdsetup=t_sdsetup+MPI_Wtime()
2612 t_sdsetup=t_sdsetup+tcpu()-tt0
2616 c-------------------------------------------------------------
2617 subroutine eigtransf1(n,ndim,ab,d,c)
2620 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2626 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2632 c-------------------------------------------------------------
2633 subroutine eigtransf(n,ndim,a,b,d,c)
2636 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2642 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2648 c-------------------------------------------------------------
2649 subroutine sd_verlet1
2650 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2652 include 'DIMENSIONS'
2653 include 'COMMON.CONTROL'
2654 include 'COMMON.VAR'
2657 include 'COMMON.LAGRANGE.5diag'
2659 include 'COMMON.LAGRANGE'
2662 include 'COMMON.LANGEVIN'
2665 include 'COMMON.LANGEVIN.lang0.5diag'
2667 include 'COMMON.LANGEVIN.lang0'
2670 include 'COMMON.CHAIN'
2671 include 'COMMON.DERIV'
2672 include 'COMMON.GEO'
2673 include 'COMMON.LOCAL'
2674 include 'COMMON.INTERACT'
2675 include 'COMMON.IOUNITS'
2676 include 'COMMON.NAMES'
2677 double precision stochforcvec(MAXRES6)
2678 common /stochcalc/ stochforcvec
2679 logical lprn /.false./
2680 integer i,j,ind,inres
2682 c write (iout,*) "dc_old"
2684 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2685 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2688 dc_work(j)=dc_old(j,0)
2689 d_t_work(j)=d_t_old(j,0)
2690 d_a_work(j)=d_a_old(j,0)
2695 dc_work(ind+j)=dc_old(j,i)
2696 d_t_work(ind+j)=d_t_old(j,i)
2697 d_a_work(ind+j)=d_a_old(j,i)
2702 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2704 dc_work(ind+j)=dc_old(j,i+nres)
2705 d_t_work(ind+j)=d_t_old(j,i+nres)
2706 d_a_work(ind+j)=d_a_old(j,i+nres)
2714 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2718 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2719 & vfric_mat(i,j),afric_mat(i,j),
2720 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2728 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2729 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2730 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2731 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2733 d_t_work_new(i)=ddt1+0.5d0*ddt2
2734 d_t_work(i)=ddt1+ddt2
2739 d_t(j,0)=d_t_work(j)
2744 dc(j,i)=dc_work(ind+j)
2745 d_t(j,i)=d_t_work(ind+j)
2750 if (itype(i).ne.10) then
2753 dc(j,inres)=dc_work(ind+j)
2754 d_t(j,inres)=d_t_work(ind+j)
2761 c--------------------------------------------------------------------------
2762 subroutine sd_verlet2
2763 c Calculating the adjusted velocities for accelerations
2765 include 'DIMENSIONS'
2766 include 'COMMON.CONTROL'
2767 include 'COMMON.VAR'
2770 include 'COMMON.LAGRANGE.5diag'
2772 include 'COMMON.LAGRANGE'
2775 include 'COMMON.LANGEVIN'
2778 include 'COMMON.LANGEVIN.lang0.5diag'
2780 include 'COMMON.LANGEVIN.lang0'
2783 include 'COMMON.CHAIN'
2784 include 'COMMON.DERIV'
2785 include 'COMMON.GEO'
2786 include 'COMMON.LOCAL'
2787 include 'COMMON.INTERACT'
2788 include 'COMMON.IOUNITS'
2789 include 'COMMON.NAMES'
2790 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2791 common /stochcalc/ stochforcvec
2792 integer i,j,ind,inres
2794 c Compute the stochastic forces which contribute to velocity change
2796 call stochastic_force(stochforcvecV)
2803 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2804 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2805 & vrand_mat2(i,j)*stochforcvecV(j)
2807 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2811 d_t(j,0)=d_t_work(j)
2816 d_t(j,i)=d_t_work(ind+j)
2821 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2824 d_t(j,inres)=d_t_work(ind+j)
2831 c-----------------------------------------------------------
2832 subroutine sd_verlet_ciccotti_setup
2833 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2836 include 'DIMENSIONS'
2840 include 'COMMON.CONTROL'
2841 include 'COMMON.VAR'
2844 include 'COMMON.LAGRANGE.5diag'
2846 include 'COMMON.LAGRANGE'
2848 include 'COMMON.LANGEVIN'
2849 include 'COMMON.CHAIN'
2850 include 'COMMON.DERIV'
2851 include 'COMMON.GEO'
2852 include 'COMMON.LOCAL'
2853 include 'COMMON.INTERACT'
2854 include 'COMMON.IOUNITS'
2855 include 'COMMON.NAMES'
2856 include 'COMMON.TIME1'
2857 double precision emgdt(MAXRES6),
2858 & pterm,vterm,rho,rhoc,vsig,
2859 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2860 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2861 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2862 logical lprn /.false./
2863 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2864 double precision ktm
2872 c AL 8/17/04 Code adapted from tinker
2874 c Get the frictional and random terms for stochastic dynamics in the
2875 c eigenspace of mass-scaled UNRES friction matrix
2878 write (iout,*) "i",i," fricgam",fricgam(i)
2879 gdt = fricgam(i) * d_time
2881 c Stochastic dynamics reduces to simple MD for zero friction
2883 if (gdt .le. zero) then
2884 pfric_vec(i) = 1.0d0
2885 vfric_vec(i) = d_time
2886 afric_vec(i) = 0.5d0*d_time*d_time
2887 prand_vec(i) = afric_vec(i)
2888 vrand_vec2(i) = vfric_vec(i)
2890 c Analytical expressions when friction coefficient is large
2895 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2896 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2897 prand_vec(i) = afric_vec(i)
2898 vrand_vec2(i) = vfric_vec(i)
2900 c Compute the scaling factors of random terms for the nonzero friction case
2902 c ktm = 0.5d0*d_time/fricgam(i)
2903 c psig = dsqrt(ktm*pterm) / fricgam(i)
2904 c vsig = dsqrt(ktm*vterm)
2905 c prand_vec(i) = psig*afric_vec(i)
2906 c vrand_vec2(i) = vsig*vfric_vec(i)
2911 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2914 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2915 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2919 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2921 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2922 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2923 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2924 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2925 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2927 t_sdsetup=t_sdsetup+MPI_Wtime()
2929 t_sdsetup=t_sdsetup+tcpu()-tt0
2933 c-------------------------------------------------------------
2934 subroutine sd_verlet1_ciccotti
2935 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2937 include 'DIMENSIONS'
2941 include 'COMMON.CONTROL'
2942 include 'COMMON.VAR'
2945 include 'COMMON.LAGRANGE.5diag'
2947 include 'COMMON.LAGRANGE'
2949 include 'COMMON.LANGEVIN'
2950 include 'COMMON.CHAIN'
2951 include 'COMMON.DERIV'
2952 include 'COMMON.GEO'
2953 include 'COMMON.LOCAL'
2954 include 'COMMON.INTERACT'
2955 include 'COMMON.IOUNITS'
2956 include 'COMMON.NAMES'
2957 double precision stochforcvec(MAXRES6)
2958 common /stochcalc/ stochforcvec
2959 logical lprn /.false./
2962 c write (iout,*) "dc_old"
2964 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2965 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2968 dc_work(j)=dc_old(j,0)
2969 d_t_work(j)=d_t_old(j,0)
2970 d_a_work(j)=d_a_old(j,0)
2975 dc_work(ind+j)=dc_old(j,i)
2976 d_t_work(ind+j)=d_t_old(j,i)
2977 d_a_work(ind+j)=d_a_old(j,i)
2982 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2984 dc_work(ind+j)=dc_old(j,i+nres)
2985 d_t_work(ind+j)=d_t_old(j,i+nres)
2986 d_a_work(ind+j)=d_a_old(j,i+nres)
2995 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2999 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
3000 & vfric_mat(i,j),afric_mat(i,j),
3001 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3009 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
3010 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3011 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3012 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3014 d_t_work_new(i)=ddt1+0.5d0*ddt2
3015 d_t_work(i)=ddt1+ddt2
3020 d_t(j,0)=d_t_work(j)
3025 dc(j,i)=dc_work(ind+j)
3026 d_t(j,i)=d_t_work(ind+j)
3031 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3034 dc(j,inres)=dc_work(ind+j)
3035 d_t(j,inres)=d_t_work(ind+j)
3042 c--------------------------------------------------------------------------
3043 subroutine sd_verlet2_ciccotti
3044 c Calculating the adjusted velocities for accelerations
3046 include 'DIMENSIONS'
3047 include 'COMMON.CONTROL'
3048 include 'COMMON.VAR'
3051 include 'COMMON.LAGRANGE.5diag'
3053 include 'COMMON.LAGRANGE'
3055 include 'COMMON.LANGEVIN'
3056 include 'COMMON.CHAIN'
3057 include 'COMMON.DERIV'
3058 include 'COMMON.GEO'
3059 include 'COMMON.LOCAL'
3060 include 'COMMON.INTERACT'
3061 include 'COMMON.IOUNITS'
3062 include 'COMMON.NAMES'
3063 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
3064 common /stochcalc/ stochforcvec
3067 c Compute the stochastic forces which contribute to velocity change
3069 call stochastic_force(stochforcvecV)
3076 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3077 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3078 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3080 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3084 d_t(j,0)=d_t_work(j)
3089 d_t(j,i)=d_t_work(ind+j)
3094 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3097 d_t(j,inres)=d_t_work(ind+j)