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