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 double precision cm(3),L(3),vcm(3)
28 double precision v_work(maxres6),v_transf(maxres6)
37 if (ilen(tmpdir).gt.0)
38 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
39 & //liczba(:ilen(liczba))//'.rst')
41 if (ilen(tmpdir).gt.0)
42 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
49 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
55 c Determine the inverse of the inertia matrix.
56 call setup_MD_matrices
60 t_MDsetup = MPI_Wtime()-tt0
62 t_MDsetup = tcpu()-tt0
65 c Entering the MD loop
71 if (lang.eq.2 .or. lang.eq.3) then
74 call sd_verlet_p_setup
76 call sd_verlet_ciccotti_setup
81 pfric0_mat(i,j,0)=pfric_mat(i,j)
82 afric0_mat(i,j,0)=afric_mat(i,j)
83 vfric0_mat(i,j,0)=vfric_mat(i,j)
84 prand0_mat(i,j,0)=prand_mat(i,j)
85 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
86 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
94 else if (lang.eq.1 .or. lang.eq.4) then
98 t_langsetup=MPI_Wtime()-tt0
101 t_langsetup=tcpu()-tt0
104 do itime=1,n_timestep
105 if (ovrtim()) goto 100
107 if (lang.gt.0 .and. surfarea .and.
108 & mod(itime,reset_fricmat).eq.0) then
109 if (lang.eq.2 .or. lang.eq.3) then
112 call sd_verlet_p_setup
114 call sd_verlet_ciccotti_setup
119 pfric0_mat(i,j,0)=pfric_mat(i,j)
120 afric0_mat(i,j,0)=afric_mat(i,j)
121 vfric0_mat(i,j,0)=vfric_mat(i,j)
122 prand0_mat(i,j,0)=prand_mat(i,j)
123 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
124 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
130 flag_stoch(i)=.false.
132 else if (lang.eq.1 .or. lang.eq.4) then
135 write (iout,'(a,i10)')
136 & "Friction matrix reset based on surface area, itime",itime
138 if (reset_vel .and. tbf .and. lang.eq.0
139 & .and. mod(itime,count_reset_vel).eq.0) then
141 write(iout,'(a,f20.2)')
142 & "Velocities reset to random values, time",totT
145 d_t_old(j,i)=d_t(j,i)
149 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
153 d_t(j,0)=d_t(j,0)-vcm(j)
156 kinetic_T=2.0d0/(dimen*Rb)*EK
157 scalfac=dsqrt(T_bath/kinetic_T)
158 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
161 d_t_old(j,i)=scalfac*d_t(j,i)
167 c Time-reversible RESPA algorithm
168 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
169 call RESPA_step(itime)
171 c Variable time step algorithm.
172 call velverlet_step(itime)
175 call brown_step(itime)
178 if (mod(itime,ntwe).eq.0) call statout(itime)
191 if (itype(i).ne.10) then
194 v_work(ind)=d_t(j,i+nres)
199 write (66,'(80f10.5)')
200 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
204 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
206 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
208 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
211 if (mod(itime,ntwx).eq.0) then
212 write (tytul,'("time",f8.2)') totT
214 call pdbout(potE,tytul,ipdb)
219 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
220 open(irest2,file=rest2name,status='unknown')
221 write(irest2,*) totT,EK,potE,totE,t_bath
223 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
226 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
238 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
240 & 'MD calculations setup:',t_MDsetup,
241 & 'Energy & gradient evaluation:',t_enegrad,
242 & 'Stochastic MD setup:',t_langsetup,
243 & 'Stochastic MD step setup:',t_sdsetup,
245 write (iout,'(/28(1h=),a25,27(1h=))')
246 & ' End of MD calculation '
249 c-------------------------------------------------------------------------------
250 subroutine brown_step(itime)
251 c------------------------------------------------
252 c Perform a single Euler integration step of Brownian dynamics
253 c------------------------------------------------
254 implicit real*8 (a-h,o-z)
259 include 'COMMON.CONTROL'
263 include 'COMMON.LANGEVIN'
265 include 'COMMON.LANGEVIN.lang0'
267 include 'COMMON.CHAIN'
268 include 'COMMON.DERIV'
270 include 'COMMON.LOCAL'
271 include 'COMMON.INTERACT'
272 include 'COMMON.IOUNITS'
273 include 'COMMON.NAMES'
274 include 'COMMON.TIME1'
275 double precision zapas(MAXRES6)
276 integer ilen,rstcount
278 double precision stochforcvec(MAXRES6)
279 double precision Bmat(MAXRES6,MAXRES2),Cmat(maxres2,maxres2),
280 & Cinv(maxres2,maxres2),GBmat(MAXRES6,MAXRES2),
281 & Tmat(MAXRES6,MAXRES2),Pmat(maxres6,maxres6),Td(maxres6),
283 common /stochcalc/ stochforcvec
286 logical lprn /.false./,lprn1 /.false./
288 double precision difftol /1.0d-5/
291 if (itype(i).ne.10) nbond=nbond+1
295 write (iout,*) "Generalized inverse of fricmat"
296 call matout(dimen,dimen,MAXRES6,MAXRES6,fricmat)
308 Bmat(ind+j,ind1)=dC_norm(j,i)
313 if (itype(i).ne.10) then
316 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
322 write (iout,*) "Matrix Bmat"
323 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Bmat)
329 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
334 write (iout,*) "Matrix GBmat"
335 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Gbmat)
341 Cmat(i,j)=Cmat(i,j)+Bmat(k,i)*GBmat(k,j)
346 write (iout,*) "Matrix Cmat"
347 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
349 call matinvert(nbond,MAXRES2,Cmat,Cinv)
351 write (iout,*) "Matrix Cinv"
352 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cinv)
358 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
363 write (iout,*) "Matrix Tmat"
364 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Tmat)
374 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
379 write (iout,*) "Matrix Pmat"
380 call MATOUT(dimen,dimen,MAXRES6,MAXRES6,Pmat)
387 Td(i)=Td(i)+vbl*Tmat(i,ind)
390 if (itype(k).ne.10) then
392 Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
397 write (iout,*) "Vector Td"
399 write (iout,'(i5,f10.5)') i,Td(i)
402 call stochastic_force(stochforcvec)
404 write (iout,*) "stochforcvec"
406 write (iout,*) i,stochforcvec(i)
410 zapas(j)=-gcart(j,0)+stochforcvec(j)
412 dC_work(j)=dC_old(j,0)
418 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
419 dC_work(ind)=dC_old(j,i)
423 if (itype(i).ne.10) then
426 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
427 dC_work(ind)=dC_old(j,i+nres)
433 write (iout,*) "Initial d_t_work"
435 write (iout,*) i,d_t_work(i)
442 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
449 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
453 write (iout,*) "Final d_t_work and zapas"
455 write (iout,*) i,d_t_work(i),zapas(i)
469 dc_work(ind+j)=dc(j,i)
475 d_t(j,i+nres)=d_t_work(ind+j)
476 dc(j,i+nres)=zapas(ind+j)
477 dc_work(ind+j)=dc(j,i+nres)
483 write (iout,*) "Before correction for rotational lengthening"
484 write (iout,*) "New coordinates",
485 & " and differences between actual and standard bond lengths"
490 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
491 & i,(dC(j,i),j=1,3),xx
494 if (itype(i).ne.10) then
496 xx=vbld(i+nres)-vbldsc0(1,itype(i))
497 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
498 & i,(dC(j,i+nres),j=1,3),xx
502 c Second correction (rotational lengthening)
508 blen2 = scalar(dc(1,i),dc(1,i))
509 ppvec(ind)=2*vbl**2-blen2
510 diffbond=dabs(vbl-dsqrt(blen2))
511 if (diffbond.gt.diffmax) diffmax=diffbond
512 if (ppvec(ind).gt.0.0d0) then
513 ppvec(ind)=dsqrt(ppvec(ind))
518 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
522 if (itype(i).ne.10) then
524 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
525 ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
526 diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2))
527 if (diffbond.gt.diffmax) diffmax=diffbond
528 if (ppvec(ind).gt.0.0d0) then
529 ppvec(ind)=dsqrt(ppvec(ind))
534 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
538 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
539 if (diffmax.lt.difftol) goto 10
543 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
549 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
560 dc_work(ind+j)=zapas(ind+j)
565 if (itype(i).ne.10) then
567 dc(j,i+nres)=zapas(ind+j)
568 dc_work(ind+j)=zapas(ind+j)
573 c Building the chain from the newly calculated coordinates
576 if (large.and. mod(itime,ntwe).eq.0) then
577 write (iout,*) "Cartesian and internal coordinates: step 1"
580 write (iout,'(a)') "Potential forces"
582 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),
583 & (-gxcart(j,i),j=1,3)
585 write (iout,'(a)') "Stochastic forces"
587 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),
588 & (stochforc(j,i+nres),j=1,3)
590 write (iout,'(a)') "Velocities"
592 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
593 & (d_t(j,i+nres),j=1,3)
598 write (iout,*) "After correction for rotational lengthening"
599 write (iout,*) "New coordinates",
600 & " and differences between actual and standard bond lengths"
605 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
606 & i,(dC(j,i),j=1,3),xx
609 if (itype(i).ne.10) then
611 xx=vbld(i+nres)-vbldsc0(1,itype(i))
612 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
613 & i,(dC(j,i+nres),j=1,3),xx
618 c write (iout,*) "Too many attempts at correcting the bonds"
626 c Calculate energy and forces
628 call etotal(potEcomp)
629 potE=potEcomp(0)-potEcomp(20)
632 c Calculate the kinetic and total energy and the kinetic temperature
635 t_enegrad=t_enegrad+MPI_Wtime()-tt0
637 t_enegrad=t_enegrad+tcpu()-tt0
640 kinetic_T=2.0d0/(dimen*Rb)*EK
643 c-------------------------------------------------------------------------------
644 subroutine velverlet_step(itime)
645 c-------------------------------------------------------------------------------
646 c Perform a single velocity Verlet step; the time step can be rescaled if
647 c increments in accelerations exceed the threshold
648 c-------------------------------------------------------------------------------
649 implicit real*8 (a-h,o-z)
653 integer ierror,ierrcode
655 include 'COMMON.CONTROL'
659 include 'COMMON.LANGEVIN'
661 include 'COMMON.LANGEVIN.lang0'
663 include 'COMMON.CHAIN'
664 include 'COMMON.DERIV'
666 include 'COMMON.LOCAL'
667 include 'COMMON.INTERACT'
668 include 'COMMON.IOUNITS'
669 include 'COMMON.NAMES'
670 include 'COMMON.TIME1'
671 include 'COMMON.MUCA'
672 double precision vcm(3),incr(3)
673 double precision cm(3),L(3)
674 integer ilen,count,rstcount
677 integer maxcount_scale /20/
679 double precision stochforcvec(MAXRES6)
680 common /stochcalc/ stochforcvec
688 else if (lang.eq.2 .or. lang.eq.3) then
689 call stochastic_force(stochforcvec)
693 icount_scale=icount_scale+1
694 if (icount_scale.gt.maxcount_scale) then
696 & "ERROR: too many attempts at scaling down the time step. ",
697 & "amax=",amax,"epdrift=",epdrift,
698 & "damax=",damax,"edriftmax=",edriftmax,
702 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
706 c First step of the velocity Verlet algorithm
709 else if (lang.eq.3) then
710 call sd_verlet1_ciccotti
711 else if (lang.eq.1) then
716 c Build the chain from the newly calculated coordinates
718 if (rattle) call rattle1
720 if (large.and. mod(itime,ntwe).eq.0) then
721 write (iout,*) "Cartesian and internal coordinates: step 1"
726 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
727 & (dc(j,i+nres),j=1,3)
729 write (iout,*) "Accelerations"
731 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
732 & (d_a(j,i+nres),j=1,3)
734 write (iout,*) "Velocities, step 1"
736 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
737 & (d_t(j,i+nres),j=1,3)
746 c Calculate energy and forces
748 call etotal(potEcomp)
749 potE=potEcomp(0)-potEcomp(20)
751 c Get the new accelerations
754 t_enegrad=t_enegrad+MPI_Wtime()-tt0
756 t_enegrad=t_enegrad+tcpu()-tt0
758 c Determine maximum acceleration and scale down the timestep if needed
760 call predict_edrift(epdrift)
761 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
762 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
764 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
766 itime_scal=itime_scal+ifac_time
767 c fac_time=dmin1(damax/amax,0.5d0)
768 fac_time=0.5d0**ifac_time
769 d_time=d_time*fac_time
770 if (lang.eq.2 .or. lang.eq.3) then
771 c write (iout,*) "Calling sd_verlet_setup: 1"
772 c Rescale the stochastic forces and recalculate or restore
773 c the matrices of tinker integrator
774 if (itime_scal.gt.maxflag_stoch) then
775 if (large) write (iout,'(a,i5,a)')
776 & "Calculate matrices for stochastic step;",
777 & " itime_scal ",itime_scal
779 call sd_verlet_p_setup
781 call sd_verlet_ciccotti_setup
783 write (iout,'(2a,i3,a,i3,1h.)')
784 & "Warning: cannot store matrices for stochastic",
785 & " integration because the index",itime_scal,
786 & " is greater than",maxflag_stoch
787 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
788 & " integration Langevin algorithm for better efficiency."
789 else if (flag_stoch(itime_scal)) then
790 if (large) write (iout,'(a,i5,a,l1)')
791 & "Restore matrices for stochastic step; itime_scal ",
792 & itime_scal," flag ",flag_stoch(itime_scal)
796 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
797 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
798 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
799 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
800 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
801 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
806 if (large) write (iout,'(2a,i5,a,l1)')
807 & "Calculate & store matrices for stochastic step;",
808 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
810 call sd_verlet_p_setup
812 call sd_verlet_ciccotti_setup
814 flag_stoch(ifac_time)=.true.
818 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
819 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
820 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
821 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
822 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
823 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
828 fac_time=1.0d0/dsqrt(fac_time)
830 stochforcvec(i)=fac_time*stochforcvec(i)
832 else if (lang.eq.1) then
833 c Rescale the accelerations due to stochastic forces
834 fac_time=1.0d0/dsqrt(fac_time)
836 d_as_work(i)=d_as_work(i)*fac_time
839 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
840 & "itime",itime," Timestep scaled down to ",
841 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
843 c Second step of the velocity Verlet algorithm
846 else if (lang.eq.3) then
847 call sd_verlet2_ciccotti
848 else if (lang.eq.1) then
853 if (rattle) call rattle2
855 if (d_time.ne.d_time0) then
857 if (lang.eq.2 .or. lang.eq.3) then
858 if (large) write (iout,'(a)')
859 & "Restore original matrices for stochastic step"
860 c write (iout,*) "Calling sd_verlet_setup: 2"
861 c Restore the matrices of tinker integrator if the time step has been restored
865 pfric_mat(i,j)=pfric0_mat(i,j,0)
866 afric_mat(i,j)=afric0_mat(i,j,0)
867 vfric_mat(i,j)=vfric0_mat(i,j,0)
868 prand_mat(i,j)=prand0_mat(i,j,0)
869 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
870 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
879 c Calculate the kinetic and the total energy and the kinetic temperature
884 c write (iout,*) "step",itime," EK",EK," EK1",EK1
886 c Couple the system to Berendsen bath if needed
887 if (tbf .and. lang.eq.0) then
890 kinetic_T=2.0d0/(dimen*Rb)*EK
891 c Backup the coordinates, velocities, and accelerations
895 d_t_old(j,i)=d_t(j,i)
896 d_a_old(j,i)=d_a(j,i)
900 if (mod(itime,ntwe).eq.0 .and. large) then
901 write (iout,*) "Velocities, step 2"
903 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
904 & (d_t(j,i+nres),j=1,3)
910 c-------------------------------------------------------------------------------
911 subroutine RESPA_step(itime)
912 c-------------------------------------------------------------------------------
913 c Perform a single RESPA step.
914 c-------------------------------------------------------------------------------
915 implicit real*8 (a-h,o-z)
920 include 'COMMON.CONTROL'
924 include 'COMMON.LANGEVIN'
926 include 'COMMON.LANGEVIN.lang0'
928 include 'COMMON.CHAIN'
929 include 'COMMON.DERIV'
931 include 'COMMON.LOCAL'
932 include 'COMMON.INTERACT'
933 include 'COMMON.IOUNITS'
934 include 'COMMON.NAMES'
935 include 'COMMON.TIME1'
936 double precision energia_short(0:n_ene),
937 & energia_long(0:n_ene)
938 double precision cm(3),L(3),vcm(3),incr(3)
939 integer ilen,count,rstcount
942 integer maxcount_scale /10/
944 double precision stochforcvec(MAXRES6)
945 common /stochcalc/ stochforcvec
948 common /cipiszcze/ itt
950 large=(itime.gt.14600 .and. itime.lt.14700)
952 if (large.and. mod(itime,ntwe).eq.0) then
953 write (iout,*) "***************** RESPA itime",itime
954 write (iout,*) "Cartesian and internal coordinates: step 0"
956 call pdbout(0.0d0,"cipiszcze",iout)
958 write (iout,*) "Accelerations"
960 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
961 & (d_a(j,i+nres),j=1,3)
963 write (iout,*) "Velocities, step 0"
965 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
966 & (d_t(j,i+nres),j=1,3)
971 c Perform the initial RESPA step (increment velocities)
972 c write (iout,*) "*********************** RESPA ini"
975 if (mod(itime,ntwe).eq.0 .and. large) then
976 write (iout,*) "Velocities, end"
978 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
979 & (d_t(j,i+nres),j=1,3)
983 c Compute the short-range forces
990 call etotal_short(energia_short)
991 if (large) write (iout,*) "energia_short",energia_short(0)
995 t_enegrad=t_enegrad+MPI_Wtime()-tt0
997 t_enegrad=t_enegrad+tcpu()-tt0
1002 d_t_old(j,i)=d_t(j,i)
1003 d_a_old(j,i)=d_a(j,i)
1007 c Split the time step
1008 d_time=d_time/ntime_split
1009 c Perform the short-range RESPSA steps (velocity Verlet increments of
1010 c positions and velocities using short-range forces)
1011 c write (iout,*) "*********************** RESPA split"
1012 do itsplit=1,ntime_split
1015 else if (lang.eq.2 .or. lang.eq.3) then
1016 call stochastic_force(stochforcvec)
1018 c First step of the velocity Verlet algorithm
1021 else if (lang.eq.3) then
1022 call sd_verlet1_ciccotti
1023 else if (lang.eq.1) then
1028 c Build the chain from the newly calculated coordinates
1029 call chainbuild_cart
1030 if (rattle) call rattle1
1032 if (large.and. mod(itime,ntwe).eq.0) then
1033 write (iout,*) "Cartesian and internal coordinates: step 1"
1034 call pdbout(0.0d0,"cipiszcze",iout)
1037 write (iout,*) "Accelerations"
1039 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1040 & (d_a(j,i+nres),j=1,3)
1042 write (iout,*) "Velocities, step 1"
1044 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1045 & (d_t(j,i+nres),j=1,3)
1054 c Calculate energy and forces
1056 call etotal_short(energia_short)
1057 if (large) write (iout,*) "energia_short",energia_short(0)
1059 c Get the new accelerations
1062 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1064 t_enegrad=t_enegrad+tcpu()-tt0
1066 c Second step of the velocity Verlet algorithm
1069 else if (lang.eq.3) then
1070 call sd_verlet2_ciccotti
1071 else if (lang.eq.1) then
1076 if (rattle) call rattle2
1077 c Backup the coordinates, velocities, and accelerations
1081 d_t_old(j,i)=d_t(j,i)
1082 d_a_old(j,i)=d_a(j,i)
1086 c Restore the time step
1088 c Compute long-range forces
1095 call etotal_long(energia_long)
1096 if (large) write (iout,*) "energia_long",energia_long(0)
1100 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1102 t_enegrad=t_enegrad+tcpu()-tt0
1104 c Compute accelerations from long-range forces
1106 if (large.and. mod(itime,ntwe).eq.0) then
1107 write (iout,*) "Cartesian and internal coordinates: step 2"
1109 call pdbout(0.0d0,"cipiszcze",iout)
1111 write (iout,*) "Accelerations"
1113 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1114 & (d_a(j,i+nres),j=1,3)
1116 write (iout,*) "Velocities, step 2"
1118 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1119 & (d_t(j,i+nres),j=1,3)
1123 c Compute the final RESPA step (increment velocities)
1124 c write (iout,*) "*********************** RESPA fin"
1126 c Compute the complete potential energy
1128 potEcomp(i)=energia_short(i)+energia_long(i)
1130 potE=potEcomp(0)-potEcomp(20)
1131 c potE=energia_short(0)+energia_long(0)
1133 c Calculate the kinetic and the total energy and the kinetic temperature
1136 c Couple the system to Berendsen bath if needed
1137 if (tbf .and. lang.eq.0) then
1140 kinetic_T=2.0d0/(dimen*Rb)*EK
1141 c Backup the coordinates, velocities, and accelerations
1143 if (mod(itime,ntwe).eq.0 .and. large) then
1144 write (iout,*) "Velocities, end"
1146 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1147 & (d_t(j,i+nres),j=1,3)
1153 c---------------------------------------------------------------------
1154 subroutine RESPA_vel
1155 c First and last RESPA step (incrementing velocities using long-range
1157 implicit real*8 (a-h,o-z)
1158 include 'DIMENSIONS'
1159 include 'COMMON.CONTROL'
1160 include 'COMMON.VAR'
1162 include 'COMMON.CHAIN'
1163 include 'COMMON.DERIV'
1164 include 'COMMON.GEO'
1165 include 'COMMON.LOCAL'
1166 include 'COMMON.INTERACT'
1167 include 'COMMON.IOUNITS'
1168 include 'COMMON.NAMES'
1170 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1174 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1178 if (itype(i).ne.10) then
1181 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1187 c-----------------------------------------------------------------
1189 c Applying velocity Verlet algorithm - step 1 to coordinates
1190 implicit real*8 (a-h,o-z)
1191 include 'DIMENSIONS'
1192 include 'COMMON.CONTROL'
1193 include 'COMMON.VAR'
1195 include 'COMMON.CHAIN'
1196 include 'COMMON.DERIV'
1197 include 'COMMON.GEO'
1198 include 'COMMON.LOCAL'
1199 include 'COMMON.INTERACT'
1200 include 'COMMON.IOUNITS'
1201 include 'COMMON.NAMES'
1202 double precision adt,adt2
1205 adt=d_a_old(j,0)*d_time
1207 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1208 d_t_new(j,0)=d_t_old(j,0)+adt2
1209 d_t(j,0)=d_t_old(j,0)+adt
1213 adt=d_a_old(j,i)*d_time
1215 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1216 d_t_new(j,i)=d_t_old(j,i)+adt2
1217 d_t(j,i)=d_t_old(j,i)+adt
1221 if (itype(i).ne.10) then
1224 adt=d_a_old(j,inres)*d_time
1226 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1227 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1228 d_t(j,inres)=d_t_old(j,inres)+adt
1234 c---------------------------------------------------------------------
1236 c Step 2 of the velocity Verlet algorithm: update velocities
1237 implicit real*8 (a-h,o-z)
1238 include 'DIMENSIONS'
1239 include 'COMMON.CONTROL'
1240 include 'COMMON.VAR'
1242 include 'COMMON.CHAIN'
1243 include 'COMMON.DERIV'
1244 include 'COMMON.GEO'
1245 include 'COMMON.LOCAL'
1246 include 'COMMON.INTERACT'
1247 include 'COMMON.IOUNITS'
1248 include 'COMMON.NAMES'
1250 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1254 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1258 if (itype(i).ne.10) then
1261 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1267 c-----------------------------------------------------------------
1268 subroutine sddir_precalc
1269 c Applying velocity Verlet algorithm - step 1 to coordinates
1270 implicit real*8 (a-h,o-z)
1271 include 'DIMENSIONS'
1272 include 'COMMON.CONTROL'
1273 include 'COMMON.VAR'
1276 include 'COMMON.LANGEVIN'
1278 include 'COMMON.LANGEVIN.lang0'
1280 include 'COMMON.CHAIN'
1281 include 'COMMON.DERIV'
1282 include 'COMMON.GEO'
1283 include 'COMMON.LOCAL'
1284 include 'COMMON.INTERACT'
1285 include 'COMMON.IOUNITS'
1286 include 'COMMON.NAMES'
1287 double precision stochforcvec(MAXRES6)
1288 common /stochcalc/ stochforcvec
1290 c Compute friction and stochastic forces
1293 call stochastic_force(stochforcvec)
1295 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1296 c forces (d_as_work)
1303 d_af_work(i)=d_af_work(i)+Ginv(i,j)*fric_work(j)
1304 d_as_work(i)=d_as_work(i)+Ginv(i,j)*stochforcvec(j)
1308 call ginv_mult(fric_work, d_af_work)
1309 call ginv_mult(stochforcvec, d_as_work)
1313 c---------------------------------------------------------------------
1314 subroutine sddir_verlet1
1315 c Applying velocity Verlet algorithm - step 1 to velocities
1316 implicit real*8 (a-h,o-z)
1317 include 'DIMENSIONS'
1318 include 'COMMON.CONTROL'
1319 include 'COMMON.VAR'
1322 include 'COMMON.LANGEVIN'
1324 include 'COMMON.LANGEVIN.lang0'
1326 include 'COMMON.CHAIN'
1327 include 'COMMON.DERIV'
1328 include 'COMMON.GEO'
1329 include 'COMMON.LOCAL'
1330 include 'COMMON.INTERACT'
1331 include 'COMMON.IOUNITS'
1332 include 'COMMON.NAMES'
1333 c Revised 3/31/05 AL: correlation between random contributions to
1334 c position and velocity increments included.
1335 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1336 double precision adt,adt2
1338 c Add the contribution from BOTH friction and stochastic force to the
1339 c coordinates, but ONLY the contribution from the friction forces to velocities
1342 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1343 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1344 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1345 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1346 d_t(j,0)=d_t_old(j,0)+adt
1351 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1352 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1353 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1354 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1355 d_t(j,i)=d_t_old(j,i)+adt
1360 if (itype(i).ne.10) then
1363 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1364 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1365 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1366 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1367 d_t(j,inres)=d_t_old(j,inres)+adt
1374 c---------------------------------------------------------------------
1375 subroutine sddir_verlet2
1376 c Calculating the adjusted velocities for accelerations
1377 implicit real*8 (a-h,o-z)
1378 include 'DIMENSIONS'
1379 include 'COMMON.CONTROL'
1380 include 'COMMON.VAR'
1383 include 'COMMON.LANGEVIN'
1385 include 'COMMON.LANGEVIN.lang0'
1387 include 'COMMON.CHAIN'
1388 include 'COMMON.DERIV'
1389 include 'COMMON.GEO'
1390 include 'COMMON.LOCAL'
1391 include 'COMMON.INTERACT'
1392 include 'COMMON.IOUNITS'
1393 include 'COMMON.NAMES'
1394 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1395 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1396 c Revised 3/31/05 AL: correlation between random contributions to
1397 c position and velocity increments included.
1398 c The correlation coefficients are calculated at low-friction limit.
1399 c Also, friction forces are now not calculated with new velocities.
1401 c call friction_force
1402 call stochastic_force(stochforcvec)
1404 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1405 c forces (d_as_work)
1409 c d_af_work(i)=0.0d0
1412 c d_af_work(i)=d_af_work(i)+Ginv(i,j)*fric_work(j)
1413 d_as_work1(i)=d_as_work1(i)+Ginv(i,j)*stochforcvec(j)
1417 call ginv_mult(stochforcvec, d_as_work1)
1424 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1425 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1430 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1431 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1436 if (itype(i).ne.10) then
1439 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1440 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1441 & +cos60*d_as_work1(ind+j))*d_time
1448 c---------------------------------------------------------------------
1449 subroutine max_accel
1451 c Find the maximum difference in the accelerations of the the sites
1452 c at the beginning and the end of the time step.
1454 implicit real*8 (a-h,o-z)
1455 include 'DIMENSIONS'
1456 include 'COMMON.CONTROL'
1457 include 'COMMON.VAR'
1459 include 'COMMON.CHAIN'
1460 include 'COMMON.DERIV'
1461 include 'COMMON.GEO'
1462 include 'COMMON.LOCAL'
1463 include 'COMMON.INTERACT'
1464 include 'COMMON.IOUNITS'
1465 double precision aux(3),accel(3)
1467 aux(j)=d_a(j,0)-d_a_old(j,0)
1474 accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1475 if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1482 if (itype(i).ne.10) then
1484 accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1488 if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1491 aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1496 c---------------------------------------------------------------------
1497 subroutine predict_edrift(epdrift)
1499 c Predict the drift of the potential energy
1501 implicit real*8 (a-h,o-z)
1502 include 'DIMENSIONS'
1503 include 'COMMON.CONTROL'
1504 include 'COMMON.VAR'
1506 include 'COMMON.CHAIN'
1507 include 'COMMON.DERIV'
1508 include 'COMMON.GEO'
1509 include 'COMMON.LOCAL'
1510 include 'COMMON.INTERACT'
1511 include 'COMMON.IOUNITS'
1512 include 'COMMON.MUCA'
1513 double precision epdrift,epdriftij
1514 c Drift of the potential energy
1520 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1521 if (lmuca) epdriftij=epdriftij*factor
1522 c write (iout,*) "back",i,j,epdriftij
1523 if (epdriftij.gt.epdrift) epdrift=epdriftij
1527 if (itype(i).ne.10) then
1530 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1531 if (lmuca) epdriftij=epdriftij*factor
1532 c write (iout,*) "side",i,j,epdriftij
1533 if (epdriftij.gt.epdrift) epdrift=epdriftij
1537 epdrift=0.5d0*epdrift*d_time*d_time
1538 c write (iout,*) "epdrift",epdrift
1541 c-----------------------------------------------------------------------
1542 subroutine verlet_bath
1544 c Coupling to the thermostat by using the Berendsen algorithm
1546 implicit real*8 (a-h,o-z)
1547 include 'DIMENSIONS'
1548 include 'COMMON.CONTROL'
1549 include 'COMMON.VAR'
1551 include 'COMMON.CHAIN'
1552 include 'COMMON.DERIV'
1553 include 'COMMON.GEO'
1554 include 'COMMON.LOCAL'
1555 include 'COMMON.INTERACT'
1556 include 'COMMON.IOUNITS'
1557 include 'COMMON.NAMES'
1558 double precision T_half,fact
1560 T_half=2.0d0/(dimen*Rb)*EK
1561 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1562 c write(iout,*) "T_half", T_half
1563 c write(iout,*) "EK", EK
1564 c write(iout,*) "fact", fact
1566 d_t(j,0)=fact*d_t(j,0)
1570 d_t(j,i)=fact*d_t(j,i)
1574 if (itype(i).ne.10) then
1577 d_t(j,inres)=fact*d_t(j,inres)
1583 c---------------------------------------------------------
1585 c Set up the initial conditions of a MD simulation
1586 implicit real*8 (a-h,o-z)
1587 include 'DIMENSIONS'
1590 include 'COMMON.SETUP'
1593 include 'COMMON.CONTROL'
1594 include 'COMMON.VAR'
1597 include 'COMMON.LANGEVIN'
1599 include 'COMMON.LANGEVIN.lang0'
1601 include 'COMMON.CHAIN'
1602 include 'COMMON.DERIV'
1603 include 'COMMON.GEO'
1604 include 'COMMON.LOCAL'
1605 include 'COMMON.INTERACT'
1606 include 'COMMON.IOUNITS'
1607 include 'COMMON.NAMES'
1608 include 'COMMON.REMD'
1609 real*8 energia_long(0:n_ene),
1610 & energia_short(0:n_ene),vcm(3),incr(3)
1611 double precision cm(3),L(3),xv,sigv,lowb,highb
1612 double precision varia(maxvar)
1620 c write(iout,*) "d_time", d_time
1621 c Compute the standard deviations of stochastic forces for Langevin dynamics
1622 c if the friction coefficients do not depend on surface area
1623 if (lang.gt.0 .and. .not.surfarea) then
1625 stdforcp(i)=stdfp*dsqrt(gamp)
1628 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
1631 c Open the pdb file for snapshotshots
1634 if (ilen(tmpdir).gt.0)
1635 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1636 & liczba(:ilen(liczba))//".pdb")
1638 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1642 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1643 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1644 & liczba(:ilen(liczba))//".x")
1645 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1648 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1649 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1650 & liczba(:ilen(liczba))//".cx")
1651 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1657 if (ilen(tmpdir).gt.0)
1658 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1659 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1661 if (ilen(tmpdir).gt.0)
1662 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1663 cartname=prefix(:ilen(prefix))//"_MD.cx"
1667 write (qstr,'(256(1h ))')
1670 iq = qinfrag(i,iset)*10
1671 iw = wfrag(i,iset)/100
1673 if(me.eq.king.or..not.out1file)
1674 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1675 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1680 iq = qinpair(i,iset)*10
1681 iw = wpair(i,iset)/100
1683 if(me.eq.king.or..not.out1file)
1684 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1685 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1689 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1691 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1693 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1695 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1699 if (restart1file) then
1701 & inquire(file=mremd_rst_name,exist=file_exist)
1702 write (*,*) me," Before broadcast: file_exist",file_exist
1703 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1705 write (*,*) me," After broadcast: file_exist",file_exist
1706 c inquire(file=mremd_rst_name,exist=file_exist)
1707 if(me.eq.king.or..not.out1file)
1708 & write(iout,*) "Initial state read by master and distributed"
1710 if (ilen(tmpdir).gt.0)
1711 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1712 & //liczba(:ilen(liczba))//'.rst')
1713 inquire(file=rest2name,exist=file_exist)
1716 if(.not.restart1file) then
1717 if(me.eq.king.or..not.out1file)
1718 & write(iout,*) "Initial state will be read from file ",
1719 & rest2name(:ilen(rest2name))
1722 call rescale_weights(t_bath)
1724 if(me.eq.king.or..not.out1file)then
1725 if (restart1file) then
1726 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1729 write(iout,*) "File ",rest2name(:ilen(rest2name)),
1732 write(iout,*) "Initial velocities randomly generated"
1738 c Generate initial velocities
1739 if(me.eq.king.or..not.out1file)
1740 & write(iout,*) "Initial velocities randomly generated"
1744 c rest2name = prefix(:ilen(prefix))//'.rst'
1745 if(me.eq.king.or..not.out1file)then
1746 write(iout,*) "Initial backbone velocities"
1748 write(iout,*) (d_t(j,i),j=1,3)
1750 write(iout,*) "Initial side-chain velocities"
1752 write(iout,*) (d_t(j,i+nres),j=1,3)
1754 c Zeroing the total angular momentum of the system
1755 write(iout,*) "Calling the zero-angular
1756 & momentum subroutine"
1759 c Getting the potential energy and forces and velocities and accelerations
1761 c write (iout,*) "velocity of the center of the mass:"
1762 c write (iout,*) (vcm(j),j=1,3)
1764 d_t(j,0)=d_t(j,0)-vcm(j)
1766 c Removing the velocity of the center of mass
1768 if(me.eq.king.or..not.out1file)then
1769 write (iout,*) "vcm right after adjustment:"
1770 write (iout,*) (vcm(j),j=1,3)
1774 if(iranconf.ne.0) then
1776 print *, 'Calling OVERLAP_SC'
1777 call overlap_sc(fail)
1781 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1782 print *,'SC_move',nft_sc,etot
1783 if(me.eq.king.or..not.out1file)
1784 & write(iout,*) 'SC_move',nft_sc,etot
1788 print *, 'Calling MINIM_DC'
1789 call minim_dc(etot,iretcode,nfun)
1791 call geom_to_var(nvar,varia)
1792 print *,'Calling MINIMIZE.'
1793 call minimize(etot,varia,iretcode,nfun)
1794 call var_to_geom(nvar,varia)
1796 if(me.eq.king.or..not.out1file)
1797 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1800 call chainbuild_cart
1803 call verlet_bath(EK)
1805 kinetic_T=2.0d0/(dimen*Rb)*EK
1806 if(me.eq.king.or..not.out1file)then
1811 call etotal(potEcomp)
1816 if (amax*d_time .gt. dvmax) d_time=d_time*dvmax/amax
1817 if(me.eq.king.or..not.out1file)then
1818 write(iout,*) "Potential energy and its components"
1819 write(iout,*) (potEcomp(i),i=0,n_ene)
1821 potE=potEcomp(0)-potEcomp(20)
1824 if (ntwe.ne.0) call statout(itime)
1825 if(me.eq.king.or..not.out1file)
1826 & write (iout,*) "Initial:",
1827 & " Kinetic energy",EK," potential energy",potE,
1828 & " total energy",totE," maximum acceleration ",
1831 write (iout,*) "Initial coordinates"
1833 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1834 & (c(j,i+nres),j=1,3)
1836 write (iout,*) "Initial dC"
1838 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1839 & (dc(j,i+nres),j=1,3)
1841 write (iout,*) "Initial velocities"
1843 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1844 & (d_t(j,i+nres),j=1,3)
1846 write (iout,*) "Initial accelerations"
1848 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1849 & (d_a(j,i+nres),j=1,3)
1855 d_t_old(j,i)=d_t(j,i)
1856 d_a_old(j,i)=d_a(j,i)
1858 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1867 call etotal_long(energia_long)
1868 if (large) write (iout,*) "energia_long",energia_long(0)
1872 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1874 t_enegrad=t_enegrad+tcpu()-tt0
1876 c call etotal_short(energia_short)
1877 c write (iout,*) "energia_long",energia_long(0),
1878 c & " energia_short",energia_short(0),
1879 c & " total",energia_long(0)+energia_short(0)
1883 c-----------------------------------------------------------
1884 subroutine random_vel
1885 implicit real*8 (a-h,o-z)
1886 include 'DIMENSIONS'
1887 include 'COMMON.CONTROL'
1888 include 'COMMON.VAR'
1891 include 'COMMON.LANGEVIN'
1893 include 'COMMON.LANGEVIN.lang0'
1895 include 'COMMON.CHAIN'
1896 include 'COMMON.DERIV'
1897 include 'COMMON.GEO'
1898 include 'COMMON.LOCAL'
1899 include 'COMMON.INTERACT'
1900 include 'COMMON.IOUNITS'
1901 include 'COMMON.NAMES'
1902 include 'COMMON.TIME1'
1903 double precision xv,sigv,lowb,highb
1904 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1905 c First generate velocities in the eigenspace of the G matrix
1906 c write (iout,*) "Calling random_vel"
1909 sigv=dsqrt((Rb*t_bath)/geigen(i))
1912 d_t_work_new(i)=anorm_distr(xv,sigv,lowb,highb)
1916 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(i)**2
1918 c Transform velocities to UNRES coordinate space
1922 d_t_work(i)=d_t_work(i)+Gvec(i,j)*d_t_work_new(j)
1925 c Transfer to the d_t vector
1927 d_t(j,0)=d_t_work(j)
1933 d_t(j,i)=d_t_work(ind)
1937 if (itype(i).ne.10) then
1940 d_t(j,i+nres)=d_t_work(ind)
1945 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1946 c & 2.0d0/(dimen*Rb)*EK,2.0d0/(dimen*Rb)*EK1
1949 c-----------------------------------------------------------
1950 subroutine sd_verlet_p_setup
1951 c Sets up the parameters of stochastic Verlet algorithm
1952 implicit real*8 (a-h,o-z)
1953 include 'DIMENSIONS'
1957 include 'COMMON.CONTROL'
1958 include 'COMMON.VAR'
1961 include 'COMMON.LANGEVIN'
1963 include 'COMMON.LANGEVIN.lang0'
1965 include 'COMMON.CHAIN'
1966 include 'COMMON.DERIV'
1967 include 'COMMON.GEO'
1968 include 'COMMON.LOCAL'
1969 include 'COMMON.INTERACT'
1970 include 'COMMON.IOUNITS'
1971 include 'COMMON.NAMES'
1972 include 'COMMON.TIME1'
1973 double precision emgdt(MAXRES6),
1974 & pterm,vterm,rho,rhoc,vsig,
1975 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1976 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1977 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1978 logical lprn /.false./
1979 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1980 double precision ktm
1987 c AL 8/17/04 Code adapted from tinker
1989 c Get the frictional and random terms for stochastic dynamics in the
1990 c eigenspace of mass-scaled UNRES friction matrix
1993 gdt = fricgam(i) * d_time
1995 c Stochastic dynamics reduces to simple MD for zero friction
1997 if (gdt .le. zero) then
1998 pfric_vec(i) = 1.0d0
1999 vfric_vec(i) = d_time
2000 afric_vec(i) = 0.5d0 * d_time * d_time
2001 prand_vec(i) = 0.0d0
2002 vrand_vec1(i) = 0.0d0
2003 vrand_vec2(i) = 0.0d0
2005 c Analytical expressions when friction coefficient is large
2008 if (gdt .ge. gdt_radius) then
2011 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2012 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2013 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2014 vterm = 1.0d0 - egdt**2
2015 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2017 c Use series expansions when friction coefficient is small
2028 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2029 & - gdt5/120.0d0 + gdt6/720.0d0
2030 & - gdt7/5040.0d0 + gdt8/40320.0d0
2031 & - gdt9/362880.0d0) / fricgam(i)**2
2032 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2033 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2034 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2035 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2036 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2037 & + 127.0d0*gdt9/90720.0d0
2038 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2039 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2040 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2041 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2042 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2043 & - 17.0d0*gdt2/1280.0d0
2044 & + 17.0d0*gdt3/6144.0d0
2045 & + 40967.0d0*gdt4/34406400.0d0
2046 & - 57203.0d0*gdt5/275251200.0d0
2047 & - 1429487.0d0*gdt6/13212057600.0d0)
2050 c Compute the scaling factors of random terms for the nonzero friction case
2052 ktm = 0.5d0*d_time/fricgam(i)
2053 psig = dsqrt(ktm*pterm) / fricgam(i)
2054 vsig = dsqrt(ktm*vterm)
2055 rhoc = dsqrt(1.0d0 - rho*rho)
2057 vrand_vec1(i) = vsig * rho
2058 vrand_vec2(i) = vsig * rhoc
2063 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2066 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2067 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2071 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2074 call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
2075 call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
2076 call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
2077 call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
2078 call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec1,vrand_mat1)
2079 call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
2081 c call eigtransf1(dimen,maxres6,mt3mt2,pfric_vec,pfric_mat)
2082 c call eigtransf1(dimen,maxres6,mt3mt2,vfric_vec,vfric_mat)
2083 c call eigtransf1(dimen,maxres6,mt3mt2,afric_vec,afric_mat)
2084 c call eigtransf1(dimen,maxres6,mt3mt1,prand_vec,prand_mat)
2085 c call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec1,vrand_mat1)
2086 c call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec2,vrand_mat2)
2088 t_sdsetup=t_sdsetup+MPI_Wtime()
2090 t_sdsetup=t_sdsetup+tcpu()-tt0
2094 c-------------------------------------------------------------
2095 subroutine eigtransf1(n,ndim,ab,d,c)
2098 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2104 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2110 c-------------------------------------------------------------
2111 subroutine eigtransf(n,ndim,a,b,d,c)
2114 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2120 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2126 c-------------------------------------------------------------
2127 subroutine sd_verlet1
2128 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2129 implicit real*8 (a-h,o-z)
2130 include 'DIMENSIONS'
2131 include 'COMMON.CONTROL'
2132 include 'COMMON.VAR'
2135 include 'COMMON.LANGEVIN'
2137 include 'COMMON.LANGEVIN.lang0'
2139 include 'COMMON.CHAIN'
2140 include 'COMMON.DERIV'
2141 include 'COMMON.GEO'
2142 include 'COMMON.LOCAL'
2143 include 'COMMON.INTERACT'
2144 include 'COMMON.IOUNITS'
2145 include 'COMMON.NAMES'
2146 double precision stochforcvec(MAXRES6)
2147 common /stochcalc/ stochforcvec
2148 logical lprn /.false./
2150 c write (iout,*) "dc_old"
2152 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2153 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2156 dc_work(j)=dc_old(j,0)
2157 d_t_work(j)=d_t_old(j,0)
2158 d_a_work(j)=d_a_old(j,0)
2163 dc_work(ind+j)=dc_old(j,i)
2164 d_t_work(ind+j)=d_t_old(j,i)
2165 d_a_work(ind+j)=d_a_old(j,i)
2170 if (itype(i).ne.10) then
2172 dc_work(ind+j)=dc_old(j,i+nres)
2173 d_t_work(ind+j)=d_t_old(j,i+nres)
2174 d_a_work(ind+j)=d_a_old(j,i+nres)
2182 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2186 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2187 & vfric_mat(i,j),afric_mat(i,j),
2188 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2196 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2197 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2198 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2199 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2201 d_t_work_new(i)=ddt1+0.5d0*ddt2
2202 d_t_work(i)=ddt1+ddt2
2207 d_t(j,0)=d_t_work(j)
2212 dc(j,i)=dc_work(ind+j)
2213 d_t(j,i)=d_t_work(ind+j)
2218 if (itype(i).ne.10) then
2221 dc(j,inres)=dc_work(ind+j)
2222 d_t(j,inres)=d_t_work(ind+j)
2229 c--------------------------------------------------------------------------
2230 subroutine sd_verlet2
2231 c Calculating the adjusted velocities for accelerations
2232 implicit real*8 (a-h,o-z)
2233 include 'DIMENSIONS'
2234 include 'COMMON.CONTROL'
2235 include 'COMMON.VAR'
2238 include 'COMMON.LANGEVIN'
2240 include 'COMMON.LANGEVIN.lang0'
2242 include 'COMMON.CHAIN'
2243 include 'COMMON.DERIV'
2244 include 'COMMON.GEO'
2245 include 'COMMON.LOCAL'
2246 include 'COMMON.INTERACT'
2247 include 'COMMON.IOUNITS'
2248 include 'COMMON.NAMES'
2249 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2250 common /stochcalc/ stochforcvec
2252 c Compute the stochastic forces which contribute to velocity change
2254 call stochastic_force(stochforcvecV)
2261 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2262 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2263 & vrand_mat2(i,j)*stochforcvecV(j)
2265 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2269 d_t(j,0)=d_t_work(j)
2274 d_t(j,i)=d_t_work(ind+j)
2279 if (itype(i).ne.10) then
2282 d_t(j,inres)=d_t_work(ind+j)
2289 c-----------------------------------------------------------
2290 subroutine sd_verlet_ciccotti_setup
2291 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2293 implicit real*8 (a-h,o-z)
2294 include 'DIMENSIONS'
2298 include 'COMMON.CONTROL'
2299 include 'COMMON.VAR'
2302 include 'COMMON.LANGEVIN'
2304 include 'COMMON.LANGEVIN.lang0'
2306 include 'COMMON.CHAIN'
2307 include 'COMMON.DERIV'
2308 include 'COMMON.GEO'
2309 include 'COMMON.LOCAL'
2310 include 'COMMON.INTERACT'
2311 include 'COMMON.IOUNITS'
2312 include 'COMMON.NAMES'
2313 include 'COMMON.TIME1'
2314 double precision emgdt(MAXRES6),
2315 & pterm,vterm,rho,rhoc,vsig,
2316 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2317 & afric_vec(MAXRES6),prand_vec(MAXRES6),
2318 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2319 logical lprn /.false./
2320 double precision zero /1.0d-8/, gdt_radius /0.05d0/
2321 double precision ktm
2328 c AL 8/17/04 Code adapted from tinker
2330 c Get the frictional and random terms for stochastic dynamics in the
2331 c eigenspace of mass-scaled UNRES friction matrix
2334 write (iout,*) "i",i," fricgam",fricgam(i)
2335 gdt = fricgam(i) * d_time
2337 c Stochastic dynamics reduces to simple MD for zero friction
2339 if (gdt .le. zero) then
2340 pfric_vec(i) = 1.0d0
2341 vfric_vec(i) = d_time
2342 afric_vec(i) = 0.5d0*d_time*d_time
2343 prand_vec(i) = afric_vec(i)
2344 vrand_vec2(i) = vfric_vec(i)
2346 c Analytical expressions when friction coefficient is large
2351 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2352 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2353 prand_vec(i) = afric_vec(i)
2354 vrand_vec2(i) = vfric_vec(i)
2356 c Compute the scaling factors of random terms for the nonzero friction case
2358 c ktm = 0.5d0*d_time/fricgam(i)
2359 c psig = dsqrt(ktm*pterm) / fricgam(i)
2360 c vsig = dsqrt(ktm*vterm)
2361 c prand_vec(i) = psig*afric_vec(i)
2362 c vrand_vec2(i) = vsig*vfric_vec(i)
2367 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2370 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2371 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2375 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2377 call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
2378 call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
2379 call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
2380 call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
2381 call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
2383 t_sdsetup=t_sdsetup+MPI_Wtime()
2385 t_sdsetup=t_sdsetup+tcpu()-tt0
2389 c-------------------------------------------------------------
2390 subroutine sd_verlet1_ciccotti
2391 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2392 implicit real*8 (a-h,o-z)
2393 include 'DIMENSIONS'
2397 include 'COMMON.CONTROL'
2398 include 'COMMON.VAR'
2401 include 'COMMON.LANGEVIN'
2403 include 'COMMON.LANGEVIN.lang0'
2405 include 'COMMON.CHAIN'
2406 include 'COMMON.DERIV'
2407 include 'COMMON.GEO'
2408 include 'COMMON.LOCAL'
2409 include 'COMMON.INTERACT'
2410 include 'COMMON.IOUNITS'
2411 include 'COMMON.NAMES'
2412 double precision stochforcvec(MAXRES6)
2413 common /stochcalc/ stochforcvec
2414 logical lprn /.false./
2416 c write (iout,*) "dc_old"
2418 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2419 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2422 dc_work(j)=dc_old(j,0)
2423 d_t_work(j)=d_t_old(j,0)
2424 d_a_work(j)=d_a_old(j,0)
2429 dc_work(ind+j)=dc_old(j,i)
2430 d_t_work(ind+j)=d_t_old(j,i)
2431 d_a_work(ind+j)=d_a_old(j,i)
2436 if (itype(i).ne.10) then
2438 dc_work(ind+j)=dc_old(j,i+nres)
2439 d_t_work(ind+j)=d_t_old(j,i+nres)
2440 d_a_work(ind+j)=d_a_old(j,i+nres)
2449 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2453 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2454 & vfric_mat(i,j),afric_mat(i,j),
2455 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2463 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2464 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2465 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2466 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2468 d_t_work_new(i)=ddt1+0.5d0*ddt2
2469 d_t_work(i)=ddt1+ddt2
2474 d_t(j,0)=d_t_work(j)
2479 dc(j,i)=dc_work(ind+j)
2480 d_t(j,i)=d_t_work(ind+j)
2485 if (itype(i).ne.10) then
2488 dc(j,inres)=dc_work(ind+j)
2489 d_t(j,inres)=d_t_work(ind+j)
2496 c--------------------------------------------------------------------------
2497 subroutine sd_verlet2_ciccotti
2498 c Calculating the adjusted velocities for accelerations
2499 implicit real*8 (a-h,o-z)
2500 include 'DIMENSIONS'
2501 include 'COMMON.CONTROL'
2502 include 'COMMON.VAR'
2505 include 'COMMON.LANGEVIN'
2507 include 'COMMON.LANGEVIN.lang0'
2509 include 'COMMON.CHAIN'
2510 include 'COMMON.DERIV'
2511 include 'COMMON.GEO'
2512 include 'COMMON.LOCAL'
2513 include 'COMMON.INTERACT'
2514 include 'COMMON.IOUNITS'
2515 include 'COMMON.NAMES'
2516 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2517 common /stochcalc/ stochforcvec
2519 c Compute the stochastic forces which contribute to velocity change
2521 call stochastic_force(stochforcvecV)
2528 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2529 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2530 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2532 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2536 d_t(j,0)=d_t_work(j)
2541 d_t(j,i)=d_t_work(ind+j)
2546 if (itype(i).ne.10) then
2549 d_t(j,inres)=d_t_work(ind+j)