1 subroutine friction_force
9 include 'COMMON.INTERACT'
12 include 'COMMON.LAGRANGE.5diag'
14 include 'COMMON.LAGRANGE'
17 include 'COMMON.LANGEVIN'
20 include 'COMMON.LANGEVIN.lang0.5diag'
22 include 'COMMON.LANGEVIN.lang0'
25 include 'COMMON.IOUNITS'
27 integer iposc,ichain,n,innt,inct
28 double precision v_work(3,maxres2),vvec(maxres2_chain),rs(maxres2)
30 double precision gamvec(MAXRES6)
32 double precision vv(3),vvtot(3,maxres),v_work(MAXRES6),
33 & ginvfric(maxres2,maxres2)
34 common /przechowalnia/ ginvfric
38 logical lprn /.false./, checkmode /.false./
40 c Here accelerations due to friction forces are computed right after forces.
44 v_work(j,nnt)=d_t(j,0)
48 v_work(j,i)=v_work(j,i-1)+d_t(j,i-1)
52 if (iabs(itype(i)).ne.10 .and. iabs(itype(i)).ne.ntyp1) then
54 v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres)
59 write (iout,*) "v_work"
61 write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3)
68 iposc=iposd_chain(ichain)
69 c write (iout,*) "friction_force j",j," ichain",ichain,
70 c & " n",n," iposc",iposc,iposc+n-1
71 innt=chain_border(1,ichain)
72 inct=chain_border(2,ichain)
74 c innt=chain_border(1,1)
75 c inct=chain_border(2,1)
77 vvec(ind+1)=v_work(j,i)
79 if (iabs(itype(i)).ne.10) then
80 vvec(ind+1)=v_work(j,i+nres)
85 write (iout,*) "vvec ind",ind," n",n
86 write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
88 c write (iout,*) "chain",i," ind",ind," n",n
89 call fivediagmult(n,DMfric(iposc),DU1fric(iposc),
90 & DU2fric(iposc),vvec(iposc),rs)
93 write (iout,'(f10.5)') (rs(i),i=1,n)
96 c write (iout,*) "ichain",ichain," i",i," j",j,
97 c & "index",3*(i-1)+j,"rs",rs(i-iposc+1)
98 fric_work(3*(i-1)+j)=-rs(i-iposc+1)
103 write (iout,*) "Vector fric_work dimen3",dimen3
104 write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
119 d_t_work(ind+j)=d_t(j,i)
124 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
126 d_t_work(ind+j)=d_t(j,i+nres)
132 call fricmat_mult(d_t_work,fric_work)
134 if (.not.checkmode) return
137 write (iout,*) "d_t_work and fric_work"
139 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
143 friction(j,0)=fric_work(j)
148 friction(j,i)=fric_work(ind+j)
153 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
155 friction(j,i+nres)=fric_work(ind+j)
161 write(iout,*) "Friction backbone"
163 write(iout,'(i5,3e15.5,5x,3e15.5)')
164 & i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
166 write(iout,*) "Friction side chain"
168 write(iout,'(i5,3e15.5,5x,3e15.5)')
169 & i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
178 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
179 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
183 write (iout,*) "vvtot backbone and sidechain"
185 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),
186 & (vvtot(j,i+nres),j=1,3)
191 v_work(ind+j)=vvtot(j,i)
197 v_work(ind+j)=vvtot(j,i+nres)
201 write (iout,*) "v_work gamvec and site-based friction forces"
203 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),
204 & gamvec(i)*v_work(i)
207 c fric_work1(i)=0.0d0
209 c fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
212 c write (iout,*) "fric_work and fric_work1"
214 c write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
220 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
224 write (iout,*) "ginvfric"
226 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
228 write (iout,*) "symmetry check"
231 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
238 c-----------------------------------------------------
239 subroutine stochastic_force(stochforcvec)
244 double precision time00
247 include 'COMMON.CHAIN'
248 include 'COMMON.DERIV'
250 include 'COMMON.LOCAL'
251 include 'COMMON.INTERACT'
254 include 'COMMON.LAGRANGE.5diag'
256 include 'COMMON.LAGRANGE'
258 include 'COMMON.TIME1'
260 include 'COMMON.LANGEVIN'
263 include 'COMMON.LANGEVIN.lang0.5diag'
265 include 'COMMON.LANGEVIN.lang0'
268 include 'COMMON.IOUNITS'
270 double precision x,sig,lowb,highb,
271 & ff(3),force(3,0:MAXRES2),zeta2,lowb2,
272 & highb2,sig2,forcvec(MAXRES6),stochforcvec(MAXRES6)
273 logical lprn /.false./
274 integer i,j,ind,ii,iti
275 double precision anorm_distr
277 integer ichain,innt,inct,iposc
292 c Compute the stochastic forces acting on bodies. Store in force.
295 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
301 force(j,i)=anorm_distr(x,sig,lowb,highb)
309 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
313 write (iout,*) "Stochastic forces on sites"
315 write (iout,'(i5,2(3f10.5,5x))') i,(force(j,i),j=1,3),
316 & (force(j,i+nres),j=1,3)
320 time_fsample=time_fsample+MPI_Wtime()-time00
322 time_fsample=time_fsample+tcpu()-time00
327 innt=chain_border(1,ichain)
328 inct=chain_border(2,ichain)
329 iposc=iposd_chain(ichain)
331 c innt=chain_border(1,1)
332 c inct=chain_border(2,1)
333 c iposc=iposd_chain(1)
334 c write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
335 c & " inct",inct," iposc",iposc
337 stochforcvec(ind+j)=0.5d0*force(j,innt)
339 if (iabs(itype(innt)).eq.10) then
341 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres)
347 stochforcvec(ind+j)=force(j,innt+nres)
353 stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1))
355 if (iabs(itype(i)).eq.10) then
357 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres)
363 stochforcvec(ind+j)=force(j,i+nres)
369 stochforcvec(ind+j)=0.5d0*force(j,inct-1)
371 if (iabs(itype(inct)).eq.10) then
373 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres)
379 stochforcvec(ind+j)=force(j,inct+nres)
383 c write (iout,*) "chain",ichain," ind",ind
386 write (iout,*) "stochforcvec"
387 write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind)
390 c Compute the stochastic forces acting on virtual-bond vectors.
396 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
399 ff(j)=ff(j)+force(j,i)
401 if (itype(i+1).ne.ntyp1) then
403 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
404 ff(j)=ff(j)+force(j,i+nres+1)
409 stochforc(j,0)=ff(j)+force(j,nnt+nres)
412 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
414 stochforc(j,i+nres)=force(j,i+nres)
419 stochforcvec(j)=stochforc(j,0)
424 stochforcvec(ind+j)=stochforc(j,i)
429 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
431 stochforcvec(ind+j)=stochforc(j,i+nres)
437 write (iout,*) "stochforcvec"
439 write(iout,'(i5,e15.5)') i,stochforcvec(i)
441 write(iout,*) "Stochastic forces backbone"
443 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
445 write(iout,*) "Stochastic forces side chain"
447 write(iout,'(i5,3e15.5)')
448 & i,(stochforc(j,i+nres),j=1,3)
458 forcvec(ind+j)=force(j,i)
465 forcvec(j+ind)=force(j,i+nres)
470 write (iout,*) "forcvec"
474 write (iout,'(2i3,2f10.5)') i,j,force(j,i),
481 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),
490 c------------------------------------------------------------------
491 subroutine setup_fricmat
496 double precision time00
500 include 'COMMON.CHAIN'
501 include 'COMMON.DERIV'
503 include 'COMMON.LOCAL'
504 include 'COMMON.INTERACT'
507 include 'COMMON.LAGRANGE.5diag'
509 include 'COMMON.LAGRANGE'
511 include 'COMMON.SETUP'
512 include 'COMMON.TIME1'
513 c integer licznik /0/
516 include 'COMMON.LANGEVIN'
519 include 'COMMON.LANGEVIN.lang0.5diag'
521 include 'COMMON.LANGEVIN.lang0'
524 include 'COMMON.IOUNITS'
526 integer i,j,k,l,ind,ind1,m,ii,iti,it,nzero,innt,inct
528 logical lprn /.false./
529 double precision dtdi,gamvec(MAXRES2)
530 common /syfek/ gamvec
532 double precision ginvfric(maxres2,maxres2),Ghalf(mmaxres2),
533 & fcopy(maxres2,maxres2)
534 double precision work(8*maxres2)
535 integer iwork(maxres2)
536 common /przechowalnia/ ginvfric,Ghalf,fcopy
540 if (fg_rank.ne.king) goto 10
547 c Load the friction coefficients corresponding to side chains
549 if (lprn) write (iout,*) "m",m
556 gamvec(ii)=gamsc(iabs(iti))
558 if (surfarea) call sdarea(gamvec)
560 write (iout,*) "Vector gamvec ii",ii
562 write (iout,'(i5,f10.5)') i, gamvec(i)
571 innt=chain_border(1,ichain)
572 inct=chain_border(2,ichain)
573 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
575 DMfric(ind)=gamvec(innt-nnt+1)/4
576 if (iabs(itype(innt)).eq.10) then
577 DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1)
580 DMfric(ind+1)=gamvec(m+innt-nnt+1)
583 c write (iout,*) "DMfric init ind",ind
586 DMfric(ind)=gamvec(i-nnt+1)/2
587 if (iabs(itype(i)).eq.10) then
588 DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1)
591 DMfric(ind+1)=gamvec(m+i-nnt+1)
595 c write (iout,*) "DMfric endloop ind",ind
596 if (inct.gt.innt) then
597 DMfric(ind)=gamvec(inct-1-nnt+1)/4
598 if (iabs(itype(inct)).eq.10) then
599 DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1)
602 DMfric(ind+1)=gamvec(inct+m-nnt+1)
606 c write (iout,*) "DMfric end ind",ind
610 ind=iposd_chain(ichain)
611 innt=chain_border(1,ichain)
612 inct=chain_border(2,ichain)
614 if (iabs(itype(i)).ne.10) then
617 DU1fric(ind)=gamvec(i-nnt+1)/4
624 ind=iposd_chain(ichain)
625 innt=chain_border(1,ichain)
626 inct=chain_border(2,ichain)
628 if (iabs(itype(i)).ne.10) then
629 DU2fric(ind)=gamvec(i-nnt+1)/4
639 write(iout,*)"The upper part of the five-diagonal friction matrix"
641 write (iout,'(a,i5)') 'Chain',ichain
642 innt=iposd_chain(ichain)
643 inct=iposd_chain(ichain)+dimen_chain(ichain)-1
645 if (i.lt.inct-1) then
646 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),
648 else if (i.eq.inct-1) then
649 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i)
651 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i)
658 c The friction matrix
663 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
670 write (iout,'(//a)') "Matrix fricmat"
671 call matout2(dimen,dimen,maxres2,maxres2,fricmat)
673 if (lang.eq.2 .or. lang.eq.3) then
674 c Mass-scale the friction matrix if non-direct integration will be performed
680 Ginvfric(i,j)=Ginvfric(i,j)+
681 & Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
686 c Diagonalize the friction matrix
691 Ghalf(ind)=Ginvfric(i,j)
694 call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
697 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
698 & " mass-scaled friction matrix"
699 call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
701 c Precompute matrices for tinker stochastic integrator
708 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
709 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
715 else if (lang.eq.4) then
716 c Diagonalize the friction matrix
721 Ghalf(ind)=fricmat(i,j)
724 call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
727 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
729 call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
731 c Determine the number of zero eigenvalues of the friction matrix
732 nzero=max0(dimen-dimen1,0)
733 c do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
736 write (iout,*) "Number of zero eigenvalues:",nzero
741 fricmat(i,j)=fricmat(i,j)
742 & +fricvec(i,k)*fricvec(j,k)/fricgam(k)
747 write (iout,'(//a)') "Generalized inverse of fricmat"
748 call matout(dimen,dimen,maxres6,maxres6,fricmat)
753 if (nfgtasks.gt.1) then
754 if (fg_rank.eq.0) then
755 c The matching BROADCAST for fg processors is called in ERGASTULUM
761 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
763 time_Bcast=time_Bcast+MPI_Wtime()-time00
765 time_Bcast=time_Bcast+tcpu()-time00
767 c print *,"Processor",myrank,
768 c & " BROADCAST iorder in SETUP_FRICMAT"
771 c write (iout,*) "setup_fricmat licznik",licznik
777 c Scatter the friction matrix
778 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
779 & nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
780 & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
783 time_scatter=time_scatter+MPI_Wtime()-time00
784 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
786 time_scatter=time_scatter+tcpu()-time00
787 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
792 fricmat(j,i)=fcopy(i,j)
795 c write (iout,*) "My chunk of fricmat"
796 c call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
802 c-------------------------------------------------------------------------------
803 subroutine sdarea(gamvec)
805 c Scale the friction coefficients according to solvent accessible surface areas
806 c Code adapted from TINKER
811 include 'COMMON.CONTROL'
815 include 'COMMON.LANGEVIN'
818 include 'COMMON.LANGEVIN.lang0.5diag'
820 include 'COMMON.LANGEVIN.lang0'
823 include 'COMMON.CHAIN'
824 include 'COMMON.DERIV'
826 include 'COMMON.LOCAL'
827 include 'COMMON.INTERACT'
828 include 'COMMON.IOUNITS'
829 include 'COMMON.NAMES'
830 double precision radius(maxres2),gamvec(maxres2)
831 double precision twosix
832 parameter (twosix=1.122462048309372981d0)
833 logical lprn /.false./
835 double precision probe,area,ratio
837 c determine new friction coefficients every few SD steps
839 c set the atomic radii to estimates of sigma values
841 c print *,"Entered sdarea"
847 c Load peptide group radii
851 c Load side chain radii
854 if (iti.ne.ntyp1) radius(i+nres)=restok(iti)
857 c write (iout,*) "i",i," radius",radius(i)
860 radius(i) = radius(i) / twosix
861 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
864 c scale atomic friction coefficients by accessible area
866 if (lprn) write (iout,*)
867 & "Original gammas, surface areas, scaling factors, new gammas, ",
868 & "std's of stochastic forces"
871 if (radius(i).gt.0.0d0) then
872 call surfatom (i,area,radius)
873 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
874 if (lprn) write (iout,'(i5,3f10.5,$)')
875 & i,gamvec(ind+1),area,ratio
878 gamvec(ind) = ratio * gamvec(ind)
880 stdforcp(i)=stdfp*dsqrt(gamvec(ind))
881 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
885 if (radius(i+nres).gt.0.0d0) then
886 call surfatom (i+nres,area,radius)
887 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
888 if (lprn) write (iout,'(i5,3f10.5,$)')
889 & i,gamvec(ind+1),area,ratio
892 gamvec(ind) = ratio * gamvec(ind)
894 if (itype(i).ne.ntyp1)
895 & stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
896 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)