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),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
501 parameter (mmaxres2=(maxres2*(maxres2+1)/2))
504 include 'COMMON.CHAIN'
505 include 'COMMON.DERIV'
507 include 'COMMON.LOCAL'
508 include 'COMMON.INTERACT'
511 include 'COMMON.LAGRANGE.5diag'
513 include 'COMMON.LAGRANGE'
515 include 'COMMON.SETUP'
516 include 'COMMON.TIME1'
517 c integer licznik /0/
520 include 'COMMON.LANGEVIN'
523 include 'COMMON.LANGEVIN.lang0.5diag'
525 include 'COMMON.LANGEVIN.lang0'
528 include 'COMMON.IOUNITS'
530 integer i,j,k,l,ind,ind1,m,ii,iti,it,nzero,innt,inct
532 logical lprn /.false./
533 double precision dtdi,gamvec(MAXRES2)
534 common /syfek/ gamvec
536 double precision ginvfric(maxres2,maxres2),Ghalf(mmaxres2),
537 & fcopy(maxres2,maxres2)
538 double precision work(8*maxres2)
539 integer iwork(maxres2)
540 common /przechowalnia/ ginvfric,Ghalf,fcopy
544 if (fg_rank.ne.king) goto 10
551 c Load the friction coefficients corresponding to side chains
553 if (lprn) write (iout,*) "m",m
560 gamvec(ii)=gamsc(iabs(iti))
562 if (surfarea) call sdarea(gamvec)
564 write (iout,*) "Vector gamvec ii",ii
566 write (iout,'(i5,f10.5)') i, gamvec(i)
575 innt=chain_border(1,ichain)
576 inct=chain_border(2,ichain)
577 c write (iout,*) "ichain",ichain," innt",innt," inct",inct
579 DMfric(ind)=gamvec(innt-nnt+1)/4
580 if (iabs(itype(innt)).eq.10) then
581 DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1)
584 DMfric(ind+1)=gamvec(m+innt-nnt+1)
587 c write (iout,*) "DMfric init ind",ind
590 DMfric(ind)=gamvec(i-nnt+1)/2
591 if (iabs(itype(i)).eq.10) then
592 DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1)
595 DMfric(ind+1)=gamvec(m+i-nnt+1)
599 c write (iout,*) "DMfric endloop ind",ind
600 if (inct.gt.innt) then
601 DMfric(ind)=gamvec(inct-1-nnt+1)/4
602 if (iabs(itype(inct)).eq.10) then
603 DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1)
606 DMfric(ind+1)=gamvec(inct+m-nnt+1)
610 c write (iout,*) "DMfric end ind",ind
614 ind=iposd_chain(ichain)
615 innt=chain_border(1,ichain)
616 inct=chain_border(2,ichain)
618 if (iabs(itype(i)).ne.10) then
621 DU1fric(ind)=gamvec(i-nnt+1)/4
628 ind=iposd_chain(ichain)
629 innt=chain_border(1,ichain)
630 inct=chain_border(2,ichain)
632 if (iabs(itype(i)).ne.10) then
633 DU2fric(ind)=gamvec(i-nnt+1)/4
643 write(iout,*)"The upper part of the five-diagonal friction matrix"
645 write (iout,'(a,i5)') 'Chain',ichain
646 innt=iposd_chain(ichain)
647 inct=iposd_chain(ichain)+dimen_chain(ichain)-1
649 if (i.lt.inct-1) then
650 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),
652 else if (i.eq.inct-1) then
653 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i)
655 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i)
662 c The friction matrix
667 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
674 write (iout,'(//a)') "Matrix fricmat"
675 call matout2(dimen,dimen,maxres2,maxres2,fricmat)
677 if (lang.eq.2 .or. lang.eq.3) then
678 c Mass-scale the friction matrix if non-direct integration will be performed
684 Ginvfric(i,j)=Ginvfric(i,j)+
685 & Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
690 c Diagonalize the friction matrix
695 Ghalf(ind)=Ginvfric(i,j)
698 call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
701 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
702 & " mass-scaled friction matrix"
703 call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
705 c Precompute matrices for tinker stochastic integrator
712 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
713 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
719 else if (lang.eq.4) then
720 c Diagonalize the friction matrix
725 Ghalf(ind)=fricmat(i,j)
728 call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
731 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
733 call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
735 c Determine the number of zero eigenvalues of the friction matrix
736 nzero=max0(dimen-dimen1,0)
737 c do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
740 write (iout,*) "Number of zero eigenvalues:",nzero
745 fricmat(i,j)=fricmat(i,j)
746 & +fricvec(i,k)*fricvec(j,k)/fricgam(k)
751 write (iout,'(//a)') "Generalized inverse of fricmat"
752 call matout(dimen,dimen,maxres6,maxres6,fricmat)
757 if (nfgtasks.gt.1) then
758 if (fg_rank.eq.0) then
759 c The matching BROADCAST for fg processors is called in ERGASTULUM
765 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
767 time_Bcast=time_Bcast+MPI_Wtime()-time00
769 time_Bcast=time_Bcast+tcpu()-time00
771 c print *,"Processor",myrank,
772 c & " BROADCAST iorder in SETUP_FRICMAT"
775 c write (iout,*) "setup_fricmat licznik",licznik
781 c Scatter the friction matrix
782 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
783 & nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
784 & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
787 time_scatter=time_scatter+MPI_Wtime()-time00
788 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
790 time_scatter=time_scatter+tcpu()-time00
791 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
796 fricmat(j,i)=fcopy(i,j)
799 c write (iout,*) "My chunk of fricmat"
800 c call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
806 c-------------------------------------------------------------------------------
807 subroutine sdarea(gamvec)
809 c Scale the friction coefficients according to solvent accessible surface areas
810 c Code adapted from TINKER
815 include 'COMMON.CONTROL'
819 include 'COMMON.LANGEVIN'
822 include 'COMMON.LANGEVIN.lang0.5diag'
824 include 'COMMON.LANGEVIN.lang0'
827 include 'COMMON.CHAIN'
828 include 'COMMON.DERIV'
830 include 'COMMON.LOCAL'
831 include 'COMMON.INTERACT'
832 include 'COMMON.IOUNITS'
833 include 'COMMON.NAMES'
834 double precision radius(maxres2),gamvec(maxres2)
835 double precision twosix
836 parameter (twosix=1.122462048309372981d0)
837 logical lprn /.false./
839 double precision probe,area,ratio
841 c determine new friction coefficients every few SD steps
843 c set the atomic radii to estimates of sigma values
845 c print *,"Entered sdarea"
851 c Load peptide group radii
855 c Load side chain radii
858 if (iti.ne.ntyp1) radius(i+nres)=restok(iti)
861 c write (iout,*) "i",i," radius",radius(i)
864 radius(i) = radius(i) / twosix
865 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
868 c scale atomic friction coefficients by accessible area
870 if (lprn) write (iout,*)
871 & "Original gammas, surface areas, scaling factors, new gammas, ",
872 & "std's of stochastic forces"
875 if (radius(i).gt.0.0d0) then
876 call surfatom (i,area,radius)
877 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
878 if (lprn) write (iout,'(i5,3f10.5,$)')
879 & i,gamvec(ind+1),area,ratio
882 gamvec(ind) = ratio * gamvec(ind)
884 stdforcp(i)=stdfp*dsqrt(gamvec(ind))
885 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
889 if (radius(i+nres).gt.0.0d0) then
890 call surfatom (i+nres,area,radius)
891 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
892 if (lprn) write (iout,'(i5,3f10.5,$)')
893 & i,gamvec(ind+1),area,ratio
896 gamvec(ind) = ratio * gamvec(ind)
898 if (itype(i).ne.ntyp1)
899 & stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
900 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)