1 c---------------------------------------------------------
3 c Set up the initial conditions of a MD simulation
4 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'
28 real*8 energia_long(0:n_ene),
29 & energia_short(0:n_ene),vcm(3),incr(3)
30 double precision cm(3),L(3),xv,sigv,lowb,highb
31 double precision varia(maxvar)
39 c write(iout,*) "d_time", d_time
40 c Compute the standard deviations of stochastic forces for Langevin dynamics
41 c if the friction coefficients do not depend on surface area
42 if (lang.gt.0 .and. .not.surfarea) then
44 stdforcp(i)=stdfp*dsqrt(gamp)
47 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
50 c Open the pdb file for snapshotshots
53 if (ilen(tmpdir).gt.0)
54 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
55 & liczba(:ilen(liczba))//".pdb")
57 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
61 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
62 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
63 & liczba(:ilen(liczba))//".x")
64 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
67 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
68 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
69 & liczba(:ilen(liczba))//".cx")
70 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
76 if (ilen(tmpdir).gt.0)
77 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
78 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
80 if (ilen(tmpdir).gt.0)
81 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
82 cartname=prefix(:ilen(prefix))//"_MD.cx"
86 write (qstr,'(256(1h ))')
89 iq = qinfrag(i,iset)*10
90 iw = wfrag(i,iset)/100
92 if(me.eq.king.or..not.out1file)
93 & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
94 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
99 iq = qinpair(i,iset)*10
100 iw = wpair(i,iset)/100
102 if(me.eq.king.or..not.out1file)
103 & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
104 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
108 c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
110 c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
112 c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
114 c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
118 if (restart1file) then
120 & inquire(file=mremd_rst_name,exist=file_exist)
121 write (*,*) me," Before broadcast: file_exist",file_exist
122 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
124 write (*,*) me," After broadcast: file_exist",file_exist
125 c inquire(file=mremd_rst_name,exist=file_exist)
126 if(me.eq.king.or..not.out1file)
127 & write(iout,*) "Initial state read by master and distributed"
129 if (ilen(tmpdir).gt.0)
130 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
131 & //liczba(:ilen(liczba))//'.rst')
132 inquire(file=rest2name,exist=file_exist)
135 if(.not.restart1file) then
136 if(me.eq.king.or..not.out1file)
137 & write(iout,*) "Initial state will be read from file ",
138 & rest2name(:ilen(rest2name))
141 call rescale_weights(t_bath)
143 if(me.eq.king.or..not.out1file)then
144 if (restart1file) then
145 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
148 write(iout,*) "File ",rest2name(:ilen(rest2name)),
151 write(iout,*) "Initial velocities randomly generated"
157 c Generate initial velocities
158 if(me.eq.king.or..not.out1file)
159 & write(iout,*) "Initial velocities randomly generated"
163 c rest2name = prefix(:ilen(prefix))//'.rst'
164 if(me.eq.king.or..not.out1file)then
165 write (iout,*) "Initial velocities"
167 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
168 & (d_t(j,i+nres),j=1,3)
170 c Zeroing the total angular momentum of the system
171 write(iout,*) "Calling the zero-angular
172 & momentum subroutine"
175 c Getting the potential energy and forces and velocities and accelerations
177 c write (iout,*) "velocity of the center of the mass:"
178 c write (iout,*) (vcm(j),j=1,3)
180 d_t(j,0)=d_t(j,0)-vcm(j)
182 c Removing the velocity of the center of mass
184 if(me.eq.king.or..not.out1file)then
185 write (iout,*) "vcm right after adjustment:"
186 write (iout,*) (vcm(j),j=1,3)
190 if(iranconf.ne.0) then
192 print *, 'Calling OVERLAP_SC'
193 call overlap_sc(fail)
197 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
198 print *,'SC_move',nft_sc,etot
199 if(me.eq.king.or..not.out1file)
200 & write(iout,*) 'SC_move',nft_sc,etot
204 print *, 'Calling MINIM_DC'
205 call minim_dc(etot,iretcode,nfun)
207 call geom_to_var(nvar,varia)
208 print *,'Calling MINIMIZE.'
209 call minimize(etot,varia,iretcode,nfun)
210 call var_to_geom(nvar,varia)
212 if(me.eq.king.or..not.out1file)
213 & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
221 kinetic_T=2.0d0/(dimen3*Rb)*EK
222 if(me.eq.king.or..not.out1file)then
232 call etotal(potEcomp)
235 t_etotal=t_etotal+MPI_Wtime()-tt0
237 t_etotal=t_etotal+tcpu()-tt0
244 if (amax*d_time .gt. dvmax) then
245 d_time=d_time*dvmax/amax
246 if(me.eq.king.or..not.out1file) write (iout,*)
247 & "Time step reduced to",d_time,
248 & " because of too large initial acceleration."
250 if(me.eq.king.or..not.out1file)then
251 write(iout,*) "Potential energy and its components"
252 call enerprint(potEcomp)
253 c write(iout,*) (potEcomp(i),i=0,n_ene)
255 potE=potEcomp(0)-potEcomp(20)
258 if (ntwe.ne.0) call statout(itime)
259 if(me.eq.king.or..not.out1file)
260 & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
261 & " Kinetic energy",EK," potential energy",potE,
262 & " total energy",totE," maximum acceleration ",
265 write (iout,*) "Initial coordinates"
267 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
268 & (c(j,i+nres),j=1,3)
270 write (iout,*) "Initial dC"
272 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
273 & (dc(j,i+nres),j=1,3)
275 write (iout,*) "Initial velocities"
276 write (iout,"(13x,' backbone ',23x,' side chain')")
278 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
279 & (d_t(j,i+nres),j=1,3)
281 write (iout,*) "Initial accelerations"
283 c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
284 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
285 & (d_a(j,i+nres),j=1,3)
291 d_t_old(j,i)=d_t(j,i)
292 d_a_old(j,i)=d_a(j,i)
294 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
303 call etotal_short(energia_short)
306 t_eshort=t_eshort+MPI_Wtime()-tt0
308 t_eshort=t_eshort+tcpu()-tt0
313 if(.not.out1file .and. large) then
314 write (iout,*) "energia_long",energia_long(0),
315 & " energia_short",energia_short(0),
316 & " total",energia_long(0)+energia_short(0)
317 write (iout,*) "Initial fast-force accelerations"
319 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
320 & (d_a(j,i+nres),j=1,3)
323 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
326 d_a_short(j,i)=d_a(j,i)
335 call etotal_long(energia_long)
338 t_elong=t_elong+MPI_Wtime()-tt0
340 t_elong=t_elong+tcpu()-tt0
345 if(.not.out1file .and. large) then
346 write (iout,*) "energia_long",energia_long(0)
347 write (iout,*) "Initial slow-force accelerations"
349 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
350 & (d_a(j,i+nres),j=1,3)
354 t_enegrad=t_enegrad+MPI_Wtime()-tt0
356 t_enegrad=t_enegrad+tcpu()-tt0