1 c-------------------------------------------------------------------------------
2 subroutine brown_step(itime)
3 c------------------------------------------------
4 c Perform a single Euler integration step of Brownian dynamics
5 c------------------------------------------------
6 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
11 include 'COMMON.LANGEVIN'
12 include 'COMMON.CHAIN'
13 include 'COMMON.DERIV'
15 include 'COMMON.LOCAL'
16 include 'COMMON.INTERACT'
17 include 'COMMON.IOUNITS'
18 include 'COMMON.NAMES'
19 include 'COMMON.TIME1'
20 double precision energia(0:n_ene),zapas(MAXRES6)
23 double precision stochforcvec(MAXRES6)
24 double precision Bmat(MAXRES6,MAXRES2),Cmat(maxres2,maxres2),
25 & Cinv(maxres2,maxres2),GBmat(MAXRES6,MAXRES2),
26 & Tmat(MAXRES6,MAXRES2),Pmat(maxres6,maxres6),Td(maxres6),
28 common /stochcalc/ stochforcvec
29 common /gucio/ cm, energia
31 logical lprn /.false./,lprn1 /.false./
33 double precision difftol /1.0d-5/
36 if (itype(i).ne.10) nbond=nbond+1
40 write (iout,*) "Generalized inverse of fricmat"
41 call matout(dimen,dimen,MAXRES6,MAXRES6,fricmat)
53 Bmat(ind+j,ind1)=dC_norm(j,i)
58 if (itype(i).ne.10) then
61 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
67 write (iout,*) "Matrix Bmat"
68 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Bmat)
74 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
79 write (iout,*) "Matrix GBmat"
80 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Gbmat)
86 Cmat(i,j)=Cmat(i,j)+Bmat(k,i)*GBmat(k,j)
91 write (iout,*) "Matrix Cmat"
92 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
94 call matinvert(nbond,MAXRES2,Cmat,Cinv)
96 write (iout,*) "Matrix Cinv"
97 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cinv)
103 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
108 write (iout,*) "Matrix Tmat"
109 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Tmat)
119 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
124 write (iout,*) "Matrix Pmat"
125 call MATOUT(dimen,dimen,MAXRES6,MAXRES6,Pmat)
132 Td(i)=Td(i)+vbl*Tmat(i,ind)
135 if (itype(k).ne.10) then
137 Td(i)=Td(i)+vbldsc0(itype(k))*Tmat(i,ind)
142 write (iout,*) "Vector Td"
144 write (iout,'(i5,f10.5)') i,Td(i)
147 call stochastic_force(stochforcvec)
149 write (iout,*) "stochforcvec"
151 write (iout,*) i,stochforcvec(i)
155 zapas(j)=-gcart(j,0)+stochforcvec(j)
157 dC_work(j)=dC_old(j,0)
163 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
164 dC_work(ind)=dC_old(j,i)
168 if (itype(i).ne.10) then
171 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
172 dC_work(ind)=dC_old(j,i+nres)
178 write (iout,*) "Initial d_t_work"
180 write (iout,*) i,d_t_work(i)
187 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
194 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
198 write (iout,*) "Final d_t_work and zapas"
200 write (iout,*) i,d_t_work(i),zapas(i)
214 dc_work(ind+j)=dc(j,i)
220 d_t(j,i+nres)=d_t_work(ind+j)
221 dc(j,i+nres)=zapas(ind+j)
222 dc_work(ind+j)=dc(j,i+nres)
228 write (iout,*) "Before correction for rotational lengthening"
229 write (iout,*) "New coordinates",
230 & " and differences between actual and standard bond lengths"
235 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
236 & i,(dC(j,i),j=1,3),xx
239 if (itype(i).ne.10) then
241 xx=vbld(i+nres)-vbldsc0(itype(i))
242 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
243 & i,(dC(j,i+nres),j=1,3),xx
247 c Second correction (rotational lengthening)
253 blen2 = scalar(dc(1,i),dc(1,i))
254 ppvec(ind)=2*vbl**2-blen2
255 diffbond=dabs(vbl-dsqrt(blen2))
256 if (diffbond.gt.diffmax) diffmax=diffbond
257 if (ppvec(ind).gt.0.0d0) then
258 ppvec(ind)=dsqrt(ppvec(ind))
263 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
267 if (itype(i).ne.10) then
269 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
270 ppvec(ind)=2*vbldsc0(itype(i))**2-blen2
271 diffbond=dabs(vbldsc0(itype(i))-dsqrt(blen2))
272 if (diffbond.gt.diffmax) diffmax=diffbond
273 if (ppvec(ind).gt.0.0d0) then
274 ppvec(ind)=dsqrt(ppvec(ind))
279 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
283 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
284 if (diffmax.lt.difftol) goto 10
288 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
294 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
305 dc_work(ind+j)=zapas(ind+j)
310 if (itype(i).ne.10) then
312 dc(j,i+nres)=zapas(ind+j)
313 dc_work(ind+j)=zapas(ind+j)
318 c Building the chain from the newly calculated coordinates
320 if (large.and. mod(itime,ntwe).eq.0) then
321 write (iout,*) "Cartesian and internal coordinates: step 1"
324 write (iout,'(a)') "Potential forces"
326 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),
327 & (-gxcart(j,i),j=1,3)
329 write (iout,'(a)') "Stochastic forces"
331 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),
332 & (stochforc(j,i+nres),j=1,3)
334 write (iout,'(a)') "Velocities"
336 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
337 & (d_t(j,i+nres),j=1,3)
341 write (iout,*) "After correction for rotational lengthening"
342 write (iout,*) "New coordinates",
343 & " and differences between actual and standard bond lengths"
348 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
349 & i,(dC(j,i),j=1,3),xx
352 if (itype(i).ne.10) then
354 xx=vbld(i+nres)-vbldsc0(itype(i))
355 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
356 & i,(dC(j,i+nres),j=1,3),xx
361 c write (iout,*) "Too many attempts at correcting the bonds"
365 c Calculate energy and forces
368 potE=energia(0)-energia(20)
371 c Calculate the kinetic and total energy and the kinetic temperature
373 t_enegrad=t_enegrad+tcpu()-tt0
375 kinetic_T=2.0d0/(dimen*Rb)*EK