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 vvec(ind+1)=v_work(j,i)
76 if (iabs(itype(i)).ne.10) then
77 vvec(ind+1)=v_work(j,i+nres)
82 write (iout,*) "vvec ind",ind
83 write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
85 c write (iout,*) "chain",i," ind",ind," n",n
86 call fivediagmult(n,DMfric(iposc),DU1fric(iposc),
87 & DU2fric(iposc),vvec,rs)
89 fric_work(3*(i-1)+j)=-rs(i)
94 write (iout,*) "Vector fric_work"
95 write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
110 d_t_work(ind+j)=d_t(j,i)
115 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
117 d_t_work(ind+j)=d_t(j,i+nres)
123 call fricmat_mult(d_t_work,fric_work)
125 if (.not.checkmode) return
128 write (iout,*) "d_t_work and fric_work"
130 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
134 friction(j,0)=fric_work(j)
139 friction(j,i)=fric_work(ind+j)
144 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
146 friction(j,i+nres)=fric_work(ind+j)
152 write(iout,*) "Friction backbone"
154 write(iout,'(i5,3e15.5,5x,3e15.5)')
155 & i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
157 write(iout,*) "Friction side chain"
159 write(iout,'(i5,3e15.5,5x,3e15.5)')
160 & i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
169 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
170 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
174 write (iout,*) "vvtot backbone and sidechain"
176 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),
177 & (vvtot(j,i+nres),j=1,3)
182 v_work(ind+j)=vvtot(j,i)
188 v_work(ind+j)=vvtot(j,i+nres)
192 write (iout,*) "v_work gamvec and site-based friction forces"
194 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),
195 & gamvec(i)*v_work(i)
198 c fric_work1(i)=0.0d0
200 c fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
203 c write (iout,*) "fric_work and fric_work1"
205 c write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
211 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
215 write (iout,*) "ginvfric"
217 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
219 write (iout,*) "symmetry check"
222 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
229 c-----------------------------------------------------
230 subroutine stochastic_force(stochforcvec)
235 double precision time00
238 include 'COMMON.CHAIN'
239 include 'COMMON.DERIV'
241 include 'COMMON.LOCAL'
242 include 'COMMON.INTERACT'
245 include 'COMMON.LAGRANGE.5diag'
247 include 'COMMON.LAGRANGE'
249 include 'COMMON.TIME1'
251 include 'COMMON.LANGEVIN'
254 include 'COMMON.LANGEVIN.lang0.5diag'
256 include 'COMMON.LANGEVIN.lang0'
259 include 'COMMON.IOUNITS'
261 double precision x,sig,lowb,highb,
262 & ff(3),force(3,0:MAXRES2),zeta2,lowb2,
263 & highb2,sig2,forcvec(MAXRES6),stochforcvec(MAXRES6)
264 logical lprn /.false./
265 integer i,j,ind,ii,iti
266 double precision anorm_distr
268 integer ichain,innt,inct,iposc
283 c Compute the stochastic forces acting on bodies. Store in force.
289 force(j,i)=anorm_distr(x,sig,lowb,highb)
297 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
301 write (iout,*) "Stochastic forces on sites"
303 write (iout,'(i5,2(3f10.5,5x))') i,(force(j,i),j=1,3),
304 & (force(j,i+nres),j=1,3)
308 time_fsample=time_fsample+MPI_Wtime()-time00
310 time_fsample=time_fsample+tcpu()-time00
315 innt=chain_border(1,ichain)
316 inct=chain_border(2,ichain)
317 iposc=iposd_chain(ichain)
318 c write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
319 c & " inct",inct," iposc",iposc
321 stochforcvec(ind+j)=0.5d0*force(j,innt)
323 if (iabs(itype(innt)).eq.10) then
325 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres)
331 stochforcvec(ind+j)=force(j,innt+nres)
337 stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1))
339 if (iabs(itype(i)).eq.10) then
341 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres)
347 stochforcvec(ind+j)=force(j,i+nres)
353 stochforcvec(ind+j)=0.5d0*force(j,inct-1)
355 if (iabs(itype(inct)).eq.10) then
357 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres)
363 stochforcvec(ind+j)=force(j,inct+nres)
367 c write (iout,*) "chain",ichain," ind",ind
370 write (iout,*) "stochforcvec"
371 write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind)
374 c Compute the stochastic forces acting on virtual-bond vectors.
380 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
383 ff(j)=ff(j)+force(j,i)
385 if (itype(i+1).ne.ntyp1) then
387 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
388 ff(j)=ff(j)+force(j,i+nres+1)
393 stochforc(j,0)=ff(j)+force(j,nnt+nres)
396 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
398 stochforc(j,i+nres)=force(j,i+nres)
403 stochforcvec(j)=stochforc(j,0)
408 stochforcvec(ind+j)=stochforc(j,i)
413 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
415 stochforcvec(ind+j)=stochforc(j,i+nres)
422 write (iout,*) "stochforcvec"
424 write(iout,'(i5,e15.5)') i,stochforcvec(i)
426 write(iout,*) "Stochastic forces backbone"
428 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
430 write(iout,*) "Stochastic forces side chain"
432 write(iout,'(i5,3e15.5)')
433 & i,(stochforc(j,i+nres),j=1,3)
443 forcvec(ind+j)=force(j,i)
450 forcvec(j+ind)=force(j,i+nres)
455 write (iout,*) "forcvec"
459 write (iout,'(2i3,2f10.5)') i,j,force(j,i),
466 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),
475 c------------------------------------------------------------------
476 subroutine setup_fricmat
481 double precision time00
485 include 'COMMON.CHAIN'
486 include 'COMMON.DERIV'
488 include 'COMMON.LOCAL'
489 include 'COMMON.INTERACT'
492 include 'COMMON.LAGRANGE.5diag'
494 include 'COMMON.LAGRANGE'
496 include 'COMMON.SETUP'
497 include 'COMMON.TIME1'
498 c integer licznik /0/
501 include 'COMMON.LANGEVIN'
504 include 'COMMON.LANGEVIN.lang0.5diag'
506 include 'COMMON.LANGEVIN.lang0'
509 include 'COMMON.IOUNITS'
511 integer i,j,k,l,ind,ind1,m,ii,iti,it,nzero,innt,inct
513 logical lprn /.false./
514 double precision dtdi,gamvec(MAXRES2)
515 common /syfek/ gamvec
517 double precision ginvfric(maxres2,maxres2),Ghalf(mmaxres2),
518 & fcopy(maxres2,maxres2)
519 double precision work(8*maxres2)
520 integer iwork(maxres2)
521 common /przechowalnia/ ginvfric,Ghalf,fcopy
525 if (fg_rank.ne.king) goto 10
532 c Load the friction coefficients corresponding to side chains
534 if (lprn) write (iout,*) "m",m
541 gamvec(ii)=gamsc(iabs(iti))
543 if (surfarea) call sdarea(gamvec)
545 write (iout,*) "Vector gamvec"
547 write (iout,'(i5,f10.5)') i, gamvec(i)
556 innt=chain_border(1,ichain)
557 inct=chain_border(2,ichain)
558 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
560 DMfric(ind)=gamvec(innt-nnt+1)/4
561 if (iabs(itype(innt)).eq.10) then
562 DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1)
565 DMfric(ind+1)=gamvec(m+innt-nnt+1)
568 c write (iout,*) "DMfric init ind",ind
571 DMfric(ind)=gamvec(i-nnt+1)/2
572 if (iabs(itype(i)).eq.10) then
573 DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1)
576 DMfric(ind+1)=gamvec(m+i-nnt+1)
580 c write (iout,*) "DMfric endloop ind",ind
581 if (inct.gt.innt) then
582 DMfric(ind)=gamvec(inct-1-nnt+1)/4
583 if (iabs(itype(inct)).eq.10) then
584 DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1)
587 DMfric(ind+1)=gamvec(inct+m-nnt+1)
591 c write (iout,*) "DMfric end ind",ind
595 ind=iposd_chain(ichain)
596 innt=chain_border(1,ichain)
597 inct=chain_border(2,ichain)
599 if (iabs(itype(i)).ne.10) then
602 DU1fric(ind)=gamvec(i-nnt+1)/4
609 ind=iposd_chain(ichain)
610 innt=chain_border(1,ichain)
611 inct=chain_border(2,ichain)
613 if (iabs(itype(i)).ne.10) then
614 DU2fric(ind)=gamvec(i-nnt+1)/4
624 write(iout,*)"The upper part of the five-diagonal friction matrix"
626 write (iout,'(a,i5)') 'Chain',ichain
627 innt=iposd_chain(ichain)
628 inct=iposd_chain(ichain)+dimen_chain(ichain)-1
630 if (i.lt.inct-1) then
631 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),
633 else if (i.eq.inct-1) then
634 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i)
636 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i)
643 c The friction matrix
648 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
655 write (iout,'(//a)') "Matrix fricmat"
656 call matout2(dimen,dimen,maxres2,maxres2,fricmat)
658 if (lang.eq.2 .or. lang.eq.3) then
659 c Mass-scale the friction matrix if non-direct integration will be performed
665 Ginvfric(i,j)=Ginvfric(i,j)+
666 & Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
671 c Diagonalize the friction matrix
676 Ghalf(ind)=Ginvfric(i,j)
679 call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
682 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
683 & " mass-scaled friction matrix"
684 call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
686 c Precompute matrices for tinker stochastic integrator
693 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
694 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
700 else if (lang.eq.4) then
701 c Diagonalize the friction matrix
706 Ghalf(ind)=fricmat(i,j)
709 call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
712 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
714 call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
716 c Determine the number of zero eigenvalues of the friction matrix
717 nzero=max0(dimen-dimen1,0)
718 c do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
721 write (iout,*) "Number of zero eigenvalues:",nzero
726 fricmat(i,j)=fricmat(i,j)
727 & +fricvec(i,k)*fricvec(j,k)/fricgam(k)
732 write (iout,'(//a)') "Generalized inverse of fricmat"
733 call matout(dimen,dimen,maxres6,maxres6,fricmat)
738 if (nfgtasks.gt.1) then
739 if (fg_rank.eq.0) then
740 c The matching BROADCAST for fg processors is called in ERGASTULUM
746 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
748 time_Bcast=time_Bcast+MPI_Wtime()-time00
750 time_Bcast=time_Bcast+tcpu()-time00
752 c print *,"Processor",myrank,
753 c & " BROADCAST iorder in SETUP_FRICMAT"
756 c write (iout,*) "setup_fricmat licznik",licznik
762 c Scatter the friction matrix
763 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
764 & nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
765 & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
768 time_scatter=time_scatter+MPI_Wtime()-time00
769 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
771 time_scatter=time_scatter+tcpu()-time00
772 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
777 fricmat(j,i)=fcopy(i,j)
780 c write (iout,*) "My chunk of fricmat"
781 c call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
787 c-------------------------------------------------------------------------------
788 subroutine sdarea(gamvec)
790 c Scale the friction coefficients according to solvent accessible surface areas
791 c Code adapted from TINKER
796 include 'COMMON.CONTROL'
800 include 'COMMON.LANGEVIN'
803 include 'COMMON.LANGEVIN.lang0.5diag'
805 include 'COMMON.LANGEVIN.lang0'
808 include 'COMMON.CHAIN'
809 include 'COMMON.DERIV'
811 include 'COMMON.LOCAL'
812 include 'COMMON.INTERACT'
813 include 'COMMON.IOUNITS'
814 include 'COMMON.NAMES'
815 double precision radius(maxres2),gamvec(maxres2)
816 double precision twosix
817 parameter (twosix=1.122462048309372981d0)
818 logical lprn /.false./
820 double precision probe,area,ratio
822 c determine new friction coefficients every few SD steps
824 c set the atomic radii to estimates of sigma values
826 c print *,"Entered sdarea"
832 c Load peptide group radii
836 c Load side chain radii
839 if (iti.ne.ntyp1) radius(i+nres)=restok(iti)
842 c write (iout,*) "i",i," radius",radius(i)
845 radius(i) = radius(i) / twosix
846 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
849 c scale atomic friction coefficients by accessible area
851 if (lprn) write (iout,*)
852 & "Original gammas, surface areas, scaling factors, new gammas, ",
853 & "std's of stochastic forces"
856 if (radius(i).gt.0.0d0) then
857 call surfatom (i,area,radius)
858 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
859 if (lprn) write (iout,'(i5,3f10.5,$)')
860 & i,gamvec(ind+1),area,ratio
863 gamvec(ind) = ratio * gamvec(ind)
865 stdforcp(i)=stdfp*dsqrt(gamvec(ind))
866 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
870 if (radius(i+nres).gt.0.0d0) then
871 call surfatom (i+nres,area,radius)
872 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
873 if (lprn) write (iout,'(i5,3f10.5,$)')
874 & i,gamvec(ind+1),area,ratio
877 gamvec(ind) = ratio * gamvec(ind)
879 if (itype(i).ne.ntyp1)
880 & stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
881 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)