2 c------------------------------------------------
3 c The driver for molecular dynamics subroutines
4 c------------------------------------------------
5 implicit real*8 (a-h,o-z)
11 include 'COMMON.SETUP'
12 include 'COMMON.CONTROL'
16 include 'COMMON.LANGEVIN'
18 include 'COMMON.LANGEVIN.lang0'
20 include 'COMMON.CHAIN'
21 include 'COMMON.DERIV'
23 include 'COMMON.LOCAL'
24 include 'COMMON.INTERACT'
25 include 'COMMON.IOUNITS'
26 include 'COMMON.NAMES'
27 include 'COMMON.TIME1'
28 double precision cm(3),L(3),vcm(3)
30 double precision v_work(maxres6),v_transf(maxres6)
39 if (ilen(tmpdir).gt.0)
40 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
41 & //liczba(:ilen(liczba))//'.rst')
43 if (ilen(tmpdir).gt.0)
44 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
51 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
57 c Determine the inverse of the inertia matrix.
58 call setup_MD_matrices
62 t_MDsetup = MPI_Wtime()-tt0
64 t_MDsetup = tcpu()-tt0
67 c Entering the MD loop
73 if (lang.eq.2 .or. lang.eq.3) then
77 call sd_verlet_p_setup
79 call sd_verlet_ciccotti_setup
83 pfric0_mat(i,j,0)=pfric_mat(i,j)
84 afric0_mat(i,j,0)=afric_mat(i,j)
85 vfric0_mat(i,j,0)=vfric_mat(i,j)
86 prand0_mat(i,j,0)=prand_mat(i,j)
87 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
88 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
97 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
99 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
103 else if (lang.eq.1 .or. lang.eq.4) then
107 t_langsetup=MPI_Wtime()-tt0
110 t_langsetup=tcpu()-tt0
113 do itime=1,n_timestep
115 if (lang.gt.0 .and. surfarea .and.
116 & mod(itime,reset_fricmat).eq.0) then
117 if (lang.eq.2 .or. lang.eq.3) then
121 call sd_verlet_p_setup
123 call sd_verlet_ciccotti_setup
127 pfric0_mat(i,j,0)=pfric_mat(i,j)
128 afric0_mat(i,j,0)=afric_mat(i,j)
129 vfric0_mat(i,j,0)=vfric_mat(i,j)
130 prand0_mat(i,j,0)=prand_mat(i,j)
131 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
132 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
137 flag_stoch(i)=.false.
140 else if (lang.eq.1 .or. lang.eq.4) then
143 write (iout,'(a,i10)')
144 & "Friction matrix reset based on surface area, itime",itime
146 if (reset_vel .and. tbf .and. lang.eq.0
147 & .and. mod(itime,count_reset_vel).eq.0) then
149 write(iout,'(a,f20.2)')
150 & "Velocities reset to random values, time",totT
153 d_t_old(j,i)=d_t(j,i)
157 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
161 d_t(j,0)=d_t(j,0)-vcm(j)
164 kinetic_T=2.0d0/(dimen3*Rb)*EK
165 scalfac=dsqrt(T_bath/kinetic_T)
166 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
169 d_t_old(j,i)=scalfac*d_t(j,i)
175 c Time-reversible RESPA algorithm
176 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
177 call RESPA_step(itime)
179 c Variable time step algorithm.
180 call velverlet_step(itime)
184 call brown_step(itime)
186 print *,"Brown dynamics not here!"
188 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
194 if (mod(itime,ntwe).eq.0) call statout(itime)
207 if (itype(i).ne.10) then
210 v_work(ind)=d_t(j,i+nres)
215 write (66,'(80f10.5)')
216 & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
220 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
222 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
224 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
227 if (mod(itime,ntwx).eq.0) then
228 write (tytul,'("time",f8.2)') totT
230 call pdbout(potE,tytul,ipdb)
235 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
236 open(irest2,file=rest2name,status='unknown')
237 write(irest2,*) totT,EK,potE,totE,t_bath
239 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
242 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
253 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
255 & 'MD calculations setup:',t_MDsetup,
256 & 'Energy & gradient evaluation:',t_enegrad,
257 & 'Stochastic MD setup:',t_langsetup,
258 & 'Stochastic MD step setup:',t_sdsetup,
260 write (iout,'(/28(1h=),a25,27(1h=))')
261 & ' End of MD calculation '
263 write (iout,*) "time for etotal",t_etotal," elong",t_elong,
265 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
266 & " time_fricmatmult",time_fricmatmult," time_fsample ",