+++ /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(maxres6)
- 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