Added src_Eshel (decoy processing for threading)
[unres.git] / source / unres / src_Eshel / SRC-SURPLUS / 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) 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) 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.21) 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) 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) 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       do i=nnt,nct
369         ind=ind+1
370         ii = ind+m
371         iti=itype(i)
372         gamvec(ii)=gamsc(iti)
373       enddo
374       if (surfarea) call sdarea(gamvec)
375 c      if (lprn) then
376 c        write (iout,*) "Matrix A and vector gamma"
377 c        do i=1,dimen1
378 c          write (iout,'(i2,$)') i
379 c          do j=1,dimen
380 c            write (iout,'(f4.1,$)') A(i,j)
381 c          enddo
382 c          write (iout,'(f8.3)') gamvec(i)
383 c        enddo
384 c      endif
385       if (lprn) then
386         write (iout,*) "Vector gamvec"
387         do i=1,dimen1
388           write (iout,'(i5,f10.5)') i, gamvec(i)
389         enddo
390       endif
391         
392 c The friction matrix       
393       do k=1,dimen
394        do i=1,dimen
395          dtdi=0.0d0
396          do j=1,dimen1
397            dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
398          enddo
399          fricmat(k,i)=dtdi
400        enddo
401       enddo 
402
403       if (lprn) then
404         write (iout,'(//a)') "Matrix fricmat"
405         call matout2(dimen,dimen,maxres2,maxres2,fricmat)
406       endif 
407       if (lang.eq.2 .or. lang.eq.3) then
408 c Mass-scale the friction matrix if non-direct integration will be performed
409       do i=1,dimen
410         do j=1,dimen
411           Ginvfric(i,j)=0.0d0
412           do k=1,dimen
413             do l=1,dimen
414               Ginvfric(i,j)=Ginvfric(i,j)+
415      &          Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
416             enddo
417           enddo
418         enddo
419       enddo
420 c Diagonalize the friction matrix
421       ind=0
422       do i=1,dimen
423         do j=1,i
424           ind=ind+1
425           Ghalf(ind)=Ginvfric(i,j)
426         enddo
427       enddo
428       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
429      &  ierr,iwork)
430       if (lprn) then
431         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
432      &    " mass-scaled friction matrix"
433         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
434       endif
435 c Precompute matrices for tinker stochastic integrator
436 #ifndef LANG0
437       do i=1,dimen
438         do j=1,dimen
439           mt1(i,j)=0.0d0
440           mt2(i,j)=0.0d0
441           do k=1,dimen
442             mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
443             mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)              
444           enddo
445           mt3(j,i)=mt1(i,j)
446         enddo
447       enddo
448 #endif
449       else if (lang.eq.4) then
450 c Diagonalize the friction matrix
451       ind=0
452       do i=1,dimen
453         do j=1,i
454           ind=ind+1
455           Ghalf(ind)=fricmat(i,j)
456         enddo
457       enddo
458       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
459      &  ierr,iwork)
460       if (lprn) then
461         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
462      &    " friction matrix"
463         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
464       endif
465 c Determine the number of zero eigenvalues of the friction matrix
466       nzero=max0(dimen-dimen1,0)
467 c      do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
468 c        nzero=nzero+1
469 c      enddo
470       write (iout,*) "Number of zero eigenvalues:",nzero
471       do i=1,dimen
472         do j=1,dimen
473           fricmat(i,j)=0.0d0
474           do k=nzero+1,dimen
475             fricmat(i,j)=fricmat(i,j)
476      &        +fricvec(i,k)*fricvec(j,k)/fricgam(k) 
477           enddo
478         enddo
479       enddo
480       if (lprn) then
481         write (iout,'(//a)') "Generalized inverse of fricmat"
482         call matout(dimen,dimen,maxres6,maxres6,fricmat)
483       endif 
484       endif
485 #ifdef MPI
486   10  continue
487       if (nfgtasks.gt.1) then
488         if (fg_rank.eq.0) then
489 c The matching BROADCAST for fg processors is called in ERGASTULUM
490 #ifdef MPI
491           time00=MPI_Wtime()
492 #else
493           time00=tcpu()
494 #endif
495           call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
496 #ifdef MPI
497           time_Bcast=time_Bcast+MPI_Wtime()-time00
498 #else
499           time_Bcast=time_Bcast+tcpu()-time00
500 #endif
501 c          print *,"Processor",myrank,
502 c     &       " BROADCAST iorder in SETUP_FRICMAT"
503         endif
504 c       licznik=licznik+1
505 c        write (iout,*) "setup_fricmat licznik",licznik
506 #ifdef MPI
507         time00=MPI_Wtime()
508 #else
509         time00=tcpu()
510 #endif
511 c Scatter the friction matrix
512         call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
513      &    nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
514      &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
515 #ifdef TIMING
516 #ifdef MPI
517         time_scatter=time_scatter+MPI_Wtime()-time00
518         time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
519 #else
520         time_scatter=time_scatter+tcpu()-time00
521         time_scatter_fmat=time_scatter_fmat+tcpu()-time00
522 #endif
523 #endif
524         do i=1,dimen
525           do j=1,2*my_ng_count
526             fricmat(j,i)=fcopy(i,j)
527           enddo
528         enddo
529 c        write (iout,*) "My chunk of fricmat"
530 c        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
531       endif
532 #endif
533       return
534       end
535 c-------------------------------------------------------------------------------
536       subroutine sdarea(gamvec)
537 c
538 c Scale the friction coefficients according to solvent accessible surface areas
539 c Code adapted from TINKER
540 c AL 9/3/04
541 c
542       implicit real*8 (a-h,o-z)
543       include 'DIMENSIONS'
544       include 'COMMON.CONTROL'
545       include 'COMMON.VAR'
546       include 'COMMON.MD'
547 #ifndef LANG0
548       include 'COMMON.LANGEVIN'
549 #else
550       include 'COMMON.LANGEVIN.lang0'
551 #endif
552       include 'COMMON.CHAIN'
553       include 'COMMON.DERIV'
554       include 'COMMON.GEO'
555       include 'COMMON.LOCAL'
556       include 'COMMON.INTERACT'
557       include 'COMMON.IOUNITS'
558       include 'COMMON.NAMES'
559       double precision radius(maxres2),gamvec(maxres2)
560       parameter (twosix=1.122462048309372981d0)
561       logical lprn /.false./
562 c
563 c     determine new friction coefficients every few SD steps
564 c
565 c     set the atomic radii to estimates of sigma values
566 c
567 c      print *,"Entered sdarea"
568       probe = 0.0d0
569       
570       do i=1,2*nres
571         radius(i)=0.0d0
572       enddo
573 c  Load peptide group radii
574       do i=nnt,nct-1
575         radius(i)=pstok
576       enddo
577 c  Load side chain radii
578       do i=nnt,nct
579         iti=itype(i)
580         radius(i+nres)=restok(iti)
581       enddo
582 c      do i=1,2*nres
583 c        write (iout,*) "i",i," radius",radius(i) 
584 c      enddo
585       do i = 1, 2*nres
586          radius(i) = radius(i) / twosix
587          if (radius(i) .ne. 0.0d0)  radius(i) = radius(i) + probe
588       end do
589 c
590 c     scale atomic friction coefficients by accessible area
591 c
592       if (lprn) write (iout,*) 
593      &  "Original gammas, surface areas, scaling factors, new gammas, ",
594      &  "std's of stochastic forces"
595       ind=0
596       do i=nnt,nct-1
597         if (radius(i).gt.0.0d0) then
598           call surfatom (i,area,radius)
599           ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
600           if (lprn) write (iout,'(i5,3f10.5,$)') 
601      &      i,gamvec(ind+1),area,ratio
602           do j=1,3
603             ind=ind+1
604             gamvec(ind) = ratio * gamvec(ind)
605           enddo
606           stdforcp(i)=stdfp*dsqrt(gamvec(ind))
607           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
608         endif
609       enddo
610       do i=nnt,nct
611         if (radius(i+nres).gt.0.0d0) then
612           call surfatom (i+nres,area,radius)
613           ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
614           if (lprn) write (iout,'(i5,3f10.5,$)') 
615      &      i,gamvec(ind+1),area,ratio
616           do j=1,3
617             ind=ind+1 
618             gamvec(ind) = ratio * gamvec(ind)
619           enddo
620           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
621           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
622         endif
623       enddo
624
625       return
626       end