1 subroutine friction_force
2 implicit real*8 (a-h,o-z)
9 include 'COMMON.INTERACT'
12 include 'COMMON.LANGEVIN'
14 include 'COMMON.LANGEVIN.lang0'
16 include 'COMMON.IOUNITS'
17 double precision gamvec(MAXRES6)
19 double precision vv(3),vvtot(3,maxres),v_work(MAXRES6),
20 & ginvfric(maxres2,maxres2)
21 common /przechowalnia/ ginvfric
23 logical lprn /.false./, checkmode /.false./
37 d_t_work(ind+j)=d_t(j,i)
42 if (itype(i).ne.10 .and. itype(i).ne.21) then
44 d_t_work(ind+j)=d_t(j,i+nres)
50 call fricmat_mult(d_t_work,fric_work)
52 if (.not.checkmode) return
55 write (iout,*) "d_t_work and fric_work"
57 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
61 friction(j,0)=fric_work(j)
66 friction(j,i)=fric_work(ind+j)
71 if (itype(i).ne.10 .and. itype(i).ne.21) then
73 friction(j,i+nres)=fric_work(ind+j)
79 write(iout,*) "Friction backbone"
81 write(iout,'(i5,3e15.5,5x,3e15.5)')
82 & i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
84 write(iout,*) "Friction side chain"
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)
96 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
97 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
101 write (iout,*) "vvtot backbone and sidechain"
103 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),
104 & (vvtot(j,i+nres),j=1,3)
109 v_work(ind+j)=vvtot(j,i)
115 v_work(ind+j)=vvtot(j,i+nres)
119 write (iout,*) "v_work gamvec and site-based friction forces"
121 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),
122 & gamvec(i)*v_work(i)
125 c fric_work1(i)=0.0d0
127 c fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
130 c write (iout,*) "fric_work and fric_work1"
132 c write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
138 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
142 write (iout,*) "ginvfric"
144 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
146 write (iout,*) "symmetry check"
149 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
155 c-----------------------------------------------------
156 subroutine stochastic_force(stochforcvec)
157 implicit real*8 (a-h,o-z)
163 include 'COMMON.CHAIN'
164 include 'COMMON.DERIV'
166 include 'COMMON.LOCAL'
167 include 'COMMON.INTERACT'
169 include 'COMMON.TIME1'
171 include 'COMMON.LANGEVIN'
173 include 'COMMON.LANGEVIN.lang0'
175 include 'COMMON.IOUNITS'
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./
193 c Compute the stochastic forces acting on bodies. Store in force.
199 force(j,i)=anorm_distr(x,sig,lowb,highb)
207 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
211 time_fsample=time_fsample+MPI_Wtime()-time00
213 time_fsample=time_fsample+tcpu()-time00
215 c Compute the stochastic forces acting on virtual-bond vectors.
221 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
224 ff(j)=ff(j)+force(j,i)
226 if (itype(i+1).ne.21) then
228 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
229 ff(j)=ff(j)+force(j,i+nres+1)
234 stochforc(j,0)=ff(j)+force(j,nnt+nres)
237 if (itype(i).ne.10 .and. itype(i).ne.21) then
239 stochforc(j,i+nres)=force(j,i+nres)
245 stochforcvec(j)=stochforc(j,0)
250 stochforcvec(ind+j)=stochforc(j,i)
255 if (itype(i).ne.10 .and. itype(i).ne.21) then
257 stochforcvec(ind+j)=stochforc(j,i+nres)
263 write (iout,*) "stochforcvec"
265 write(iout,'(i5,e15.5)') i,stochforcvec(i)
267 write(iout,*) "Stochastic forces backbone"
269 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
271 write(iout,*) "Stochastic forces side chain"
273 write(iout,'(i5,3e15.5)')
274 & i,(stochforc(j,i+nres),j=1,3)
284 forcvec(ind+j)=force(j,i)
291 forcvec(j+ind)=force(j,i+nres)
296 write (iout,*) "forcvec"
300 write (iout,'(2i3,2f10.5)') i,j,force(j,i),
307 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),
317 c------------------------------------------------------------------
318 subroutine setup_fricmat
319 implicit real*8 (a-h,o-z)
325 include 'COMMON.CHAIN'
326 include 'COMMON.DERIV'
328 include 'COMMON.LOCAL'
329 include 'COMMON.INTERACT'
331 include 'COMMON.SETUP'
332 include 'COMMON.TIME1'
333 c integer licznik /0/
336 include 'COMMON.LANGEVIN'
338 include 'COMMON.LANGEVIN.lang0'
340 include 'COMMON.IOUNITS'
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
351 if (fg_rank.ne.king) goto 10
353 c Zeroing out fricmat
359 c Load the friction coefficients corresponding to peptide groups
365 c Load the friction coefficients corresponding to side chains
373 gamvec(ii)=gamsc(iti)
375 if (surfarea) call sdarea(gamvec)
377 c write (iout,*) "Matrix A and vector gamma"
379 c write (iout,'(i2,$)') i
381 c write (iout,'(f4.1,$)') A(i,j)
383 c write (iout,'(f8.3)') gamvec(i)
387 write (iout,*) "Vector gamvec"
389 write (iout,'(i5,f10.5)') i, gamvec(i)
393 c The friction matrix
398 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
405 write (iout,'(//a)') "Matrix fricmat"
406 call matout2(dimen,dimen,maxres2,maxres2,fricmat)
408 if (lang.eq.2 .or. lang.eq.3) then
409 c Mass-scale the friction matrix if non-direct integration will be performed
415 Ginvfric(i,j)=Ginvfric(i,j)+
416 & Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
421 c Diagonalize the friction matrix
426 Ghalf(ind)=Ginvfric(i,j)
429 call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
432 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
433 & " mass-scaled friction matrix"
434 call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
436 c Precompute matrices for tinker stochastic integrator
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)
450 else if (lang.eq.4) then
451 c Diagonalize the friction matrix
456 Ghalf(ind)=fricmat(i,j)
459 call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
462 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
464 call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
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)
471 write (iout,*) "Number of zero eigenvalues:",nzero
476 fricmat(i,j)=fricmat(i,j)
477 & +fricvec(i,k)*fricvec(j,k)/fricgam(k)
482 write (iout,'(//a)') "Generalized inverse of fricmat"
483 call matout(dimen,dimen,maxres6,maxres6,fricmat)
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
496 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
497 time_Bcast=time_Bcast+MPI_Wtime()-time00
498 c print *,"Processor",myrank,
499 c & " BROADCAST iorder in SETUP_FRICMAT"
502 c write (iout,*) "setup_fricmat licznik",licznik
504 c Scatter the friction matrix
505 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
506 & nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
507 & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
508 time_scatter=time_scatter+MPI_Wtime()-time00
511 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
513 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
518 fricmat(j,i)=fcopy(i,j)
521 c write (iout,*) "My chunk of fricmat"
522 c call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
527 c-------------------------------------------------------------------------------
528 subroutine sdarea(gamvec)
530 c Scale the friction coefficients according to solvent accessible surface areas
531 c Code adapted from TINKER
534 implicit real*8 (a-h,o-z)
536 include 'COMMON.CONTROL'
540 include 'COMMON.LANGEVIN'
542 include 'COMMON.LANGEVIN.lang0'
544 include 'COMMON.CHAIN'
545 include 'COMMON.DERIV'
547 include 'COMMON.LOCAL'
548 include 'COMMON.INTERACT'
549 include 'COMMON.IOUNITS'
550 include 'COMMON.NAMES'
551 double precision radius(maxres2),gamvec(maxres2)
552 parameter (twosix=1.122462048309372981d0)
553 logical lprn /.true./
555 c determine new friction coefficients every few SD steps
557 c set the atomic radii to estimates of sigma values
559 c print *,"Entered sdarea"
565 c Load peptide group radii
569 c Load side chain radii
572 radius(i+nres)=restok(iti)
575 c write (iout,*) "i",i," radius",radius(i)
578 radius(i) = radius(i) / twosix
579 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
582 c scale atomic friction coefficients by accessible area
584 if (lprn) write (iout,*)
585 & "Original gammas, surface areas, scaling factors, new gammas, ",
586 & "std's of stochastic forces"
589 if (radius(i).gt.0.0d0) then
590 call surfatom (i,area,radius)
591 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
592 if (lprn) write (iout,'(i5,3f10.5,$)')
593 & i,gamvec(ind+1),area,ratio
596 gamvec(ind) = ratio * gamvec(ind)
598 stdforcp(i)=stdfp*dsqrt(gamvec(ind))
599 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
603 if (radius(i+nres).gt.0.0d0) then
604 call surfatom (i+nres,area,radius)
605 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
606 if (lprn) write (iout,'(i5,3f10.5,$)')
607 & i,gamvec(ind+1),area,ratio
610 gamvec(ind) = ratio * gamvec(ind)
612 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
613 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)