2 c------------------------------------------------
3 c The driver for molecular dynamics subroutines
4 c------------------------------------------------
5 implicit real*8 (a-h,o-z)
10 include 'COMMON.CONTROL'
14 include 'COMMON.LANGEVIN'
16 include 'COMMON.LANGEVIN.lang0'
18 include 'COMMON.CHAIN'
19 include 'COMMON.DERIV'
21 include 'COMMON.LOCAL'
22 include 'COMMON.INTERACT'
23 include 'COMMON.IOUNITS'
24 include 'COMMON.NAMES'
25 include 'COMMON.TIME1'
26 include 'COMMON.HAIRPIN'
27 double precision cm(3),L(3),vcm(3)
29 double precision v_work(maxres6),v_transf(maxres6)
38 if (ilen(tmpdir).gt.0)
39 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
40 & //liczba(:ilen(liczba))//'.rst')
42 if (ilen(tmpdir).gt.0)
43 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
50 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
56 c Determine the inverse of the inertia matrix.
57 call setup_MD_matrices
61 t_MDsetup = MPI_Wtime()-tt0
63 t_MDsetup = tcpu()-tt0
66 c Entering the MD loop
72 if (lang.eq.2 .or. lang.eq.3) then
75 call sd_verlet_p_setup
77 call sd_verlet_ciccotti_setup
82 pfric0_mat(i,j,0)=pfric_mat(i,j)
83 afric0_mat(i,j,0)=afric_mat(i,j)
84 vfric0_mat(i,j,0)=vfric_mat(i,j)
85 prand0_mat(i,j,0)=prand_mat(i,j)
86 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
87 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
95 else if (lang.eq.1 .or. lang.eq.4) then
99 t_langsetup=MPI_Wtime()-tt0
102 t_langsetup=tcpu()-tt0
105 do itime=1,n_timestep
106 if (ovrtim()) goto 100
108 if (lang.gt.0 .and. surfarea .and.
109 & mod(itime,reset_fricmat).eq.0) then
110 if (lang.eq.2 .or. lang.eq.3) then
113 call sd_verlet_p_setup
115 call sd_verlet_ciccotti_setup
120 pfric0_mat(i,j,0)=pfric_mat(i,j)
121 afric0_mat(i,j,0)=afric_mat(i,j)
122 vfric0_mat(i,j,0)=vfric_mat(i,j)
123 prand0_mat(i,j,0)=prand_mat(i,j)
124 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
125 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
131 flag_stoch(i)=.false.
133 else if (lang.eq.1 .or. lang.eq.4) then
136 write (iout,'(a,i10)')
137 & "Friction matrix reset based on surface area, itime",itime
139 if (reset_vel .and. tbf .and. lang.eq.0
140 & .and. mod(itime,count_reset_vel).eq.0) then
142 write(iout,'(a,f20.2)')
143 & "Velocities reset to random values, time",totT
146 d_t_old(j,i)=d_t(j,i)
150 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
154 d_t(j,0)=d_t(j,0)-vcm(j)
157 kinetic_T=2.0d0/(dimen*Rb)*EK
158 scalfac=dsqrt(T_bath/kinetic_T)
159 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
162 d_t_old(j,i)=scalfac*d_t(j,i)
168 c Time-reversible RESPA algorithm
169 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
170 call RESPA_step(itime)
172 c Variable time step algorithm.
173 call velverlet_step(itime)
176 call brown_step(itime)
179 if (mod(itime,ntwe).eq.0) then
195 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
198 v_work(ind)=d_t(j,i+nres)
203 write (66,'(80f10.5)')
204 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
208 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
210 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
212 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
215 if (mod(itime,ntwx).eq.0) then
216 write (tytul,'("time",f8.2)') totT
218 call hairpin(.true.,nharp,iharp)
219 call secondary2(.true.)
220 call pdbout(potE,tytul,ipdb)
225 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
226 open(irest2,file=rest2name,status='unknown')
227 write(irest2,*) totT,EK,potE,totE,t_bath
229 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
232 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
244 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
246 & 'MD calculations setup:',t_MDsetup,
247 & 'Energy & gradient evaluation:',t_enegrad,
248 & 'Stochastic MD setup:',t_langsetup,
249 & 'Stochastic MD step setup:',t_sdsetup,
251 write (iout,'(/28(1h=),a25,27(1h=))')
252 & ' End of MD calculation '
255 c-------------------------------------------------------------------------------
256 subroutine brown_step(itime)
257 c------------------------------------------------
258 c Perform a single Euler integration step of Brownian dynamics
259 c------------------------------------------------
260 implicit real*8 (a-h,o-z)
265 include 'COMMON.CONTROL'
269 include 'COMMON.LANGEVIN'
271 include 'COMMON.LANGEVIN.lang0'
273 include 'COMMON.CHAIN'
274 include 'COMMON.DERIV'
276 include 'COMMON.LOCAL'
277 include 'COMMON.INTERACT'
278 include 'COMMON.IOUNITS'
279 include 'COMMON.NAMES'
280 include 'COMMON.TIME1'
281 double precision zapas(MAXRES6)
282 integer ilen,rstcount
284 double precision stochforcvec(MAXRES6)
285 double precision Bmat(MAXRES6,MAXRES2),Cmat(maxres2,maxres2),
286 & Cinv(maxres2,maxres2),GBmat(MAXRES6,MAXRES2),
287 & Tmat(MAXRES6,MAXRES2),Pmat(maxres6,maxres6),Td(maxres6),
289 common /stochcalc/ stochforcvec
292 logical lprn /.false./,lprn1 /.false./
294 double precision difftol /1.0d-5/
297 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) nbond=nbond+1
301 write (iout,*) "Generalized inverse of fricmat"
302 call matout(dimen,dimen,MAXRES6,MAXRES6,fricmat)
314 Bmat(ind+j,ind1)=dC_norm(j,i)
319 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
322 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
328 write (iout,*) "Matrix Bmat"
329 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Bmat)
335 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
340 write (iout,*) "Matrix GBmat"
341 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Gbmat)
347 Cmat(i,j)=Cmat(i,j)+Bmat(k,i)*GBmat(k,j)
352 write (iout,*) "Matrix Cmat"
353 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
355 call matinvert(nbond,MAXRES2,Cmat,Cinv)
357 write (iout,*) "Matrix Cinv"
358 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cinv)
364 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
369 write (iout,*) "Matrix Tmat"
370 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Tmat)
380 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
385 write (iout,*) "Matrix Pmat"
386 call MATOUT(dimen,dimen,MAXRES6,MAXRES6,Pmat)
393 Td(i)=Td(i)+vbl*Tmat(i,ind)
396 if (itype(k).ne.10 .and. itype(i).ne.ntyp1) then
398 Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
403 write (iout,*) "Vector Td"
405 write (iout,'(i5,f10.5)') i,Td(i)
408 call stochastic_force(stochforcvec)
410 write (iout,*) "stochforcvec"
412 write (iout,*) i,stochforcvec(i)
416 zapas(j)=-gcart(j,0)+stochforcvec(j)
418 dC_work(j)=dC_old(j,0)
424 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
425 dC_work(ind)=dC_old(j,i)
429 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
432 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
433 dC_work(ind)=dC_old(j,i+nres)
439 write (iout,*) "Initial d_t_work"
441 write (iout,*) i,d_t_work(i)
448 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
455 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
459 write (iout,*) "Final d_t_work and zapas"
461 write (iout,*) i,d_t_work(i),zapas(i)
475 dc_work(ind+j)=dc(j,i)
481 d_t(j,i+nres)=d_t_work(ind+j)
482 dc(j,i+nres)=zapas(ind+j)
483 dc_work(ind+j)=dc(j,i+nres)
489 write (iout,*) "Before correction for rotational lengthening"
490 write (iout,*) "New coordinates",
491 & " and differences between actual and standard bond lengths"
496 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
497 & i,(dC(j,i),j=1,3),xx
500 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
502 xx=vbld(i+nres)-vbldsc0(1,itype(i))
503 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
504 & i,(dC(j,i+nres),j=1,3),xx
508 c Second correction (rotational lengthening)
514 blen2 = scalar(dc(1,i),dc(1,i))
515 ppvec(ind)=2*vbl**2-blen2
516 diffbond=dabs(vbl-dsqrt(blen2))
517 if (diffbond.gt.diffmax) diffmax=diffbond
518 if (ppvec(ind).gt.0.0d0) then
519 ppvec(ind)=dsqrt(ppvec(ind))
524 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
528 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
530 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
531 ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
532 diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2))
533 if (diffbond.gt.diffmax) diffmax=diffbond
534 if (ppvec(ind).gt.0.0d0) then
535 ppvec(ind)=dsqrt(ppvec(ind))
540 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
544 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
545 if (diffmax.lt.difftol) goto 10
549 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
555 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
566 dc_work(ind+j)=zapas(ind+j)
571 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
573 dc(j,i+nres)=zapas(ind+j)
574 dc_work(ind+j)=zapas(ind+j)
579 c Building the chain from the newly calculated coordinates
582 if (large.and. mod(itime,ntwe).eq.0) then
583 write (iout,*) "Cartesian and internal coordinates: step 1"
586 write (iout,'(a)') "Potential forces"
588 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),
589 & (-gxcart(j,i),j=1,3)
591 write (iout,'(a)') "Stochastic forces"
593 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),
594 & (stochforc(j,i+nres),j=1,3)
596 write (iout,'(a)') "Velocities"
598 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
599 & (d_t(j,i+nres),j=1,3)
604 write (iout,*) "After correction for rotational lengthening"
605 write (iout,*) "New coordinates",
606 & " and differences between actual and standard bond lengths"
611 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
612 & i,(dC(j,i),j=1,3),xx
615 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
617 xx=vbld(i+nres)-vbldsc0(1,itype(i))
618 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
619 & i,(dC(j,i+nres),j=1,3),xx
624 c write (iout,*) "Too many attempts at correcting the bonds"
632 c Calculate energy and forces
634 call etotal(potEcomp)
635 potE=potEcomp(0)-potEcomp(20)
638 c Calculate the kinetic and total energy and the kinetic temperature
641 t_enegrad=t_enegrad+MPI_Wtime()-tt0
643 t_enegrad=t_enegrad+tcpu()-tt0
646 kinetic_T=2.0d0/(dimen*Rb)*EK
649 c-------------------------------------------------------------------------------
650 subroutine velverlet_step(itime)
651 c-------------------------------------------------------------------------------
652 c Perform a single velocity Verlet step; the time step can be rescaled if
653 c increments in accelerations exceed the threshold
654 c-------------------------------------------------------------------------------
655 implicit real*8 (a-h,o-z)
659 integer ierror,ierrcode
661 include 'COMMON.CONTROL'
665 include 'COMMON.LANGEVIN'
667 include 'COMMON.LANGEVIN.lang0'
669 include 'COMMON.CHAIN'
670 include 'COMMON.DERIV'
672 include 'COMMON.LOCAL'
673 include 'COMMON.INTERACT'
674 include 'COMMON.IOUNITS'
675 include 'COMMON.NAMES'
676 include 'COMMON.TIME1'
677 include 'COMMON.MUCA'
678 double precision vcm(3),incr(3)
679 double precision cm(3),L(3)
680 integer ilen,count,rstcount
683 integer maxcount_scale /20/
685 double precision stochforcvec(MAXRES6)
686 common /stochcalc/ stochforcvec
694 else if (lang.eq.2 .or. lang.eq.3) then
695 call stochastic_force(stochforcvec)
699 icount_scale=icount_scale+1
700 if (icount_scale.gt.maxcount_scale) then
702 & "ERROR: too many attempts at scaling down the time step. ",
703 & "amax=",amax,"epdrift=",epdrift,
704 & "damax=",damax,"edriftmax=",edriftmax,
708 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
712 c First step of the velocity Verlet algorithm
715 else if (lang.eq.3) then
716 call sd_verlet1_ciccotti
717 else if (lang.eq.1) then
722 c Build the chain from the newly calculated coordinates
724 if (rattle) call rattle1
726 if (large.and. mod(itime,ntwe).eq.0) then
727 write (iout,*) "Cartesian and internal coordinates: step 1"
732 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
733 & (dc(j,i+nres),j=1,3)
735 write (iout,*) "Accelerations"
737 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
738 & (d_a(j,i+nres),j=1,3)
740 write (iout,*) "Velocities, step 1"
742 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
743 & (d_t(j,i+nres),j=1,3)
752 c Calculate energy and forces
754 call etotal(potEcomp)
755 potE=potEcomp(0)-potEcomp(20)
757 c Get the new accelerations
760 t_enegrad=t_enegrad+MPI_Wtime()-tt0
762 t_enegrad=t_enegrad+tcpu()-tt0
764 c Determine maximum acceleration and scale down the timestep if needed
766 call predict_edrift(epdrift)
767 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
768 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
770 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
772 itime_scal=itime_scal+ifac_time
773 c fac_time=dmin1(damax/amax,0.5d0)
774 fac_time=0.5d0**ifac_time
775 d_time=d_time*fac_time
776 if (lang.eq.2 .or. lang.eq.3) then
777 c write (iout,*) "Calling sd_verlet_setup: 1"
778 c Rescale the stochastic forces and recalculate or restore
779 c the matrices of tinker integrator
780 if (itime_scal.gt.maxflag_stoch) then
781 if (large) write (iout,'(a,i5,a)')
782 & "Calculate matrices for stochastic step;",
783 & " itime_scal ",itime_scal
785 call sd_verlet_p_setup
787 call sd_verlet_ciccotti_setup
789 write (iout,'(2a,i3,a,i3,1h.)')
790 & "Warning: cannot store matrices for stochastic",
791 & " integration because the index",itime_scal,
792 & " is greater than",maxflag_stoch
793 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
794 & " integration Langevin algorithm for better efficiency."
795 else if (flag_stoch(itime_scal)) then
796 if (large) write (iout,'(a,i5,a,l1)')
797 & "Restore matrices for stochastic step; itime_scal ",
798 & itime_scal," flag ",flag_stoch(itime_scal)
802 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
803 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
804 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
805 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
806 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
807 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
812 if (large) write (iout,'(2a,i5,a,l1)')
813 & "Calculate & store matrices for stochastic step;",
814 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
816 call sd_verlet_p_setup
818 call sd_verlet_ciccotti_setup
820 flag_stoch(ifac_time)=.true.
824 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
825 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
826 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
827 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
828 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
829 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
834 fac_time=1.0d0/dsqrt(fac_time)
836 stochforcvec(i)=fac_time*stochforcvec(i)
838 else if (lang.eq.1) then
839 c Rescale the accelerations due to stochastic forces
840 fac_time=1.0d0/dsqrt(fac_time)
842 d_as_work(i)=d_as_work(i)*fac_time
845 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
846 & "itime",itime," Timestep scaled down to ",
847 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
849 c Second step of the velocity Verlet algorithm
852 else if (lang.eq.3) then
853 call sd_verlet2_ciccotti
854 else if (lang.eq.1) then
859 if (rattle) call rattle2
861 if (d_time.ne.d_time0) then
863 if (lang.eq.2 .or. lang.eq.3) then
864 if (large) write (iout,'(a)')
865 & "Restore original matrices for stochastic step"
866 c write (iout,*) "Calling sd_verlet_setup: 2"
867 c Restore the matrices of tinker integrator if the time step has been restored
871 pfric_mat(i,j)=pfric0_mat(i,j,0)
872 afric_mat(i,j)=afric0_mat(i,j,0)
873 vfric_mat(i,j)=vfric0_mat(i,j,0)
874 prand_mat(i,j)=prand0_mat(i,j,0)
875 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
876 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
885 c Calculate the kinetic and the total energy and the kinetic temperature
890 c write (iout,*) "step",itime," EK",EK," EK1",EK1
892 c Couple the system to Berendsen bath if needed
893 if (tbf .and. lang.eq.0) then
896 kinetic_T=2.0d0/(dimen*Rb)*EK
897 c Backup the coordinates, velocities, and accelerations
901 d_t_old(j,i)=d_t(j,i)
902 d_a_old(j,i)=d_a(j,i)
906 if (mod(itime,ntwe).eq.0 .and. large) then
907 write (iout,*) "Velocities, step 2"
909 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
910 & (d_t(j,i+nres),j=1,3)
916 c-------------------------------------------------------------------------------
917 subroutine RESPA_step(itime)
918 c-------------------------------------------------------------------------------
919 c Perform a single RESPA step.
920 c-------------------------------------------------------------------------------
921 implicit real*8 (a-h,o-z)
926 include 'COMMON.CONTROL'
930 include 'COMMON.LANGEVIN'
932 include 'COMMON.LANGEVIN.lang0'
934 include 'COMMON.CHAIN'
935 include 'COMMON.DERIV'
937 include 'COMMON.LOCAL'
938 include 'COMMON.INTERACT'
939 include 'COMMON.IOUNITS'
940 include 'COMMON.NAMES'
941 include 'COMMON.TIME1'
942 double precision energia_short(0:n_ene),
943 & energia_long(0:n_ene)
944 double precision cm(3),L(3),vcm(3),incr(3)
945 integer ilen,count,rstcount
948 integer maxcount_scale /10/
950 double precision stochforcvec(MAXRES6)
951 common /stochcalc/ stochforcvec
954 common /cipiszcze/ itt
956 large=(itime.gt.14600 .and. itime.lt.14700)
958 if (large.and. mod(itime,ntwe).eq.0) then
959 write (iout,*) "***************** RESPA itime",itime
960 write (iout,*) "Cartesian and internal coordinates: step 0"
962 call pdbout(0.0d0,"cipiszcze",iout)
964 write (iout,*) "Accelerations"
966 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
967 & (d_a(j,i+nres),j=1,3)
969 write (iout,*) "Velocities, step 0"
971 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
972 & (d_t(j,i+nres),j=1,3)
977 c Perform the initial RESPA step (increment velocities)
978 c write (iout,*) "*********************** RESPA ini"
981 if (mod(itime,ntwe).eq.0 .and. large) then
982 write (iout,*) "Velocities, end"
984 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
985 & (d_t(j,i+nres),j=1,3)
989 c Compute the short-range forces
996 call etotal_short(energia_short)
997 if (large) write (iout,*) "energia_short",energia_short(0)
1001 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1003 t_enegrad=t_enegrad+tcpu()-tt0
1008 d_t_old(j,i)=d_t(j,i)
1009 d_a_old(j,i)=d_a(j,i)
1013 c Split the time step
1014 d_time=d_time/ntime_split
1015 c Perform the short-range RESPSA steps (velocity Verlet increments of
1016 c positions and velocities using short-range forces)
1017 c write (iout,*) "*********************** RESPA split"
1018 do itsplit=1,ntime_split
1021 else if (lang.eq.2 .or. lang.eq.3) then
1022 call stochastic_force(stochforcvec)
1024 c First step of the velocity Verlet algorithm
1027 else if (lang.eq.3) then
1028 call sd_verlet1_ciccotti
1029 else if (lang.eq.1) then
1034 c Build the chain from the newly calculated coordinates
1035 call chainbuild_cart
1036 if (rattle) call rattle1
1038 if (large.and. mod(itime,ntwe).eq.0) then
1039 write (iout,*) "Cartesian and internal coordinates: step 1"
1040 call pdbout(0.0d0,"cipiszcze",iout)
1043 write (iout,*) "Accelerations"
1045 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1046 & (d_a(j,i+nres),j=1,3)
1048 write (iout,*) "Velocities, step 1"
1050 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1051 & (d_t(j,i+nres),j=1,3)
1060 c Calculate energy and forces
1062 call etotal_short(energia_short)
1063 if (large) write (iout,*) "energia_short",energia_short(0)
1065 c Get the new accelerations
1068 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1070 t_enegrad=t_enegrad+tcpu()-tt0
1072 c Second step of the velocity Verlet algorithm
1075 else if (lang.eq.3) then
1076 call sd_verlet2_ciccotti
1077 else if (lang.eq.1) then
1082 if (rattle) call rattle2
1083 c Backup the coordinates, velocities, and accelerations
1087 d_t_old(j,i)=d_t(j,i)
1088 d_a_old(j,i)=d_a(j,i)
1092 c Restore the time step
1094 c Compute long-range forces
1101 call etotal_long(energia_long)
1102 if (large) write (iout,*) "energia_long",energia_long(0)
1106 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1108 t_enegrad=t_enegrad+tcpu()-tt0
1110 c Compute accelerations from long-range forces
1112 if (large.and. mod(itime,ntwe).eq.0) then
1113 write (iout,*) "Cartesian and internal coordinates: step 2"
1115 call pdbout(0.0d0,"cipiszcze",iout)
1117 write (iout,*) "Accelerations"
1119 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1120 & (d_a(j,i+nres),j=1,3)
1122 write (iout,*) "Velocities, step 2"
1124 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1125 & (d_t(j,i+nres),j=1,3)
1129 c Compute the final RESPA step (increment velocities)
1130 c write (iout,*) "*********************** RESPA fin"
1132 c Compute the complete potential energy
1134 potEcomp(i)=energia_short(i)+energia_long(i)
1136 potE=potEcomp(0)-potEcomp(20)
1137 c potE=energia_short(0)+energia_long(0)
1139 c Calculate the kinetic and the total energy and the kinetic temperature
1142 c Couple the system to Berendsen bath if needed
1143 if (tbf .and. lang.eq.0) then
1146 kinetic_T=2.0d0/(dimen*Rb)*EK
1147 c Backup the coordinates, velocities, and accelerations
1149 if (mod(itime,ntwe).eq.0 .and. large) then
1150 write (iout,*) "Velocities, end"
1152 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1153 & (d_t(j,i+nres),j=1,3)
1159 c---------------------------------------------------------------------
1160 subroutine RESPA_vel
1161 c First and last RESPA step (incrementing velocities using long-range
1163 implicit real*8 (a-h,o-z)
1164 include 'DIMENSIONS'
1165 include 'COMMON.CONTROL'
1166 include 'COMMON.VAR'
1168 include 'COMMON.CHAIN'
1169 include 'COMMON.DERIV'
1170 include 'COMMON.GEO'
1171 include 'COMMON.LOCAL'
1172 include 'COMMON.INTERACT'
1173 include 'COMMON.IOUNITS'
1174 include 'COMMON.NAMES'
1176 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1180 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1184 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1187 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1193 c-----------------------------------------------------------------
1195 c Applying velocity Verlet algorithm - step 1 to coordinates
1196 implicit real*8 (a-h,o-z)
1197 include 'DIMENSIONS'
1198 include 'COMMON.CONTROL'
1199 include 'COMMON.VAR'
1201 include 'COMMON.CHAIN'
1202 include 'COMMON.DERIV'
1203 include 'COMMON.GEO'
1204 include 'COMMON.LOCAL'
1205 include 'COMMON.INTERACT'
1206 include 'COMMON.IOUNITS'
1207 include 'COMMON.NAMES'
1208 double precision adt,adt2
1211 adt=d_a_old(j,0)*d_time
1213 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1214 d_t_new(j,0)=d_t_old(j,0)+adt2
1215 d_t(j,0)=d_t_old(j,0)+adt
1219 adt=d_a_old(j,i)*d_time
1221 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1222 d_t_new(j,i)=d_t_old(j,i)+adt2
1223 d_t(j,i)=d_t_old(j,i)+adt
1227 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1230 adt=d_a_old(j,inres)*d_time
1232 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1233 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1234 d_t(j,inres)=d_t_old(j,inres)+adt
1240 c---------------------------------------------------------------------
1242 c Step 2 of the velocity Verlet algorithm: update velocities
1243 implicit real*8 (a-h,o-z)
1244 include 'DIMENSIONS'
1245 include 'COMMON.CONTROL'
1246 include 'COMMON.VAR'
1248 include 'COMMON.CHAIN'
1249 include 'COMMON.DERIV'
1250 include 'COMMON.GEO'
1251 include 'COMMON.LOCAL'
1252 include 'COMMON.INTERACT'
1253 include 'COMMON.IOUNITS'
1254 include 'COMMON.NAMES'
1256 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1260 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1264 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1267 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1273 c-----------------------------------------------------------------
1274 subroutine sddir_precalc
1275 c Applying velocity Verlet algorithm - step 1 to coordinates
1276 implicit real*8 (a-h,o-z)
1277 include 'DIMENSIONS'
1278 include 'COMMON.CONTROL'
1279 include 'COMMON.VAR'
1282 include 'COMMON.LANGEVIN'
1284 include 'COMMON.LANGEVIN.lang0'
1286 include 'COMMON.CHAIN'
1287 include 'COMMON.DERIV'
1288 include 'COMMON.GEO'
1289 include 'COMMON.LOCAL'
1290 include 'COMMON.INTERACT'
1291 include 'COMMON.IOUNITS'
1292 include 'COMMON.NAMES'
1293 double precision stochforcvec(MAXRES6)
1294 common /stochcalc/ stochforcvec
1296 c Compute friction and stochastic forces
1299 call stochastic_force(stochforcvec)
1301 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1302 c forces (d_as_work)
1309 d_af_work(i)=d_af_work(i)+Ginv(i,j)*fric_work(j)
1310 d_as_work(i)=d_as_work(i)+Ginv(i,j)*stochforcvec(j)
1314 call ginv_mult(fric_work, d_af_work)
1315 call ginv_mult(stochforcvec, d_as_work)
1319 c---------------------------------------------------------------------
1320 subroutine sddir_verlet1
1321 c Applying velocity Verlet algorithm - step 1 to velocities
1322 implicit real*8 (a-h,o-z)
1323 include 'DIMENSIONS'
1324 include 'COMMON.CONTROL'
1325 include 'COMMON.VAR'
1328 include 'COMMON.LANGEVIN'
1330 include 'COMMON.LANGEVIN.lang0'
1332 include 'COMMON.CHAIN'
1333 include 'COMMON.DERIV'
1334 include 'COMMON.GEO'
1335 include 'COMMON.LOCAL'
1336 include 'COMMON.INTERACT'
1337 include 'COMMON.IOUNITS'
1338 include 'COMMON.NAMES'
1339 c Revised 3/31/05 AL: correlation between random contributions to
1340 c position and velocity increments included.
1341 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1342 double precision adt,adt2
1344 c Add the contribution from BOTH friction and stochastic force to the
1345 c coordinates, but ONLY the contribution from the friction forces to velocities
1348 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1349 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1350 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1351 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1352 d_t(j,0)=d_t_old(j,0)+adt
1357 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1358 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1359 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1360 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1361 d_t(j,i)=d_t_old(j,i)+adt
1366 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1369 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1370 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1371 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1372 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1373 d_t(j,inres)=d_t_old(j,inres)+adt
1380 c---------------------------------------------------------------------
1381 subroutine sddir_verlet2
1382 c Calculating the adjusted velocities for accelerations
1383 implicit real*8 (a-h,o-z)
1384 include 'DIMENSIONS'
1385 include 'COMMON.CONTROL'
1386 include 'COMMON.VAR'
1389 include 'COMMON.LANGEVIN'
1391 include 'COMMON.LANGEVIN.lang0'
1393 include 'COMMON.CHAIN'
1394 include 'COMMON.DERIV'
1395 include 'COMMON.GEO'
1396 include 'COMMON.LOCAL'
1397 include 'COMMON.INTERACT'
1398 include 'COMMON.IOUNITS'
1399 include 'COMMON.NAMES'
1400 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1401 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1402 c Revised 3/31/05 AL: correlation between random contributions to
1403 c position and velocity increments included.
1404 c The correlation coefficients are calculated at low-friction limit.
1405 c Also, friction forces are now not calculated with new velocities.
1407 c call friction_force
1408 call stochastic_force(stochforcvec)
1410 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1411 c forces (d_as_work)
1415 c d_af_work(i)=0.0d0
1418 c d_af_work(i)=d_af_work(i)+Ginv(i,j)*fric_work(j)
1419 d_as_work1(i)=d_as_work1(i)+Ginv(i,j)*stochforcvec(j)
1423 call ginv_mult(stochforcvec, d_as_work1)
1430 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1431 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1436 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1437 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1442 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1445 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1446 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1447 & +cos60*d_as_work1(ind+j))*d_time
1454 c---------------------------------------------------------------------
1455 subroutine max_accel
1457 c Find the maximum difference in the accelerations of the the sites
1458 c at the beginning and the end of the time step.
1460 implicit real*8 (a-h,o-z)
1461 include 'DIMENSIONS'
1462 include 'COMMON.CONTROL'
1463 include 'COMMON.VAR'
1465 include 'COMMON.CHAIN'
1466 include 'COMMON.DERIV'
1467 include 'COMMON.GEO'
1468 include 'COMMON.LOCAL'
1469 include 'COMMON.INTERACT'
1470 include 'COMMON.IOUNITS'
1471 double precision aux(3),accel(3)
1473 aux(j)=d_a(j,0)-d_a_old(j,0)
1480 accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1481 if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1488 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1490 accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1494 if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1497 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
1507 implicit real*8 (a-h,o-z)
1508 include 'DIMENSIONS'
1509 include 'COMMON.CONTROL'
1510 include 'COMMON.VAR'
1512 include 'COMMON.CHAIN'
1513 include 'COMMON.DERIV'
1514 include 'COMMON.GEO'
1515 include 'COMMON.LOCAL'
1516 include 'COMMON.INTERACT'
1517 include 'COMMON.IOUNITS'
1518 include 'COMMON.MUCA'
1519 double precision epdrift,epdriftij
1520 c Drift of the potential energy
1526 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1527 if (lmuca) epdriftij=epdriftij*factor
1528 c write (iout,*) "back",i,j,epdriftij
1529 if (epdriftij.gt.epdrift) epdrift=epdriftij
1533 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1536 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1537 if (lmuca) epdriftij=epdriftij*factor
1538 c write (iout,*) "side",i,j,epdriftij
1539 if (epdriftij.gt.epdrift) epdrift=epdriftij
1543 epdrift=0.5d0*epdrift*d_time*d_time
1544 c write (iout,*) "epdrift",epdrift
1547 c-----------------------------------------------------------------------
1548 subroutine verlet_bath
1550 c Coupling to the thermostat by using the Berendsen algorithm
1552 implicit real*8 (a-h,o-z)
1553 include 'DIMENSIONS'
1554 include 'COMMON.CONTROL'
1555 include 'COMMON.VAR'
1557 include 'COMMON.CHAIN'
1558 include 'COMMON.DERIV'
1559 include 'COMMON.GEO'
1560 include 'COMMON.LOCAL'
1561 include 'COMMON.INTERACT'
1562 include 'COMMON.IOUNITS'
1563 include 'COMMON.NAMES'
1564 double precision T_half,fact
1566 T_half=2.0d0/(dimen*Rb)*EK
1567 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1568 c write(iout,*) "T_half", T_half
1569 c write(iout,*) "EK", EK
1570 c write(iout,*) "fact", fact
1572 d_t(j,0)=fact*d_t(j,0)
1576 d_t(j,i)=fact*d_t(j,i)
1580 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1583 d_t(j,inres)=fact*d_t(j,inres)
1589 c---------------------------------------------------------
1591 c Set up the initial conditions of a MD simulation
1592 implicit real*8 (a-h,o-z)
1593 include 'DIMENSIONS'
1596 include 'COMMON.SETUP'
1599 include 'COMMON.CONTROL'
1600 include 'COMMON.VAR'
1603 include 'COMMON.LANGEVIN'
1605 include 'COMMON.LANGEVIN.lang0'
1607 include 'COMMON.CHAIN'
1608 include 'COMMON.DERIV'
1609 include 'COMMON.GEO'
1610 include 'COMMON.LOCAL'
1611 include 'COMMON.INTERACT'
1612 include 'COMMON.IOUNITS'
1613 include 'COMMON.NAMES'
1614 include 'COMMON.REMD'
1615 real*8 energia_long(0:n_ene),
1616 & energia_short(0:n_ene),vcm(3),incr(3)
1617 double precision cm(3),L(3),xv,sigv,lowb,highb
1618 double precision varia(maxvar)
1626 c write(iout,*) "d_time", d_time
1627 c Compute the standard deviations of stochastic forces for Langevin dynamics
1628 c if the friction coefficients do not depend on surface area
1629 if (lang.gt.0 .and. .not.surfarea) then
1631 stdforcp(i)=stdfp*dsqrt(gamp)
1634 if (itype(i).ne.ntyp1)
1635 & stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(iabs(itype(i))))
1638 c Open the pdb file for snapshotshots
1641 if (ilen(tmpdir).gt.0)
1642 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1643 & liczba(:ilen(liczba))//".pdb")
1645 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1649 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1650 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1651 & liczba(:ilen(liczba))//".x")
1652 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1655 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1656 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1657 & liczba(:ilen(liczba))//".cx")
1658 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1664 if (ilen(tmpdir).gt.0)
1665 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1666 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1668 if (ilen(tmpdir).gt.0)
1669 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1670 cartname=prefix(:ilen(prefix))//"_MD.cx"
1674 write (qstr,'(256(1h ))')
1677 iq = qinfrag(i,iset)*10
1678 iw = wfrag(i,iset)/100
1680 if(me.eq.king.or..not.out1file)
1681 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1682 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1687 iq = qinpair(i,iset)*10
1688 iw = wpair(i,iset)/100
1690 if(me.eq.king.or..not.out1file)
1691 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1692 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1696 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1698 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1700 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1702 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1706 if (restart1file) then
1708 & inquire(file=mremd_rst_name,exist=file_exist)
1709 write (*,*) me," Before broadcast: file_exist",file_exist
1710 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1712 write (*,*) me," After broadcast: file_exist",file_exist
1713 c inquire(file=mremd_rst_name,exist=file_exist)
1714 if(me.eq.king.or..not.out1file)
1715 & write(iout,*) "Initial state read by master and distributed"
1717 if (ilen(tmpdir).gt.0)
1718 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1719 & //liczba(:ilen(liczba))//'.rst')
1720 inquire(file=rest2name,exist=file_exist)
1723 if(.not.restart1file) then
1724 if(me.eq.king.or..not.out1file)
1725 & write(iout,*) "Initial state will be read from file ",
1726 & rest2name(:ilen(rest2name))
1729 call rescale_weights(t_bath)
1731 if(me.eq.king.or..not.out1file)then
1732 if (restart1file) then
1733 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1736 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1739 write(iout,*) "Initial velocities randomly generated"
1745 c Generate initial velocities
1746 if(me.eq.king.or..not.out1file)
1747 & write(iout,*) "Initial velocities randomly generated"
1751 c rest2name = prefix(:ilen(prefix))//'.rst'
1752 if(me.eq.king.or..not.out1file)then
1753 write(iout,*) "Initial backbone velocities"
1755 write(iout,*) (d_t(j,i),j=1,3)
1757 write(iout,*) "Initial side-chain velocities"
1759 write(iout,*) (d_t(j,i+nres),j=1,3)
1761 c Zeroing the total angular momentum of the system
1762 write(iout,*) "Calling the zero-angular
1763 & momentum subroutine"
1766 c Getting the potential energy and forces and velocities and accelerations
1768 c write (iout,*) "velocity of the center of the mass:"
1769 c write (iout,*) (vcm(j),j=1,3)
1771 d_t(j,0)=d_t(j,0)-vcm(j)
1773 c Removing the velocity of the center of mass
1775 if(me.eq.king.or..not.out1file)then
1776 write (iout,*) "vcm right after adjustment:"
1777 write (iout,*) (vcm(j),j=1,3)
1781 if(iranconf.ne.0) then
1783 print *, 'Calling OVERLAP_SC'
1784 call overlap_sc(fail)
1788 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1789 print *,'SC_move',nft_sc,etot
1790 if(me.eq.king.or..not.out1file)
1791 & write(iout,*) 'SC_move',nft_sc,etot
1795 print *, 'Calling MINIM_DC'
1796 call minim_dc(etot,iretcode,nfun)
1798 call geom_to_var(nvar,varia)
1799 print *,'Calling MINIMIZE.'
1800 call minimize(etot,varia,iretcode,nfun)
1801 call var_to_geom(nvar,varia)
1803 if(me.eq.king.or..not.out1file)
1804 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1807 call chainbuild_cart
1810 call verlet_bath(EK)
1812 kinetic_T=2.0d0/(dimen*Rb)*EK
1813 if(me.eq.king.or..not.out1file)then
1818 call etotal(potEcomp)
1823 if (amax*d_time .gt. dvmax) d_time=d_time*dvmax/amax
1824 if(me.eq.king.or..not.out1file)then
1825 write(iout,*) "Potential energy and its components"
1826 write(iout,*) (potEcomp(i),i=0,n_ene)
1828 potE=potEcomp(0)-potEcomp(20)
1831 if (ntwe.ne.0) call statout(itime)
1832 if(me.eq.king.or..not.out1file)
1833 & write (iout,*) "Initial:",
1834 & " Kinetic energy",EK," potential energy",potE,
1835 & " total energy",totE," maximum acceleration ",
1838 write (iout,*) "Initial coordinates"
1840 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1841 & (c(j,i+nres),j=1,3)
1843 write (iout,*) "Initial dC"
1845 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1846 & (dc(j,i+nres),j=1,3)
1848 write (iout,*) "Initial velocities"
1850 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1851 & (d_t(j,i+nres),j=1,3)
1853 write (iout,*) "Initial accelerations"
1855 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1856 & (d_a(j,i+nres),j=1,3)
1862 d_t_old(j,i)=d_t(j,i)
1863 d_a_old(j,i)=d_a(j,i)
1865 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1874 call etotal_long(energia_long)
1875 if (large) write (iout,*) "energia_long",energia_long(0)
1879 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1881 t_enegrad=t_enegrad+tcpu()-tt0
1883 c call etotal_short(energia_short)
1884 c write (iout,*) "energia_long",energia_long(0),
1885 c & " energia_short",energia_short(0),
1886 c & " total",energia_long(0)+energia_short(0)
1890 c-----------------------------------------------------------
1891 subroutine random_vel
1892 implicit real*8 (a-h,o-z)
1893 include 'DIMENSIONS'
1894 include 'COMMON.CONTROL'
1895 include 'COMMON.VAR'
1898 include 'COMMON.LANGEVIN'
1900 include 'COMMON.LANGEVIN.lang0'
1902 include 'COMMON.CHAIN'
1903 include 'COMMON.DERIV'
1904 include 'COMMON.GEO'
1905 include 'COMMON.LOCAL'
1906 include 'COMMON.INTERACT'
1907 include 'COMMON.IOUNITS'
1908 include 'COMMON.NAMES'
1909 include 'COMMON.TIME1'
1910 double precision xv,sigv,lowb,highb
1911 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1912 c First generate velocities in the eigenspace of the G matrix
1913 c write (iout,*) "Calling random_vel"
1916 sigv=dsqrt((Rb*t_bath)/geigen(i))
1919 d_t_work_new(i)=anorm_distr(xv,sigv,lowb,highb)
1923 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(i)**2
1925 c Transform velocities to UNRES coordinate space
1929 d_t_work(i)=d_t_work(i)+Gvec(i,j)*d_t_work_new(j)
1932 c Transfer to the d_t vector
1934 d_t(j,0)=d_t_work(j)
1940 if (itype(i).ne.ntyp1 .and. itype(i+1).ne.ntyp1) then
1941 d_t(j,i)=d_t_work(ind)
1948 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1951 d_t(j,i+nres)=d_t_work(ind)
1956 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1957 c & 2.0d0/(dimen*Rb)*EK,2.0d0/(dimen*Rb)*EK1
1960 c-----------------------------------------------------------
1961 subroutine sd_verlet_p_setup
1962 c Sets up the parameters of stochastic Verlet algorithm
1963 implicit real*8 (a-h,o-z)
1964 include 'DIMENSIONS'
1968 include 'COMMON.CONTROL'
1969 include 'COMMON.VAR'
1972 include 'COMMON.LANGEVIN'
1974 include 'COMMON.LANGEVIN.lang0'
1976 include 'COMMON.CHAIN'
1977 include 'COMMON.DERIV'
1978 include 'COMMON.GEO'
1979 include 'COMMON.LOCAL'
1980 include 'COMMON.INTERACT'
1981 include 'COMMON.IOUNITS'
1982 include 'COMMON.NAMES'
1983 include 'COMMON.TIME1'
1984 double precision emgdt(MAXRES6),
1985 & pterm,vterm,rho,rhoc,vsig,
1986 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1987 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1988 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1989 logical lprn /.false./
1990 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1991 double precision ktm
1998 c AL 8/17/04 Code adapted from tinker
2000 c Get the frictional and random terms for stochastic dynamics in the
2001 c eigenspace of mass-scaled UNRES friction matrix
2004 gdt = fricgam(i) * d_time
2006 c Stochastic dynamics reduces to simple MD for zero friction
2008 if (gdt .le. zero) then
2009 pfric_vec(i) = 1.0d0
2010 vfric_vec(i) = d_time
2011 afric_vec(i) = 0.5d0 * d_time * d_time
2012 prand_vec(i) = 0.0d0
2013 vrand_vec1(i) = 0.0d0
2014 vrand_vec2(i) = 0.0d0
2016 c Analytical expressions when friction coefficient is large
2019 if (gdt .ge. gdt_radius) then
2022 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2023 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2024 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2025 vterm = 1.0d0 - egdt**2
2026 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2028 c Use series expansions when friction coefficient is small
2039 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2040 & - gdt5/120.0d0 + gdt6/720.0d0
2041 & - gdt7/5040.0d0 + gdt8/40320.0d0
2042 & - gdt9/362880.0d0) / fricgam(i)**2
2043 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2044 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2045 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2046 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2047 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2048 & + 127.0d0*gdt9/90720.0d0
2049 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2050 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2051 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2052 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2053 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2054 & - 17.0d0*gdt2/1280.0d0
2055 & + 17.0d0*gdt3/6144.0d0
2056 & + 40967.0d0*gdt4/34406400.0d0
2057 & - 57203.0d0*gdt5/275251200.0d0
2058 & - 1429487.0d0*gdt6/13212057600.0d0)
2061 c Compute the scaling factors of random terms for the nonzero friction case
2063 ktm = 0.5d0*d_time/fricgam(i)
2064 psig = dsqrt(ktm*pterm) / fricgam(i)
2065 vsig = dsqrt(ktm*vterm)
2066 rhoc = dsqrt(1.0d0 - rho*rho)
2068 vrand_vec1(i) = vsig * rho
2069 vrand_vec2(i) = vsig * rhoc
2074 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2077 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2078 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2082 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2085 call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
2086 call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
2087 call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
2088 call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
2089 call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec1,vrand_mat1)
2090 call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
2092 c call eigtransf1(dimen,maxres6,mt3mt2,pfric_vec,pfric_mat)
2093 c call eigtransf1(dimen,maxres6,mt3mt2,vfric_vec,vfric_mat)
2094 c call eigtransf1(dimen,maxres6,mt3mt2,afric_vec,afric_mat)
2095 c call eigtransf1(dimen,maxres6,mt3mt1,prand_vec,prand_mat)
2096 c call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec1,vrand_mat1)
2097 c call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec2,vrand_mat2)
2099 t_sdsetup=t_sdsetup+MPI_Wtime()
2101 t_sdsetup=t_sdsetup+tcpu()-tt0
2105 c-------------------------------------------------------------
2106 subroutine eigtransf1(n,ndim,ab,d,c)
2109 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2115 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2121 c-------------------------------------------------------------
2122 subroutine eigtransf(n,ndim,a,b,d,c)
2125 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2131 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2137 c-------------------------------------------------------------
2138 subroutine sd_verlet1
2139 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2140 implicit real*8 (a-h,o-z)
2141 include 'DIMENSIONS'
2142 include 'COMMON.CONTROL'
2143 include 'COMMON.VAR'
2146 include 'COMMON.LANGEVIN'
2148 include 'COMMON.LANGEVIN.lang0'
2150 include 'COMMON.CHAIN'
2151 include 'COMMON.DERIV'
2152 include 'COMMON.GEO'
2153 include 'COMMON.LOCAL'
2154 include 'COMMON.INTERACT'
2155 include 'COMMON.IOUNITS'
2156 include 'COMMON.NAMES'
2157 double precision stochforcvec(MAXRES6)
2158 common /stochcalc/ stochforcvec
2159 logical lprn /.false./
2161 c write (iout,*) "dc_old"
2163 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2164 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2167 dc_work(j)=dc_old(j,0)
2168 d_t_work(j)=d_t_old(j,0)
2169 d_a_work(j)=d_a_old(j,0)
2174 dc_work(ind+j)=dc_old(j,i)
2175 d_t_work(ind+j)=d_t_old(j,i)
2176 d_a_work(ind+j)=d_a_old(j,i)
2181 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2183 dc_work(ind+j)=dc_old(j,i+nres)
2184 d_t_work(ind+j)=d_t_old(j,i+nres)
2185 d_a_work(ind+j)=d_a_old(j,i+nres)
2193 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2197 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2198 & vfric_mat(i,j),afric_mat(i,j),
2199 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2207 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2208 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2209 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2210 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2212 d_t_work_new(i)=ddt1+0.5d0*ddt2
2213 d_t_work(i)=ddt1+ddt2
2218 d_t(j,0)=d_t_work(j)
2223 dc(j,i)=dc_work(ind+j)
2224 d_t(j,i)=d_t_work(ind+j)
2229 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2232 dc(j,inres)=dc_work(ind+j)
2233 d_t(j,inres)=d_t_work(ind+j)
2240 c--------------------------------------------------------------------------
2241 subroutine sd_verlet2
2242 c Calculating the adjusted velocities for accelerations
2243 implicit real*8 (a-h,o-z)
2244 include 'DIMENSIONS'
2245 include 'COMMON.CONTROL'
2246 include 'COMMON.VAR'
2249 include 'COMMON.LANGEVIN'
2251 include 'COMMON.LANGEVIN.lang0'
2253 include 'COMMON.CHAIN'
2254 include 'COMMON.DERIV'
2255 include 'COMMON.GEO'
2256 include 'COMMON.LOCAL'
2257 include 'COMMON.INTERACT'
2258 include 'COMMON.IOUNITS'
2259 include 'COMMON.NAMES'
2260 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2261 common /stochcalc/ stochforcvec
2263 c Compute the stochastic forces which contribute to velocity change
2265 call stochastic_force(stochforcvecV)
2272 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2273 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2274 & vrand_mat2(i,j)*stochforcvecV(j)
2276 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2280 d_t(j,0)=d_t_work(j)
2285 d_t(j,i)=d_t_work(ind+j)
2290 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2293 d_t(j,inres)=d_t_work(ind+j)
2300 c-----------------------------------------------------------
2301 subroutine sd_verlet_ciccotti_setup
2302 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2304 implicit real*8 (a-h,o-z)
2305 include 'DIMENSIONS'
2309 include 'COMMON.CONTROL'
2310 include 'COMMON.VAR'
2313 include 'COMMON.LANGEVIN'
2315 include 'COMMON.LANGEVIN.lang0'
2317 include 'COMMON.CHAIN'
2318 include 'COMMON.DERIV'
2319 include 'COMMON.GEO'
2320 include 'COMMON.LOCAL'
2321 include 'COMMON.INTERACT'
2322 include 'COMMON.IOUNITS'
2323 include 'COMMON.NAMES'
2324 include 'COMMON.TIME1'
2325 double precision emgdt(MAXRES6),
2326 & pterm,vterm,rho,rhoc,vsig,
2327 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2328 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2329 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2330 logical lprn /.false./
2331 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2332 double precision ktm
2339 c AL 8/17/04 Code adapted from tinker
2341 c Get the frictional and random terms for stochastic dynamics in the
2342 c eigenspace of mass-scaled UNRES friction matrix
2345 write (iout,*) "i",i," fricgam",fricgam(i)
2346 gdt = fricgam(i) * d_time
2348 c Stochastic dynamics reduces to simple MD for zero friction
2350 if (gdt .le. zero) then
2351 pfric_vec(i) = 1.0d0
2352 vfric_vec(i) = d_time
2353 afric_vec(i) = 0.5d0*d_time*d_time
2354 prand_vec(i) = afric_vec(i)
2355 vrand_vec2(i) = vfric_vec(i)
2357 c Analytical expressions when friction coefficient is large
2362 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2363 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2364 prand_vec(i) = afric_vec(i)
2365 vrand_vec2(i) = vfric_vec(i)
2367 c Compute the scaling factors of random terms for the nonzero friction case
2369 c ktm = 0.5d0*d_time/fricgam(i)
2370 c psig = dsqrt(ktm*pterm) / fricgam(i)
2371 c vsig = dsqrt(ktm*vterm)
2372 c prand_vec(i) = psig*afric_vec(i)
2373 c vrand_vec2(i) = vsig*vfric_vec(i)
2378 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2381 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2382 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2386 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2388 call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
2389 call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
2390 call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
2391 call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
2392 call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
2394 t_sdsetup=t_sdsetup+MPI_Wtime()
2396 t_sdsetup=t_sdsetup+tcpu()-tt0
2400 c-------------------------------------------------------------
2401 subroutine sd_verlet1_ciccotti
2402 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2403 implicit real*8 (a-h,o-z)
2404 include 'DIMENSIONS'
2408 include 'COMMON.CONTROL'
2409 include 'COMMON.VAR'
2412 include 'COMMON.LANGEVIN'
2414 include 'COMMON.LANGEVIN.lang0'
2416 include 'COMMON.CHAIN'
2417 include 'COMMON.DERIV'
2418 include 'COMMON.GEO'
2419 include 'COMMON.LOCAL'
2420 include 'COMMON.INTERACT'
2421 include 'COMMON.IOUNITS'
2422 include 'COMMON.NAMES'
2423 double precision stochforcvec(MAXRES6)
2424 common /stochcalc/ stochforcvec
2425 logical lprn /.false./
2427 c write (iout,*) "dc_old"
2429 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2430 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2433 dc_work(j)=dc_old(j,0)
2434 d_t_work(j)=d_t_old(j,0)
2435 d_a_work(j)=d_a_old(j,0)
2440 dc_work(ind+j)=dc_old(j,i)
2441 d_t_work(ind+j)=d_t_old(j,i)
2442 d_a_work(ind+j)=d_a_old(j,i)
2447 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2449 dc_work(ind+j)=dc_old(j,i+nres)
2450 d_t_work(ind+j)=d_t_old(j,i+nres)
2451 d_a_work(ind+j)=d_a_old(j,i+nres)
2460 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2464 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2465 & vfric_mat(i,j),afric_mat(i,j),
2466 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2474 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2475 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2476 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2477 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2479 d_t_work_new(i)=ddt1+0.5d0*ddt2
2480 d_t_work(i)=ddt1+ddt2
2485 d_t(j,0)=d_t_work(j)
2490 dc(j,i)=dc_work(ind+j)
2491 d_t(j,i)=d_t_work(ind+j)
2496 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2499 dc(j,inres)=dc_work(ind+j)
2500 d_t(j,inres)=d_t_work(ind+j)
2507 c--------------------------------------------------------------------------
2508 subroutine sd_verlet2_ciccotti
2509 c Calculating the adjusted velocities for accelerations
2510 implicit real*8 (a-h,o-z)
2511 include 'DIMENSIONS'
2512 include 'COMMON.CONTROL'
2513 include 'COMMON.VAR'
2516 include 'COMMON.LANGEVIN'
2518 include 'COMMON.LANGEVIN.lang0'
2520 include 'COMMON.CHAIN'
2521 include 'COMMON.DERIV'
2522 include 'COMMON.GEO'
2523 include 'COMMON.LOCAL'
2524 include 'COMMON.INTERACT'
2525 include 'COMMON.IOUNITS'
2526 include 'COMMON.NAMES'
2527 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2528 common /stochcalc/ stochforcvec
2530 c Compute the stochastic forces which contribute to velocity change
2532 call stochastic_force(stochforcvecV)
2539 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2540 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2541 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2543 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2547 d_t(j,0)=d_t_work(j)
2552 d_t(j,i)=d_t_work(ind+j)
2557 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2560 d_t(j,inres)=d_t_work(ind+j)