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)
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
148 if (amax*d_time .gt. dvmax) d_time=d_time*dvmax/amax
149 write(iout,*) "Potential energy"
150 write(iout,*) (energia(i),i=0,n_ene)
151 potE=energia(0)-energia(20)
154 write (iout,*) "Initial:",
155 & " Kinetic energy",EK," potential energy",potE,
156 & " total energy",totE," maximum acceleration ",
159 write (iout,*) "Initial velocities"
161 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
162 & (d_t(j,i+nres),j=1,3)
164 write (iout,*) "Initial accelerations"
166 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
167 & (d_a(j,i+nres),j=1,3)
173 d_t_old(j,i)=d_t(j,i)
174 d_a_old(j,i)=d_a(j,i)
176 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
180 call etotal_long(energia_long)
183 c call etotal_short(energia_short)
184 c write (iout,*) "energia_long",energia_long(0),
185 c & " energia_short",energia_short(0),
186 c & " total",energia_long(0)+energia_short(0)