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)
11 include 'COMMON.CONTROL'
15 include 'COMMON.LAGRANGE.5diag'
17 include 'COMMON.LAGRANGE'
20 include 'COMMON.LANGEVIN'
23 include 'COMMON.LANGEVIN.lang0.5diag'
25 include 'COMMON.LANGEVIN.lang0'
28 include 'COMMON.CHAIN'
29 include 'COMMON.DERIV'
31 include 'COMMON.LOCAL'
32 include 'COMMON.INTERACT'
33 include 'COMMON.IOUNITS'
34 include 'COMMON.NAMES'
35 include 'COMMON.TIME1'
36 double precision zapas(MAXRES6)
39 double precision stochforcvec(MAXRES6)
40 double precision Bmat(MAXRES6,MAXRES2),Cmat(maxres2,maxres2),
41 & Cinv(maxres2,maxres2),GBmat(MAXRES6,MAXRES2),
42 & Tmat(MAXRES6,MAXRES2),Pmat(maxres6,maxres6),Td(maxres6),
44 common /stochcalc/ stochforcvec
47 logical lprn /.false./,lprn1 /.false./
49 double precision difftol /1.0d-5/
52 if (itype(i).ne.10) nbond=nbond+1
56 write (iout,*) "Generalized inverse of fricmat"
57 call matout(dimen,dimen,MAXRES6,MAXRES6,fricmat)
69 Bmat(ind+j,ind1)=dC_norm(j,i)
74 if (itype(i).ne.10) then
77 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
83 write (iout,*) "Matrix Bmat"
84 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Bmat)
90 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
95 write (iout,*) "Matrix GBmat"
96 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Gbmat)
102 Cmat(i,j)=Cmat(i,j)+Bmat(k,i)*GBmat(k,j)
107 write (iout,*) "Matrix Cmat"
108 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
110 call matinvert(nbond,MAXRES2,Cmat,Cinv)
112 write (iout,*) "Matrix Cinv"
113 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cinv)
119 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
124 write (iout,*) "Matrix Tmat"
125 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Tmat)
135 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
140 write (iout,*) "Matrix Pmat"
141 call MATOUT(dimen,dimen,MAXRES6,MAXRES6,Pmat)
148 Td(i)=Td(i)+vbl*Tmat(i,ind)
151 if (itype(k).ne.10) then
153 Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
158 write (iout,*) "Vector Td"
160 write (iout,'(i5,f10.5)') i,Td(i)
163 call stochastic_force(stochforcvec)
165 write (iout,*) "stochforcvec"
167 write (iout,*) i,stochforcvec(i)
171 zapas(j)=-gcart(j,0)+stochforcvec(j)
173 dC_work(j)=dC_old(j,0)
179 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
180 dC_work(ind)=dC_old(j,i)
184 if (itype(i).ne.10) then
187 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
188 dC_work(ind)=dC_old(j,i+nres)
194 write (iout,*) "Initial d_t_work"
196 write (iout,*) i,d_t_work(i)
203 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
210 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
214 write (iout,*) "Final d_t_work and zapas"
216 write (iout,*) i,d_t_work(i),zapas(i)
230 dc_work(ind+j)=dc(j,i)
236 d_t(j,i+nres)=d_t_work(ind+j)
237 dc(j,i+nres)=zapas(ind+j)
238 dc_work(ind+j)=dc(j,i+nres)
244 write (iout,*) "Before correction for rotational lengthening"
245 write (iout,*) "New coordinates",
246 & " and differences between actual and standard bond lengths"
251 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
252 & i,(dC(j,i),j=1,3),xx
255 if (itype(i).ne.10) then
257 xx=vbld(i+nres)-vbldsc0(1,itype(i))
258 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
259 & i,(dC(j,i+nres),j=1,3),xx
263 c Second correction (rotational lengthening)
269 blen2 = scalar(dc(1,i),dc(1,i))
270 ppvec(ind)=2*vbl**2-blen2
271 diffbond=dabs(vbl-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 (itype(i).ne.10) then
285 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
286 ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
287 diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2))
288 if (diffbond.gt.diffmax) diffmax=diffbond
289 if (ppvec(ind).gt.0.0d0) then
290 ppvec(ind)=dsqrt(ppvec(ind))
295 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
299 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
300 if (diffmax.lt.difftol) goto 10
304 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
310 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
321 dc_work(ind+j)=zapas(ind+j)
326 if (itype(i).ne.10) then
328 dc(j,i+nres)=zapas(ind+j)
329 dc_work(ind+j)=zapas(ind+j)
334 c Building the chain from the newly calculated coordinates
337 if (large.and. mod(itime,ntwe).eq.0) then
338 write (iout,*) "Cartesian and internal coordinates: step 1"
341 write (iout,'(a)') "Potential forces"
343 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),
344 & (-gxcart(j,i),j=1,3)
346 write (iout,'(a)') "Stochastic forces"
348 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),
349 & (stochforc(j,i+nres),j=1,3)
351 write (iout,'(a)') "Velocities"
353 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
354 & (d_t(j,i+nres),j=1,3)
359 write (iout,*) "After correction for rotational lengthening"
360 write (iout,*) "New coordinates",
361 & " and differences between actual and standard bond lengths"
366 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
367 & i,(dC(j,i),j=1,3),xx
370 if (itype(i).ne.10) then
372 xx=vbld(i+nres)-vbldsc0(1,itype(i))
373 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
374 & i,(dC(j,i+nres),j=1,3),xx
379 c write (iout,*) "Too many attempts at correcting the bonds"
387 c Calculate energy and forces
389 call etotal(potEcomp)
390 potE=potEcomp(0)-potEcomp(27)
394 c Calculate the kinetic and total energy and the kinetic temperature
397 t_enegrad=t_enegrad+MPI_Wtime()-tt0
399 t_enegrad=t_enegrad+tcpu()-tt0
402 kinetic_T=2.0d0/(dimen*Rb)*EK