1 c---------------------------------------------------------
3 c Set up the initial conditions of a MD simulation
4 implicit real*8 (a-h,o-z)
13 include 'COMMON.CONTROL'
16 include 'COMMON.LANGEVIN'
17 include 'COMMON.CHAIN'
18 include 'COMMON.DERIV'
20 include 'COMMON.LOCAL'
21 include 'COMMON.INTERACT'
22 include 'COMMON.IOUNITS'
23 include 'COMMON.NAMES'
24 real*8 energia(0:n_ene),energia_long(0:n_ene),
25 & energia_short(0:n_ene),vcm(3),incr(3),E_short
26 double precision cm(3),L(3),xv,sigv,lowb,highb
33 write(iout,*) "d_time", d_time
34 c Compute the standard deviations of stochastic forces for Langevin dynamics
35 c if the friction coefficients do not depend on surface area
36 if (lang.gt.0 .and. .not.surfarea) then
38 stdforcp(i)=stdfp*dsqrt(gamp)
41 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
44 c Open the pdb file for snapshotshots
49 npos = dlog10(dfloat(nprocs-1))+1
52 write (liczba,'(i1)') npos
53 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
55 write (liczba,form) myrank
58 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
61 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
66 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
68 cartname=prefix(:ilen(prefix))//"_MD.cx"
72 write (qstr,'(256(1h ))')
78 write (iout,*) "Frag",qinfrag(i),wfrag(i),iq,iw
79 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
87 write (iout,*) "Pair",i,qinpair(i),wpair(i),iq,iw
88 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
92 pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
93 cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
94 statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
98 write(iout,*) "Initial state will be read from file ",
99 & rest2name(:ilen(rest2name))
102 c Generate initial velocities
103 write(iout,*) "Initial velocities randomly generated"
107 c rest2name = prefix(:ilen(prefix))//'.rst'
108 write(iout,*) "Initial backbone velocities"
110 write(iout,*) (d_t(j,i),j=1,3)
112 write(iout,*) "Initial side-chain velocities"
114 write(iout,*) (d_t(j,i+nres),j=1,3)
116 c Zeroing the total angular momentum of the system
117 write(iout,*) "Calling the zero-angular
118 & momentum subroutine"
120 c Getting the potential energy and forces and velocities and accelerations
122 c write (iout,*) "velocity of the center of the mass:"
123 c write (iout,*) (vcm(j),j=1,3)
125 d_t(j,0)=d_t(j,0)-vcm(j)
127 c Removing the velocity of the center of mass
129 write (iout,*) "vcm right after adjustment:"
130 write (iout,*) (vcm(j),j=1,3)
139 kinetic_T=2.0d0/(dimen*Rb)*EK
146 if(tnp .or. tnp1) then
149 HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
151 write(iout,*) 'H0= ',H0
155 HNose1=Hnose_nh(EK,potE)
157 write (iout,*) 'H0= ',H0
164 if (amax*d_time .gt. dvmax) d_time=d_time*dvmax/amax
165 write(iout,*) "Potential energy"
166 write(iout,*) (energia(i),i=0,n_ene)
167 potE=energia(0)-energia(20)
171 write (iout,*) "Initial:",
172 & " Kinetic energy",EK," potential energy",potE,
173 & " total energy",totE," maximum acceleration ",
176 write (iout,*) "Initial velocities"
178 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
179 & (d_t(j,i+nres),j=1,3)
181 write (iout,*) "Initial accelerations"
183 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
184 & (d_a(j,i+nres),j=1,3)
190 d_t_old(j,i)=d_t(j,i)
191 d_a_old(j,i)=d_a(j,i)
193 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
197 call etotal_long(energia_long)
198 E_long=energia_long(0)
199 if(tnp .or. tnp1) then
200 call etotal_short(energia_short)
201 E_short=energia_short(0)
202 HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen)
205 c_new_var_csplit Csplit=H0-E_long
206 c Csplit = H0-energia_short(0)
207 write(iout,*) 'Csplit= ',Csplit
211 c call etotal_short(energia_short)
212 c write (iout,*) "energia_long",energia_long(0),
213 c & " energia_short",energia_short(0),
214 c & " total",energia_long(0)+energia_short(0)