wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-M / stochfric.F
1       subroutine friction_force
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'COMMON.VAR'
5       include 'COMMON.CHAIN'
6       include 'COMMON.DERIV'
7       include 'COMMON.GEO'
8       include 'COMMON.LOCAL'
9       include 'COMMON.INTERACT'
10       include 'COMMON.MD'
11 #ifndef LANG0
12       include 'COMMON.LANGEVIN'
13 #else
14       include 'COMMON.LANGEVIN.lang0'
15 #endif
16       include 'COMMON.IOUNITS'
17       double precision gamvec(MAXRES6)
18       common /syfek/ gamvec
19       double precision vv(3),vvtot(3,maxres),v_work(MAXRES6),
20      & ginvfric(maxres2,maxres2)
21       common /przechowalnia/ ginvfric
22       
23       logical lprn /.false./, checkmode /.false./
24
25       do i=0,MAXRES2
26         do j=1,3
27           friction(j,i)=0.0d0
28         enddo
29       enddo
30   
31       do j=1,3
32         d_t_work(j)=d_t(j,0)
33       enddo
34       ind=3
35       do i=nnt,nct-1
36         do j=1,3
37           d_t_work(ind+j)=d_t(j,i)
38         enddo
39         ind=ind+3
40       enddo
41       do i=nnt,nct
42         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
43           do j=1,3
44             d_t_work(ind+j)=d_t(j,i+nres)
45           enddo
46           ind=ind+3
47         endif
48       enddo
49
50       call fricmat_mult(d_t_work,fric_work)
51       
52       if (.not.checkmode) return
53
54       if (lprn) then
55         write (iout,*) "d_t_work and fric_work"
56         do i=1,3*dimen
57           write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
58         enddo
59       endif
60       do j=1,3
61         friction(j,0)=fric_work(j)
62       enddo
63       ind=3
64       do i=nnt,nct-1
65         do j=1,3
66           friction(j,i)=fric_work(ind+j)
67         enddo
68         ind=ind+3
69       enddo
70       do i=nnt,nct
71         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
72           do j=1,3
73             friction(j,i+nres)=fric_work(ind+j)
74           enddo
75           ind=ind+3
76         endif
77       enddo
78       if (lprn) then
79         write(iout,*) "Friction backbone"
80         do i=0,nct-1
81           write(iout,'(i5,3e15.5,5x,3e15.5)') 
82      &     i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
83         enddo
84         write(iout,*) "Friction side chain"
85         do i=nnt,nct
86           write(iout,'(i5,3e15.5,5x,3e15.5)') 
87      &     i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
88         enddo   
89       endif
90       if (lprn) then
91         do j=1,3
92           vv(j)=d_t(j,0)
93         enddo
94         do i=nnt,nct
95           do j=1,3
96             vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
97             vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
98             vv(j)=vv(j)+d_t(j,i)
99           enddo
100         enddo
101         write (iout,*) "vvtot backbone and sidechain"
102         do i=nnt,nct
103           write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),
104      &     (vvtot(j,i+nres),j=1,3)
105         enddo
106         ind=0
107         do i=nnt,nct-1
108           do j=1,3
109             v_work(ind+j)=vvtot(j,i)
110           enddo
111           ind=ind+3
112         enddo
113         do i=nnt,nct
114           do j=1,3
115             v_work(ind+j)=vvtot(j,i+nres)
116           enddo
117           ind=ind+3
118         enddo
119         write (iout,*) "v_work gamvec and site-based friction forces"
120         do i=1,dimen1
121           write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),
122      &      gamvec(i)*v_work(i) 
123         enddo
124 c        do i=1,dimen
125 c          fric_work1(i)=0.0d0
126 c          do j=1,dimen1
127 c            fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
128 c          enddo
129 c        enddo  
130 c        write (iout,*) "fric_work and fric_work1"
131 c        do i=1,dimen
132 c          write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
133 c        enddo 
134         do i=1,dimen
135           do j=1,dimen
136             ginvfric(i,j)=0.0d0
137             do k=1,dimen
138               ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
139             enddo
140           enddo
141         enddo
142         write (iout,*) "ginvfric"
143         do i=1,dimen
144           write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
145         enddo
146         write (iout,*) "symmetry check"
147         do i=1,dimen
148           do j=1,i-1
149             write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
150           enddo   
151         enddo
152       endif 
153       return
154       end
155 c-----------------------------------------------------
156       subroutine stochastic_force(stochforcvec)
157       implicit real*8 (a-h,o-z)
158       include 'DIMENSIONS'
159 #ifdef MPI
160       include 'mpif.h'
161 #endif
162       include 'COMMON.VAR'
163       include 'COMMON.CHAIN'
164       include 'COMMON.DERIV'
165       include 'COMMON.GEO'
166       include 'COMMON.LOCAL'
167       include 'COMMON.INTERACT'
168       include 'COMMON.MD'
169       include 'COMMON.TIME1'
170 #ifndef LANG0
171       include 'COMMON.LANGEVIN'
172 #else
173       include 'COMMON.LANGEVIN.lang0'
174 #endif
175       include 'COMMON.IOUNITS'
176       
177       double precision x,sig,lowb,highb,
178      & ff(3),force(3,0:MAXRES2),zeta2,lowb2,
179      & highb2,sig2,forcvec(MAXRES6),stochforcvec(MAXRES6)
180       logical lprn /.false./
181       do i=0,MAXRES2
182         do j=1,3
183           stochforc(j,i)=0.0d0
184         enddo
185       enddo
186       x=0.0d0   
187
188 #ifdef MPI
189       time00=MPI_Wtime()
190 #else
191       time00=tcpu()
192 #endif
193 c Compute the stochastic forces acting on bodies. Store in force.
194       do i=nnt,nct-1
195         sig=stdforcp(i)
196         lowb=-5*sig
197         highb=5*sig
198         do j=1,3
199           force(j,i)=anorm_distr(x,sig,lowb,highb)
200         enddo
201       enddo
202       do i=nnt,nct
203         sig2=stdforcsc(i)
204         lowb2=-5*sig2
205         highb2=5*sig2
206         do j=1,3
207           force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
208         enddo
209       enddo
210 #ifdef MPI
211       time_fsample=time_fsample+MPI_Wtime()-time00
212 #else
213       time_fsample=time_fsample+tcpu()-time00
214 #endif
215 c Compute the stochastic forces acting on virtual-bond vectors.
216       do j=1,3
217         ff(j)=0.0d0
218       enddo
219       do i=nct-1,nnt,-1
220         do j=1,3
221           stochforc(j,i)=ff(j)+0.5d0*force(j,i)
222         enddo
223         do j=1,3
224           ff(j)=ff(j)+force(j,i)
225         enddo
226         if (itype(i+1).ne.ntyp1) then
227           do j=1,3
228             stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
229             ff(j)=ff(j)+force(j,i+nres+1)
230           enddo
231         endif
232       enddo 
233       do j=1,3
234         stochforc(j,0)=ff(j)+force(j,nnt+nres)
235       enddo
236       do i=nnt,nct
237         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
238           do j=1,3
239             stochforc(j,i+nres)=force(j,i+nres)
240           enddo
241         endif
242       enddo 
243
244       do j=1,3
245         stochforcvec(j)=stochforc(j,0)
246       enddo
247       ind=3
248       do i=nnt,nct-1
249         do j=1,3 
250           stochforcvec(ind+j)=stochforc(j,i)
251         enddo
252         ind=ind+3
253       enddo
254       do i=nnt,nct
255         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
256           do j=1,3
257             stochforcvec(ind+j)=stochforc(j,i+nres)
258           enddo
259           ind=ind+3
260         endif
261       enddo
262       if (lprn) then
263         write (iout,*) "stochforcvec"
264         do i=1,3*dimen
265           write(iout,'(i5,e15.5)') i,stochforcvec(i)
266         enddo
267         write(iout,*) "Stochastic forces backbone"
268         do i=0,nct-1
269           write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
270         enddo
271         write(iout,*) "Stochastic forces side chain"
272         do i=nnt,nct
273           write(iout,'(i5,3e15.5)') 
274      &      i,(stochforc(j,i+nres),j=1,3)
275         enddo   
276       endif
277
278       if (lprn) then
279
280       ind=0
281       do i=nnt,nct-1
282         write (iout,*) i,ind
283         do j=1,3
284           forcvec(ind+j)=force(j,i)
285         enddo
286         ind=ind+3
287       enddo
288       do i=nnt,nct
289         write (iout,*) i,ind
290         do j=1,3
291           forcvec(j+ind)=force(j,i+nres)
292         enddo
293         ind=ind+3
294       enddo 
295
296       write (iout,*) "forcvec"
297       ind=0
298       do i=nnt,nct-1
299         do j=1,3
300           write (iout,'(2i3,2f10.5)') i,j,force(j,i),
301      &      forcvec(ind+j)
302         enddo
303         ind=ind+3
304       enddo
305       do i=nnt,nct
306         do j=1,3
307           write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),
308      &     forcvec(ind+j)
309         enddo
310         ind=ind+3
311       enddo
312
313       endif
314
315       return
316       end
317 c------------------------------------------------------------------
318       subroutine setup_fricmat
319       implicit real*8 (a-h,o-z)
320 #ifdef MPI
321       include 'mpif.h'
322 #endif
323       include 'DIMENSIONS'
324       include 'COMMON.VAR'
325       include 'COMMON.CHAIN'
326       include 'COMMON.DERIV'
327       include 'COMMON.GEO'
328       include 'COMMON.LOCAL'
329       include 'COMMON.INTERACT'
330       include 'COMMON.MD'
331       include 'COMMON.SETUP'
332       include 'COMMON.TIME1'
333 c      integer licznik /0/
334 c      save licznik
335 #ifndef LANG0
336       include 'COMMON.LANGEVIN'
337 #else
338       include 'COMMON.LANGEVIN.lang0'
339 #endif
340       include 'COMMON.IOUNITS'
341       integer IERROR
342       integer i,j,ind,ind1,m
343       logical lprn /.false./
344       double precision dtdi,gamvec(MAXRES2),
345      &  ginvfric(maxres2,maxres2),Ghalf(mmaxres2),fcopy(maxres2,maxres2)
346       common /syfek/ gamvec
347       double precision work(8*maxres2)
348       integer iwork(maxres2)
349       common /przechowalnia/ ginvfric,Ghalf,fcopy
350 #ifdef MPI
351       if (fg_rank.ne.king) goto 10
352 #endif
353 c  Zeroing out fricmat
354       do i=1,dimen
355         do j=1,dimen
356           fricmat(i,j)=0.0d0
357         enddo    
358       enddo
359 c  Load the friction coefficients corresponding to peptide groups
360       ind1=0
361       do i=nnt,nct-1
362         ind1=ind1+1
363         gamvec(ind1)=gamp
364       enddo
365 c  Load the friction coefficients corresponding to side chains
366       m=nct-nnt
367       ind=0
368 C      gamsc(ntyp1)=1.0d0
369       do i=nnt,nct
370         ind=ind+1
371         ii = ind+m
372         iti=itype(i)
373         gamvec(ii)=gamsc(iabs(iti))
374       enddo
375       if (surfarea) call sdarea(gamvec)
376 c      if (lprn) then
377 c        write (iout,*) "Matrix A and vector gamma"
378 c        do i=1,dimen1
379 c          write (iout,'(i2,$)') i
380 c          do j=1,dimen
381 c            write (iout,'(f4.1,$)') A(i,j)
382 c          enddo
383 c          write (iout,'(f8.3)') gamvec(i)
384 c        enddo
385 c      endif
386       if (lprn) then
387         write (iout,*) "Vector gamvec"
388         do i=1,dimen1
389           write (iout,'(i5,f10.5)') i, gamvec(i)
390         enddo
391       endif
392         
393 c The friction matrix       
394       do k=1,dimen
395        do i=1,dimen
396          dtdi=0.0d0
397          do j=1,dimen1
398            dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
399          enddo
400          fricmat(k,i)=dtdi
401        enddo
402       enddo 
403
404       if (lprn) then
405         write (iout,'(//a)') "Matrix fricmat"
406         call matout2(dimen,dimen,maxres2,maxres2,fricmat)
407       endif 
408       if (lang.eq.2 .or. lang.eq.3) then
409 c Mass-scale the friction matrix if non-direct integration will be performed
410       do i=1,dimen
411         do j=1,dimen
412           Ginvfric(i,j)=0.0d0
413           do k=1,dimen
414             do l=1,dimen
415               Ginvfric(i,j)=Ginvfric(i,j)+
416      &          Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
417             enddo
418           enddo
419         enddo
420       enddo
421 c Diagonalize the friction matrix
422       ind=0
423       do i=1,dimen
424         do j=1,i
425           ind=ind+1
426           Ghalf(ind)=Ginvfric(i,j)
427         enddo
428       enddo
429       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
430      &  ierr,iwork)
431       if (lprn) then
432         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
433      &    " mass-scaled friction matrix"
434         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
435       endif
436 c Precompute matrices for tinker stochastic integrator
437 #ifndef LANG0
438       do i=1,dimen
439         do j=1,dimen
440           mt1(i,j)=0.0d0
441           mt2(i,j)=0.0d0
442           do k=1,dimen
443             mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
444             mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)              
445           enddo
446           mt3(j,i)=mt1(i,j)
447         enddo
448       enddo
449 #endif
450       else if (lang.eq.4) then
451 c Diagonalize the friction matrix
452       ind=0
453       do i=1,dimen
454         do j=1,i
455           ind=ind+1
456           Ghalf(ind)=fricmat(i,j)
457         enddo
458       enddo
459       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
460      &  ierr,iwork)
461       if (lprn) then
462         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
463      &    " friction matrix"
464         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
465       endif
466 c Determine the number of zero eigenvalues of the friction matrix
467       nzero=max0(dimen-dimen1,0)
468 c      do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
469 c        nzero=nzero+1
470 c      enddo
471       write (iout,*) "Number of zero eigenvalues:",nzero
472       do i=1,dimen
473         do j=1,dimen
474           fricmat(i,j)=0.0d0
475           do k=nzero+1,dimen
476             fricmat(i,j)=fricmat(i,j)
477      &        +fricvec(i,k)*fricvec(j,k)/fricgam(k) 
478           enddo
479         enddo
480       enddo
481       if (lprn) then
482         write (iout,'(//a)') "Generalized inverse of fricmat"
483         call matout(dimen,dimen,maxres6,maxres6,fricmat)
484       endif 
485       endif
486 #ifdef MPI
487   10  continue
488       if (nfgtasks.gt.1) then
489         if (fg_rank.eq.0) then
490 c The matching BROADCAST for fg processors is called in ERGASTULUM
491 #ifdef MPI
492           time00=MPI_Wtime()
493 #else
494           time00=tcpu()
495 #endif
496           call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
497 #ifdef MPI
498           time_Bcast=time_Bcast+MPI_Wtime()-time00
499 #else
500           time_Bcast=time_Bcast+tcpu()-time00
501 #endif
502 c          print *,"Processor",myrank,
503 c     &       " BROADCAST iorder in SETUP_FRICMAT"
504         endif
505 c       licznik=licznik+1
506 c        write (iout,*) "setup_fricmat licznik",licznik
507 #ifdef MPI
508         time00=MPI_Wtime()
509 #else
510         time00=tcpu()
511 #endif
512 c Scatter the friction matrix
513         call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
514      &    nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
515      &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
516 #ifdef TIMING
517 #ifdef MPI
518         time_scatter=time_scatter+MPI_Wtime()-time00
519         time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
520 #else
521         time_scatter=time_scatter+tcpu()-time00
522         time_scatter_fmat=time_scatter_fmat+tcpu()-time00
523 #endif
524 #endif
525         do i=1,dimen
526           do j=1,2*my_ng_count
527             fricmat(j,i)=fcopy(i,j)
528           enddo
529         enddo
530 c        write (iout,*) "My chunk of fricmat"
531 c        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
532       endif
533 #endif
534       return
535       end
536 c-------------------------------------------------------------------------------
537       subroutine sdarea(gamvec)
538 c
539 c Scale the friction coefficients according to solvent accessible surface areas
540 c Code adapted from TINKER
541 c AL 9/3/04
542 c
543       implicit real*8 (a-h,o-z)
544       include 'DIMENSIONS'
545       include 'COMMON.CONTROL'
546       include 'COMMON.VAR'
547       include 'COMMON.MD'
548 #ifndef LANG0
549       include 'COMMON.LANGEVIN'
550 #else
551       include 'COMMON.LANGEVIN.lang0'
552 #endif
553       include 'COMMON.CHAIN'
554       include 'COMMON.DERIV'
555       include 'COMMON.GEO'
556       include 'COMMON.LOCAL'
557       include 'COMMON.INTERACT'
558       include 'COMMON.IOUNITS'
559       include 'COMMON.NAMES'
560       double precision radius(maxres2),gamvec(maxres2)
561       parameter (twosix=1.122462048309372981d0)
562       logical lprn /.false./
563 c
564 c     determine new friction coefficients every few SD steps
565 c
566 c     set the atomic radii to estimates of sigma values
567 c
568 c      print *,"Entered sdarea"
569       probe = 0.0d0
570       
571       do i=1,2*nres
572         radius(i)=0.0d0
573       enddo
574 c  Load peptide group radii
575       do i=nnt,nct-1
576         radius(i)=pstok
577       enddo
578 c  Load side chain radii
579       do i=nnt,nct
580         iti=itype(i)
581         radius(i+nres)=restok(iti)
582       enddo
583 c      do i=1,2*nres
584 c        write (iout,*) "i",i," radius",radius(i) 
585 c      enddo
586       do i = 1, 2*nres
587          radius(i) = radius(i) / twosix
588          if (radius(i) .ne. 0.0d0)  radius(i) = radius(i) + probe
589       end do
590 c
591 c     scale atomic friction coefficients by accessible area
592 c
593       if (lprn) write (iout,*) 
594      &  "Original gammas, surface areas, scaling factors, new gammas, ",
595      &  "std's of stochastic forces"
596       ind=0
597       do i=nnt,nct-1
598         if (radius(i).gt.0.0d0) then
599           call surfatom (i,area,radius)
600           ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
601           if (lprn) write (iout,'(i5,3f10.5,$)') 
602      &      i,gamvec(ind+1),area,ratio
603           do j=1,3
604             ind=ind+1
605             gamvec(ind) = ratio * gamvec(ind)
606           enddo
607           stdforcp(i)=stdfp*dsqrt(gamvec(ind))
608           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
609         endif
610       enddo
611       do i=nnt,nct
612         if (radius(i+nres).gt.0.0d0) then
613           call surfatom (i+nres,area,radius)
614           ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
615           if (lprn) write (iout,'(i5,3f10.5,$)') 
616      &      i,gamvec(ind+1),area,ratio
617           do j=1,3
618             ind=ind+1 
619             gamvec(ind) = ratio * gamvec(ind)
620           enddo
621           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
622           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
623         endif
624       enddo
625
626       return
627       end