--- /dev/null
+ subroutine friction_force
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.IOUNITS'
+ double precision gamvec(MAXRES6)
+ common /syfek/ gamvec
+ double precision vv(3),vvtot(3,maxres),v_work(MAXRES6),
+ & ginvfric(maxres2,maxres2)
+ common /przechowalnia/ ginvfric
+
+ logical lprn /.false./, checkmode /.false./
+
+ do i=0,MAXRES2
+ do j=1,3
+ friction(j,i)=0.0d0
+ enddo
+ enddo
+
+ do j=1,3
+ d_t_work(j)=d_t(j,0)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ d_t_work(ind+j)=d_t(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+ d_t_work(ind+j)=d_t(j,i+nres)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+
+ call fricmat_mult(d_t_work,fric_work)
+
+ if (.not.checkmode) return
+
+ if (lprn) then
+ write (iout,*) "d_t_work and fric_work"
+ do i=1,3*dimen
+ write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
+ enddo
+ endif
+ do j=1,3
+ friction(j,0)=fric_work(j)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ friction(j,i)=fric_work(ind+j)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+ friction(j,i+nres)=fric_work(ind+j)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ if (lprn) then
+ write(iout,*) "Friction backbone"
+ do i=0,nct-1
+ write(iout,'(i5,3e15.5,5x,3e15.5)')
+ & i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
+ enddo
+ write(iout,*) "Friction side chain"
+ do i=nnt,nct
+ write(iout,'(i5,3e15.5,5x,3e15.5)')
+ & i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
+ enddo
+ endif
+ if (lprn) then
+ do j=1,3
+ vv(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct
+ do j=1,3
+ vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
+ vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
+ vv(j)=vv(j)+d_t(j,i)
+ enddo
+ enddo
+ write (iout,*) "vvtot backbone and sidechain"
+ do i=nnt,nct
+ write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),
+ & (vvtot(j,i+nres),j=1,3)
+ enddo
+ ind=0
+ do i=nnt,nct-1
+ do j=1,3
+ v_work(ind+j)=vvtot(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ do j=1,3
+ v_work(ind+j)=vvtot(j,i+nres)
+ enddo
+ ind=ind+3
+ enddo
+ write (iout,*) "v_work gamvec and site-based friction forces"
+ do i=1,dimen1
+ write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),
+ & gamvec(i)*v_work(i)
+ enddo
+c do i=1,dimen
+c fric_work1(i)=0.0d0
+c do j=1,dimen1
+c fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
+c enddo
+c enddo
+c write (iout,*) "fric_work and fric_work1"
+c do i=1,dimen
+c write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
+c enddo
+ do i=1,dimen
+ do j=1,dimen
+ ginvfric(i,j)=0.0d0
+ do k=1,dimen
+ ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
+ enddo
+ enddo
+ enddo
+ write (iout,*) "ginvfric"
+ do i=1,dimen
+ write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
+ enddo
+ write (iout,*) "symmetry check"
+ do i=1,dimen
+ do j=1,i-1
+ write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
+ enddo
+ enddo
+ endif
+ return
+ end
+c-----------------------------------------------------
+ subroutine stochastic_force(stochforcvec)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.MD'
+ include 'COMMON.TIME1'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.IOUNITS'
+
+ double precision x,sig,lowb,highb,
+ & ff(3),force(3,0:MAXRES2),zeta2,lowb2,
+ & highb2,sig2,forcvec(MAXRES6),stochforcvec(MAXRES6)
+ logical lprn /.false./
+ do i=0,MAXRES2
+ do j=1,3
+ stochforc(j,i)=0.0d0
+ enddo
+ enddo
+ x=0.0d0
+
+#ifdef MPI
+ time00=MPI_Wtime()
+#else
+ time00=tcpu()
+#endif
+c Compute the stochastic forces acting on bodies. Store in force.
+ do i=nnt,nct-1
+ sig=stdforcp(i)
+ lowb=-5*sig
+ highb=5*sig
+ do j=1,3
+ force(j,i)=anorm_distr(x,sig,lowb,highb)
+ enddo
+ enddo
+ do i=nnt,nct
+ sig2=stdforcsc(i)
+ lowb2=-5*sig2
+ highb2=5*sig2
+ do j=1,3
+ force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
+ enddo
+ enddo
+#ifdef MPI
+ time_fsample=time_fsample+MPI_Wtime()-time00
+#else
+ time_fsample=time_fsample+tcpu()-time00
+#endif
+c Compute the stochastic forces acting on virtual-bond vectors.
+ do j=1,3
+ ff(j)=0.0d0
+ enddo
+ do i=nct-1,nnt,-1
+ do j=1,3
+ stochforc(j,i)=ff(j)+0.5d0*force(j,i)
+ enddo
+ do j=1,3
+ ff(j)=ff(j)+force(j,i)
+ enddo
+ if (itype(i+1).ne.21) then
+ do j=1,3
+ stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
+ ff(j)=ff(j)+force(j,i+nres+1)
+ enddo
+ endif
+ enddo
+ do j=1,3
+ stochforc(j,0)=ff(j)+force(j,nnt+nres)
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+ stochforc(j,i+nres)=force(j,i+nres)
+ enddo
+ endif
+ enddo
+
+ do j=1,3
+ stochforcvec(j)=stochforc(j,0)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ stochforcvec(ind+j)=stochforc(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+ stochforcvec(ind+j)=stochforc(j,i+nres)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ if (lprn) then
+ write (iout,*) "stochforcvec"
+ do i=1,3*dimen
+ write(iout,'(i5,e15.5)') i,stochforcvec(i)
+ enddo
+ write(iout,*) "Stochastic forces backbone"
+ do i=0,nct-1
+ write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
+ enddo
+ write(iout,*) "Stochastic forces side chain"
+ do i=nnt,nct
+ write(iout,'(i5,3e15.5)')
+ & i,(stochforc(j,i+nres),j=1,3)
+ enddo
+ endif
+
+ if (lprn) then
+
+ ind=0
+ do i=nnt,nct-1
+ write (iout,*) i,ind
+ do j=1,3
+ forcvec(ind+j)=force(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ write (iout,*) i,ind
+ do j=1,3
+ forcvec(j+ind)=force(j,i+nres)
+ enddo
+ ind=ind+3
+ enddo
+
+ write (iout,*) "forcvec"
+ ind=0
+ do i=nnt,nct-1
+ do j=1,3
+ write (iout,'(2i3,2f10.5)') i,j,force(j,i),
+ & forcvec(ind+j)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ do j=1,3
+ write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),
+ & forcvec(ind+j)
+ enddo
+ ind=ind+3
+ enddo
+
+ endif
+
+ return
+ end
+c------------------------------------------------------------------
+ subroutine setup_fricmat
+ implicit real*8 (a-h,o-z)
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.MD'
+ include 'COMMON.SETUP'
+ include 'COMMON.TIME1'
+c integer licznik /0/
+c save licznik
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.IOUNITS'
+ integer IERROR
+ integer i,j,ind,ind1,m
+ logical lprn /.false./
+ double precision dtdi,gamvec(MAXRES2),
+ & ginvfric(maxres2,maxres2),Ghalf(mmaxres2),fcopy(maxres2,maxres2)
+ common /syfek/ gamvec
+ double precision work(8*maxres2)
+ integer iwork(maxres2)
+ common /przechowalnia/ ginvfric,Ghalf,fcopy
+#ifdef MPI
+ if (fg_rank.ne.king) goto 10
+#endif
+c Zeroing out fricmat
+ do i=1,dimen
+ do j=1,dimen
+ fricmat(i,j)=0.0d0
+ enddo
+ enddo
+c Load the friction coefficients corresponding to peptide groups
+ ind1=0
+ do i=nnt,nct-1
+ ind1=ind1+1
+ gamvec(ind1)=gamp
+ enddo
+c Load the friction coefficients corresponding to side chains
+ m=nct-nnt
+ ind=0
+ do i=nnt,nct
+ ind=ind+1
+ ii = ind+m
+ iti=itype(i)
+ gamvec(ii)=gamsc(iti)
+ enddo
+ if (surfarea) call sdarea(gamvec)
+c if (lprn) then
+c write (iout,*) "Matrix A and vector gamma"
+c do i=1,dimen1
+c write (iout,'(i2,$)') i
+c do j=1,dimen
+c write (iout,'(f4.1,$)') A(i,j)
+c enddo
+c write (iout,'(f8.3)') gamvec(i)
+c enddo
+c endif
+ if (lprn) then
+ write (iout,*) "Vector gamvec"
+ do i=1,dimen1
+ write (iout,'(i5,f10.5)') i, gamvec(i)
+ enddo
+ endif
+
+c The friction matrix
+ do k=1,dimen
+ do i=1,dimen
+ dtdi=0.0d0
+ do j=1,dimen1
+ dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
+ enddo
+ fricmat(k,i)=dtdi
+ enddo
+ enddo
+
+ if (lprn) then
+ write (iout,'(//a)') "Matrix fricmat"
+ call matout2(dimen,dimen,maxres2,maxres2,fricmat)
+ endif
+ if (lang.eq.2 .or. lang.eq.3) then
+c Mass-scale the friction matrix if non-direct integration will be performed
+ do i=1,dimen
+ do j=1,dimen
+ Ginvfric(i,j)=0.0d0
+ do k=1,dimen
+ do l=1,dimen
+ Ginvfric(i,j)=Ginvfric(i,j)+
+ & Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
+ enddo
+ enddo
+ enddo
+ enddo
+c Diagonalize the friction matrix
+ ind=0
+ do i=1,dimen
+ do j=1,i
+ ind=ind+1
+ Ghalf(ind)=Ginvfric(i,j)
+ enddo
+ enddo
+ call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
+ & ierr,iwork)
+ if (lprn) then
+ write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
+ & " mass-scaled friction matrix"
+ call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
+ endif
+c Precompute matrices for tinker stochastic integrator
+#ifndef LANG0
+ do i=1,dimen
+ do j=1,dimen
+ mt1(i,j)=0.0d0
+ mt2(i,j)=0.0d0
+ do k=1,dimen
+ mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
+ mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
+ enddo
+ mt3(j,i)=mt1(i,j)
+ enddo
+ enddo
+#endif
+ else if (lang.eq.4) then
+c Diagonalize the friction matrix
+ ind=0
+ do i=1,dimen
+ do j=1,i
+ ind=ind+1
+ Ghalf(ind)=fricmat(i,j)
+ enddo
+ enddo
+ call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
+ & ierr,iwork)
+ if (lprn) then
+ write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
+ & " friction matrix"
+ call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
+ endif
+c Determine the number of zero eigenvalues of the friction matrix
+ nzero=max0(dimen-dimen1,0)
+c do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
+c nzero=nzero+1
+c enddo
+ write (iout,*) "Number of zero eigenvalues:",nzero
+ do i=1,dimen
+ do j=1,dimen
+ fricmat(i,j)=0.0d0
+ do k=nzero+1,dimen
+ fricmat(i,j)=fricmat(i,j)
+ & +fricvec(i,k)*fricvec(j,k)/fricgam(k)
+ enddo
+ enddo
+ enddo
+ if (lprn) then
+ write (iout,'(//a)') "Generalized inverse of fricmat"
+ call matout(dimen,dimen,maxres6,maxres6,fricmat)
+ endif
+ endif
+#ifdef MPI
+ 10 continue
+ if (nfgtasks.gt.1) then
+ if (fg_rank.eq.0) then
+c The matching BROADCAST for fg processors is called in ERGASTULUM
+#ifdef MPI
+ time00=MPI_Wtime()
+#else
+ time00=tcpu()
+#endif
+ call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
+#ifdef MPI
+ time_Bcast=time_Bcast+MPI_Wtime()-time00
+#else
+ time_Bcast=time_Bcast+tcpu()-time00
+#endif
+c print *,"Processor",myrank,
+c & " BROADCAST iorder in SETUP_FRICMAT"
+ endif
+c licznik=licznik+1
+c write (iout,*) "setup_fricmat licznik",licznik
+#ifdef MPI
+ time00=MPI_Wtime()
+#else
+ time00=tcpu()
+#endif
+c Scatter the friction matrix
+ call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
+ & nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
+ & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
+#ifdef TIMING
+#ifdef MPI
+ time_scatter=time_scatter+MPI_Wtime()-time00
+ time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
+#else
+ time_scatter=time_scatter+tcpu()-time00
+ time_scatter_fmat=time_scatter_fmat+tcpu()-time00
+#endif
+#endif
+ do i=1,dimen
+ do j=1,2*my_ng_count
+ fricmat(j,i)=fcopy(i,j)
+ enddo
+ enddo
+c write (iout,*) "My chunk of fricmat"
+c call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
+ endif
+#endif
+ return
+ end
+c-------------------------------------------------------------------------------
+ subroutine sdarea(gamvec)
+c
+c Scale the friction coefficients according to solvent accessible surface areas
+c Code adapted from TINKER
+c AL 9/3/04
+c
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ double precision radius(maxres2),gamvec(maxres2)
+ parameter (twosix=1.122462048309372981d0)
+ logical lprn /.false./
+c
+c determine new friction coefficients every few SD steps
+c
+c set the atomic radii to estimates of sigma values
+c
+c print *,"Entered sdarea"
+ probe = 0.0d0
+
+ do i=1,2*nres
+ radius(i)=0.0d0
+ enddo
+c Load peptide group radii
+ do i=nnt,nct-1
+ radius(i)=pstok
+ enddo
+c Load side chain radii
+ do i=nnt,nct
+ iti=itype(i)
+ radius(i+nres)=restok(iti)
+ enddo
+c do i=1,2*nres
+c write (iout,*) "i",i," radius",radius(i)
+c enddo
+ do i = 1, 2*nres
+ radius(i) = radius(i) / twosix
+ if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
+ end do
+c
+c scale atomic friction coefficients by accessible area
+c
+ if (lprn) write (iout,*)
+ & "Original gammas, surface areas, scaling factors, new gammas, ",
+ & "std's of stochastic forces"
+ ind=0
+ do i=nnt,nct-1
+ if (radius(i).gt.0.0d0) then
+ call surfatom (i,area,radius)
+ ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
+ if (lprn) write (iout,'(i5,3f10.5,$)')
+ & i,gamvec(ind+1),area,ratio
+ do j=1,3
+ ind=ind+1
+ gamvec(ind) = ratio * gamvec(ind)
+ enddo
+ stdforcp(i)=stdfp*dsqrt(gamvec(ind))
+ if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
+ endif
+ enddo
+ do i=nnt,nct
+ if (radius(i+nres).gt.0.0d0) then
+ call surfatom (i+nres,area,radius)
+ ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
+ if (lprn) write (iout,'(i5,3f10.5,$)')
+ & i,gamvec(ind+1),area,ratio
+ do j=1,3
+ ind=ind+1
+ gamvec(ind) = ratio * gamvec(ind)
+ enddo
+ stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
+ if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
+ endif
+ enddo
+
+ return
+ end
projects / unres.git / blobdiff