1 subroutine friction_force
11 include 'COMMON.LOCAL'
12 include 'COMMON.INTERACT'
15 include 'COMMON.LAGRANGE.5diag'
17 include 'COMMON.LAGRANGE'
20 include 'COMMON.LANGEVIN'
23 include 'COMMON.LANGEVIN.lang0.5diag'
25 include 'COMMON.LANGEVIN.lang0'
28 include 'COMMON.IOUNITS'
30 integer iposc,ichain,n,innt,inct
31 double precision v_work(3,maxres2),vvec(maxres2),rs(maxres2)
33 double precision gamvec(MAXRES6)
35 double precision vv(3),vvtot(3,maxres),v_work(MAXRES6),
36 & ginvfric(maxres2,maxres2)
37 common /przechowalnia/ ginvfric
41 logical lprn /.false./, checkmode /.false./
44 include 'COMMON.TIME1'
45 double precision time01
47 c Here accelerations due to friction forces are computed right after forces.
48 d_t_work(:6*nres)=0.0d0
51 v_work(j,nnt)=d_t(j,0)
55 v_work(j,i)=v_work(j,i-1)+d_t(j,i-1)
59 if (iabs(itype(i)).ne.10 .and. iabs(itype(i)).ne.ntyp1) then
61 v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres)
66 write (iout,*) "v_work"
68 write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3)
75 iposc=iposd_chain(ichain)
76 c write (iout,*) "friction_force j",j," ichain",ichain,
77 c & " n",n," iposc",iposc,iposc+n-1
78 innt=chain_border(1,ichain)
79 inct=chain_border(2,ichain)
81 c innt=chain_border(1,1)
82 c inct=chain_border(2,1)
84 vvec(ind+1)=v_work(j,i)
86 if (iabs(itype(i)).ne.10) then
87 vvec(ind+1)=v_work(j,i+nres)
92 write (iout,*) "vvec ind",ind," n",n
93 write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
95 c write (iout,*) "chain",i," ind",ind," n",n
99 call fivediagmult(n,DMfric(iposc),DU1fric(iposc),
100 & DU2fric(iposc),vvec(iposc),rs)
102 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
106 write (iout,'(f10.5)') (rs(i),i=1,n)
109 c write (iout,*) "ichain",ichain," i",i," j",j,
110 c & "index",3*(i-1)+j,"rs",rs(i-iposc+1)
111 fric_work(3*(i-1)+j)=-rs(i-iposc+1)
116 write (iout,*) "Vector fric_work dimen3",dimen3
117 write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
132 d_t_work(ind+j)=d_t(j,i)
137 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
139 d_t_work(ind+j)=d_t(j,i+nres)
145 call fricmat_mult(d_t_work,fric_work)
147 if (.not.checkmode) return
150 write (iout,*) "d_t_work and fric_work"
152 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
156 friction(j,0)=fric_work(j)
161 friction(j,i)=fric_work(ind+j)
166 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
168 friction(j,i+nres)=fric_work(ind+j)
174 write(iout,*) "Friction backbone"
176 write(iout,'(i5,3e15.5,5x,3e15.5)')
177 & i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
179 write(iout,*) "Friction side chain"
181 write(iout,'(i5,3e15.5,5x,3e15.5)')
182 & i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
191 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
192 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
196 write (iout,*) "vvtot backbone and sidechain"
198 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),
199 & (vvtot(j,i+nres),j=1,3)
204 v_work(ind+j)=vvtot(j,i)
210 v_work(ind+j)=vvtot(j,i+nres)
214 write (iout,*) "v_work gamvec and site-based friction forces"
216 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),
217 & gamvec(i)*v_work(i)
220 c fric_work1(i)=0.0d0
222 c fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
225 c write (iout,*) "fric_work and fric_work1"
227 c write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
233 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
237 write (iout,*) "ginvfric"
239 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
241 write (iout,*) "symmetry check"
244 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
251 c-----------------------------------------------------
252 subroutine stochastic_force(stochforcvec)
257 double precision time00
260 include 'COMMON.CHAIN'
261 include 'COMMON.DERIV'
263 include 'COMMON.LOCAL'
264 include 'COMMON.INTERACT'
267 include 'COMMON.LAGRANGE.5diag'
269 include 'COMMON.LAGRANGE'
271 include 'COMMON.TIME1'
273 include 'COMMON.LANGEVIN'
276 include 'COMMON.LANGEVIN.lang0.5diag'
278 include 'COMMON.LANGEVIN.lang0'
281 include 'COMMON.IOUNITS'
283 double precision x,sig,lowb,highb,
284 & ff(3),force(3,0:MAXRES2),zeta2,lowb2,
285 & highb2,sig2,forcvec(MAXRES6),stochforcvec(MAXRES6)
286 logical lprn /.false./
287 integer i,j,ind,ii,iti
288 double precision anorm_distr
290 integer ichain,innt,inct,iposc
305 c Compute the stochastic forces acting on bodies. Store in force.
308 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
314 force(j,i)=anorm_distr(x,sig,lowb,highb)
322 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
326 write (iout,*) "Stochastic forces on sites"
328 write (iout,'(i5,2(3f10.5,5x))') i,(force(j,i),j=1,3),
329 & (force(j,i+nres),j=1,3)
333 time_fsample=time_fsample+MPI_Wtime()-time00
335 time_fsample=time_fsample+tcpu()-time00
340 innt=chain_border(1,ichain)
341 inct=chain_border(2,ichain)
342 iposc=iposd_chain(ichain)
344 c innt=chain_border(1,1)
345 c inct=chain_border(2,1)
346 c iposc=iposd_chain(1)
347 c write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
348 c & " inct",inct," iposc",iposc
350 stochforcvec(ind+j)=0.5d0*force(j,innt)
352 if (iabs(itype(innt)).eq.10) then
354 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres)
360 stochforcvec(ind+j)=force(j,innt+nres)
366 stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1))
368 if (iabs(itype(i)).eq.10) then
370 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres)
376 stochforcvec(ind+j)=force(j,i+nres)
382 stochforcvec(ind+j)=0.5d0*force(j,inct-1)
384 if (iabs(itype(inct)).eq.10) then
386 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres)
392 stochforcvec(ind+j)=force(j,inct+nres)
396 c write (iout,*) "chain",ichain," ind",ind
399 write (iout,*) "stochforcvec"
400 write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind)
403 c Compute the stochastic forces acting on virtual-bond vectors.
409 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
412 ff(j)=ff(j)+force(j,i)
414 if (itype(i+1).ne.ntyp1) then
416 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
417 ff(j)=ff(j)+force(j,i+nres+1)
422 stochforc(j,0)=ff(j)+force(j,nnt+nres)
425 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
427 stochforc(j,i+nres)=force(j,i+nres)
432 stochforcvec(j)=stochforc(j,0)
437 stochforcvec(ind+j)=stochforc(j,i)
442 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
444 stochforcvec(ind+j)=stochforc(j,i+nres)
450 write (iout,*) "stochforcvec"
452 write(iout,'(i5,e15.5)') i,stochforcvec(i)
454 write(iout,*) "Stochastic forces backbone"
456 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
458 write(iout,*) "Stochastic forces side chain"
460 write(iout,'(i5,3e15.5)')
461 & i,(stochforc(j,i+nres),j=1,3)
471 forcvec(ind+j)=force(j,i)
478 forcvec(j+ind)=force(j,i+nres)
483 write (iout,*) "forcvec"
487 write (iout,'(2i3,2f10.5)') i,j,force(j,i),
494 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),
503 c------------------------------------------------------------------
504 subroutine setup_fricmat
509 double precision time00
514 parameter (mmaxres2=(maxres2*(maxres2+1)/2))
517 include 'COMMON.CHAIN'
518 include 'COMMON.DERIV'
520 include 'COMMON.LOCAL'
521 include 'COMMON.INTERACT'
524 include 'COMMON.LAGRANGE.5diag'
526 include 'COMMON.LAGRANGE'
528 include 'COMMON.SETUP'
529 include 'COMMON.TIME1'
530 c integer licznik /0/
533 include 'COMMON.LANGEVIN'
536 include 'COMMON.LANGEVIN.lang0.5diag'
538 include 'COMMON.LANGEVIN.lang0'
541 include 'COMMON.IOUNITS'
543 integer i,j,k,l,ind,ind1,m,ii,iti,it,nzero,innt,inct
545 logical lprn /.false./
546 double precision dtdi,gamvec(MAXRES2)
547 common /syfek/ gamvec
549 double precision ginvfric(maxres2,maxres2),Ghalf(mmaxres2),
550 & fcopy(maxres2,maxres2)
551 double precision work(8*maxres2)
552 integer iwork(maxres2)
553 common /przechowalnia/ ginvfric,Ghalf,fcopy
557 if (fg_rank.ne.king) goto 10
564 c Load the friction coefficients corresponding to side chains
566 if (lprn) write (iout,*) "m",m
573 gamvec(ii)=gamsc(iabs(iti))
575 if (surfarea) call sdarea(gamvec)
577 write (iout,*) "Vector gamvec ii",ii
579 write (iout,'(i5,f10.5)') i, gamvec(i)
583 DMfric(:2*nres)=0.0d0
584 DU1fric(:2*nres)=0.0d0
585 DU2fric(:2*nres)=0.0d0
588 innt=chain_border(1,ichain)
589 inct=chain_border(2,ichain)
590 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
592 DMfric(ind)=gamvec(innt-nnt+1)/4
593 if (iabs(itype(innt)).eq.10) then
594 DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1)
597 DMfric(ind+1)=gamvec(m+innt-nnt+1)
600 c write (iout,*) "DMfric init ind",ind
603 DMfric(ind)=gamvec(i-nnt+1)/2
604 if (iabs(itype(i)).eq.10) then
605 DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1)
608 DMfric(ind+1)=gamvec(m+i-nnt+1)
612 c write (iout,*) "DMfric endloop ind",ind
613 if (inct.gt.innt) then
614 DMfric(ind)=gamvec(inct-1-nnt+1)/4
615 if (iabs(itype(inct)).eq.10) then
616 DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1)
619 DMfric(ind+1)=gamvec(inct+m-nnt+1)
623 c write (iout,*) "DMfric end ind",ind
627 ind=iposd_chain(ichain)
628 innt=chain_border(1,ichain)
629 inct=chain_border(2,ichain)
631 if (iabs(itype(i)).ne.10) then
634 DU1fric(ind)=gamvec(i-nnt+1)/4
641 ind=iposd_chain(ichain)
642 innt=chain_border(1,ichain)
643 inct=chain_border(2,ichain)
645 if (iabs(itype(i)).ne.10) then
646 DU2fric(ind)=gamvec(i-nnt+1)/4
656 write(iout,*)"The upper part of the five-diagonal friction matrix"
658 write (iout,'(a,i5)') 'Chain',ichain
659 innt=iposd_chain(ichain)
660 inct=iposd_chain(ichain)+dimen_chain(ichain)-1
662 if (i.lt.inct-1) then
663 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),
665 else if (i.eq.inct-1) then
666 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i)
668 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i)
675 c The friction matrix
680 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
687 write (iout,'(//a)') "Matrix fricmat"
688 call matout2(dimen,dimen,maxres2,maxres2,fricmat)
690 if (lang.eq.2 .or. lang.eq.3) then
691 c Mass-scale the friction matrix if non-direct integration will be performed
697 Ginvfric(i,j)=Ginvfric(i,j)+
698 & Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
703 c Diagonalize the friction matrix
708 Ghalf(ind)=Ginvfric(i,j)
711 call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
714 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
715 & " mass-scaled friction matrix"
716 call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
718 c Precompute matrices for tinker stochastic integrator
725 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
726 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
732 else if (lang.eq.4) then
733 c Diagonalize the friction matrix
738 Ghalf(ind)=fricmat(i,j)
741 call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
744 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
746 call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
748 c Determine the number of zero eigenvalues of the friction matrix
749 nzero=max0(dimen-dimen1,0)
750 c do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
753 write (iout,*) "Number of zero eigenvalues:",nzero
758 fricmat(i,j)=fricmat(i,j)
759 & +fricvec(i,k)*fricvec(j,k)/fricgam(k)
764 write (iout,'(//a)') "Generalized inverse of fricmat"
765 call matout(dimen,dimen,maxres6,maxres6,fricmat)
770 if (nfgtasks.gt.1) then
771 if (fg_rank.eq.0) then
772 c The matching BROADCAST for fg processors is called in ERGASTULUM
778 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
780 time_Bcast=time_Bcast+MPI_Wtime()-time00
782 time_Bcast=time_Bcast+tcpu()-time00
784 c print *,"Processor",myrank,
785 c & " BROADCAST iorder in SETUP_FRICMAT"
788 c write (iout,*) "setup_fricmat licznik",licznik
794 c Scatter the friction matrix
795 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
796 & nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
797 & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
800 time_scatter=time_scatter+MPI_Wtime()-time00
801 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
803 time_scatter=time_scatter+tcpu()-time00
804 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
809 fricmat(j,i)=fcopy(i,j)
812 c write (iout,*) "My chunk of fricmat"
813 c call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
819 c-------------------------------------------------------------------------------
820 subroutine sdarea(gamvec)
822 c Scale the friction coefficients according to solvent accessible surface areas
823 c Code adapted from TINKER
828 include 'COMMON.CONTROL'
832 include 'COMMON.LANGEVIN'
835 include 'COMMON.LANGEVIN.lang0.5diag'
837 include 'COMMON.LANGEVIN.lang0'
840 include 'COMMON.CHAIN'
841 include 'COMMON.DERIV'
843 include 'COMMON.LOCAL'
844 include 'COMMON.INTERACT'
845 include 'COMMON.IOUNITS'
846 include 'COMMON.NAMES'
847 double precision radius(maxres2),gamvec(maxres2)
848 double precision twosix
849 parameter (twosix=1.122462048309372981d0)
850 logical lprn /.false./
852 double precision probe,area,ratio
854 c determine new friction coefficients every few SD steps
856 c set the atomic radii to estimates of sigma values
858 c print *,"Entered sdarea"
864 c Load peptide group radii
868 c Load side chain radii
871 if (iti.ne.ntyp1) radius(i+nres)=restok(iti)
874 c write (iout,*) "i",i," radius",radius(i)
877 radius(i) = radius(i) / twosix
878 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
881 c scale atomic friction coefficients by accessible area
883 if (lprn) write (iout,*)
884 & "Original gammas, surface areas, scaling factors, new gammas, ",
885 & "std's of stochastic forces"
888 if (radius(i).gt.0.0d0) then
889 call surfatom (i,area,radius)
890 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
891 if (lprn) write (iout,'(i5,3f10.5,$)')
892 & i,gamvec(ind+1),area,ratio
895 gamvec(ind) = ratio * gamvec(ind)
897 stdforcp(i)=stdfp*dsqrt(gamvec(ind))
898 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
902 if (radius(i+nres).gt.0.0d0) then
903 call surfatom (i+nres,area,radius)
904 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
905 if (lprn) write (iout,'(i5,3f10.5,$)')
906 & i,gamvec(ind+1),area,ratio
909 gamvec(ind) = ratio * gamvec(ind)
911 if (itype(i).ne.ntyp1)
912 & stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
913 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)