050aac869c4290b73bfd7a1980e7c397cc119c80
[unres.git] / source / unres / src_MD / md-diff / md / brown_step.f
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)
7       include 'DIMENSIONS'
8       include 'COMMON.CONTROL'
9       include 'COMMON.VAR'
10       include 'COMMON.MD'
11       include 'COMMON.LANGEVIN'
12       include 'COMMON.CHAIN'
13       include 'COMMON.DERIV'
14       include 'COMMON.GEO'
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)
21       integer ilen,rstcount
22       external ilen
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),
27      & ppvec(maxres2)
28       common /stochcalc/ stochforcvec
29       common /gucio/ cm, energia
30       integer itime
31       logical lprn /.false./,lprn1 /.false./
32       integer maxiter /5/
33       double precision difftol /1.0d-5/
34       nbond=nct-nnt
35       do i=nnt,nct
36         if (itype(i).ne.10) nbond=nbond+1
37       enddo
38 c
39       if (lprn1) then
40         write (iout,*) "Generalized inverse of fricmat"
41         call matout(dimen,dimen,MAXRES6,MAXRES6,fricmat)
42       endif 
43       do i=1,dimen
44         do j=1,nbond
45           Bmat(i,j)=0.0d0
46         enddo
47       enddo
48       ind=3
49       ind1=0
50       do i=nnt,nct-1
51         ind1=ind1+1
52         do j=1,3
53           Bmat(ind+j,ind1)=dC_norm(j,i)
54         enddo
55         ind=ind+3
56       enddo
57       do i=nnt,nct
58         if (itype(i).ne.10) then
59           ind1=ind1+1
60           do j=1,3
61             Bmat(ind+j,ind1)=dC_norm(j,i+nres)
62           enddo
63           ind=ind+3
64         endif
65       enddo
66       if (lprn1) then 
67         write (iout,*) "Matrix Bmat"
68         call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Bmat)
69       endif
70       do i=1,dimen
71         do j=1,nbond
72           GBmat(i,j)=0.0d0
73           do k=1,dimen
74             GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
75           enddo
76         enddo
77       enddo   
78       if (lprn1) then
79         write (iout,*) "Matrix GBmat"
80         call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Gbmat)
81       endif
82       do i=1,nbond
83         do j=1,nbond
84           Cmat(i,j)=0.0d0
85           do k=1,dimen
86             Cmat(i,j)=Cmat(i,j)+Bmat(k,i)*GBmat(k,j)
87           enddo
88         enddo
89       enddo
90       if (lprn1) then
91         write (iout,*) "Matrix Cmat"
92         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
93       endif
94       call matinvert(nbond,MAXRES2,Cmat,Cinv) 
95       if (lprn1) then
96         write (iout,*) "Matrix Cinv"
97         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cinv)
98       endif
99       do i=1,dimen
100         do j=1,nbond
101           Tmat(i,j)=0.0d0
102           do k=1,nbond
103             Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
104           enddo
105         enddo
106       enddo
107       if (lprn1) then
108         write (iout,*) "Matrix Tmat"
109         call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Tmat)
110       endif
111       do i=1,dimen
112         do j=1,dimen
113           if (i.eq.j) then
114             Pmat(i,j)=1.0d0
115           else
116             Pmat(i,j)=0.0d0
117           endif
118           do k=1,nbond
119             Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
120           enddo
121         enddo
122       enddo
123       if (lprn1) then
124         write (iout,*) "Matrix Pmat"
125         call MATOUT(dimen,dimen,MAXRES6,MAXRES6,Pmat)
126       endif
127       do i=1,dimen
128         Td(i)=0.0d0
129         ind=0
130         do k=nnt,nct-1
131           ind=ind+1
132           Td(i)=Td(i)+vbl*Tmat(i,ind)
133         enddo
134         do k=nnt,nct
135           if (itype(k).ne.10) then
136             ind=ind+1
137             Td(i)=Td(i)+vbldsc0(itype(k))*Tmat(i,ind)
138           endif
139         enddo
140       enddo 
141       if (lprn1) then
142         write (iout,*) "Vector Td"
143         do i=1,dimen
144           write (iout,'(i5,f10.5)') i,Td(i)
145         enddo
146       endif
147       call stochastic_force(stochforcvec)
148       if (lprn) then
149         write (iout,*) "stochforcvec"
150         do i=1,dimen
151           write (iout,*) i,stochforcvec(i)
152         enddo
153       endif
154       do j=1,3
155         zapas(j)=-gcart(j,0)+stochforcvec(j)
156         d_t_work(j)=d_t(j,0)
157         dC_work(j)=dC_old(j,0)
158       enddo
159       ind=3      
160       do i=nnt,nct-1
161         do j=1,3
162           ind=ind+1
163           zapas(ind)=-gcart(j,i)+stochforcvec(ind)
164           dC_work(ind)=dC_old(j,i)
165         enddo
166       enddo
167       do i=nnt,nct
168         if (itype(i).ne.10) then
169           do j=1,3
170             ind=ind+1
171             zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
172             dC_work(ind)=dC_old(j,i+nres)
173           enddo
174         endif
175       enddo
176
177       if (lprn) then
178         write (iout,*) "Initial d_t_work"
179         do i=1,dimen
180           write (iout,*) i,d_t_work(i)
181         enddo
182       endif
183
184       do i=1,dimen
185         d_t_work(i)=0.0d0
186         do j=1,dimen
187           d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
188         enddo
189       enddo
190
191       do i=1,dimen
192         zapas(i)=Td(i)
193         do j=1,dimen
194           zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
195         enddo
196       enddo
197       if (lprn1) then
198         write (iout,*) "Final d_t_work and zapas"
199         do i=1,dimen
200           write (iout,*) i,d_t_work(i),zapas(i)
201         enddo
202       endif
203
204       do j=1,3
205         d_t(j,0)=d_t_work(j)
206         dc(j,0)=zapas(j)
207         dc_work(j)=dc(j,0)
208       enddo
209       ind=3
210       do i=nnt,nct-1
211         do j=1,3
212           d_t(j,i)=d_t_work(i)
213           dc(j,i)=zapas(ind+j)
214           dc_work(ind+j)=dc(j,i)
215         enddo
216         ind=ind+3
217       enddo
218       do i=nnt,nct
219         do j=1,3
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)
223         enddo
224         ind=ind+3
225       enddo
226       if (lprn) then
227         call chainbuild_cart
228         write (iout,*) "Before correction for rotational lengthening"
229         write (iout,*) "New coordinates",
230      &  " and differences between actual and standard bond lengths"
231         ind=0
232         do i=nnt,nct-1
233           ind=ind+1
234           xx=vbld(i+1)-vbl
235           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
236      &        i,(dC(j,i),j=1,3),xx
237         enddo
238         do i=nnt,nct
239           if (itype(i).ne.10) then
240             ind=ind+1
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
244           endif
245         enddo
246       endif
247 c Second correction (rotational lengthening)
248 c      do iter=1,maxiter
249       diffmax=0.0d0
250       ind=0
251       do i=nnt,nct-1
252         ind=ind+1
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))
259         else
260           ppvec(ind)=0.0d0
261         endif
262         if (lprn) then
263           write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
264         endif
265       enddo
266       do i=nnt,nct
267         if (itype(i).ne.10) then
268           ind=ind+1
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))
275           else
276             ppvec(ind)=0.0d0
277           endif
278           if (lprn) then
279             write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
280           endif
281         endif
282       enddo
283       if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
284       if (diffmax.lt.difftol) goto 10
285       do i=1,dimen
286         Td(i)=0.0d0
287         do j=1,nbond
288           Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
289         enddo
290       enddo 
291       do i=1,dimen
292         zapas(i)=Td(i)
293         do j=1,dimen
294           zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
295         enddo
296       enddo
297       do j=1,3
298         dc(j,0)=zapas(j)
299         dc_work(j)=zapas(j)
300       enddo
301       ind=3
302       do i=nnt,nct-1
303         do j=1,3
304           dc(j,i)=zapas(ind+j)
305           dc_work(ind+j)=zapas(ind+j)
306         enddo
307         ind=ind+3
308       enddo
309       do i=nnt,nct
310         if (itype(i).ne.10) then
311           do j=1,3
312             dc(j,i+nres)=zapas(ind+j)
313             dc_work(ind+j)=zapas(ind+j)
314           enddo
315           ind=ind+3
316         endif
317       enddo 
318 c   Building the chain from the newly calculated coordinates    
319       call chainbuild_cart
320       if (large.and. mod(itime,ntwe).eq.0) then
321         write (iout,*) "Cartesian and internal coordinates: step 1"
322         call cartprint
323         call intout
324         write (iout,'(a)') "Potential forces"
325         do i=0,nres
326           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),
327      &    (-gxcart(j,i),j=1,3)
328         enddo
329         write (iout,'(a)') "Stochastic forces"
330         do i=0,nres
331           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),
332      &    (stochforc(j,i+nres),j=1,3)
333         enddo
334         write (iout,'(a)') "Velocities"
335         do i=0,nres
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)
338         enddo
339       endif
340       if (lprn) then
341         write (iout,*) "After correction for rotational lengthening"
342         write (iout,*) "New coordinates",
343      &  " and differences between actual and standard bond lengths"
344         ind=0
345         do i=nnt,nct-1
346           ind=ind+1
347           xx=vbld(i+1)-vbl
348           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
349      &        i,(dC(j,i),j=1,3),xx
350         enddo
351         do i=nnt,nct
352           if (itype(i).ne.10) then
353             ind=ind+1
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
357           endif
358         enddo
359       endif
360 c      ENDDO
361 c      write (iout,*) "Too many attempts at correcting the bonds"
362 c      stop
363    10 continue
364       tt0 = tcpu()
365 c Calculate energy and forces
366       call zerograd
367       call etotal(energia)
368       potE=energia(0)-energia(20)
369       call cartgrad
370       totT=totT+d_time
371 c  Calculate the kinetic and total energy and the kinetic temperature
372       call kinetic(EK)
373       t_enegrad=t_enegrad+tcpu()-tt0
374       totE=EK+potE
375       kinetic_T=2.0d0/(dimen*Rb)*EK
376       return
377       end