subroutine friction_force
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.VAR'
include 'COMMON.CHAIN'
include 'COMMON.LOCAL'
include 'COMMON.INTERACT'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.IOUNITS'
+#ifdef FIVEDIAG
+ integer iposc,ichain,n,innt,inct
+ double precision v_work(3,maxres2),vvec(maxres2_chain),rs(maxres2)
+#else
double precision gamvec(MAXRES6)
common /syfek/ gamvec
double precision vv(3),vvtot(3,maxres),v_work(MAXRES6),
& ginvfric(maxres2,maxres2)
common /przechowalnia/ ginvfric
+#endif
+ integer i,j,k,ind
logical lprn /.false./, checkmode /.false./
-
+#ifdef FIVEDIAG
+c Here accelerations due to friction forces are computed right after forces.
+ d_t_work=0.0d0
+ do j=1,3
+ v_work(j,1)=d_t(j,0)
+ v_work(j,nnt)=d_t(j,0)
+ enddo
+ do i=nnt+1,nct
+ do j=1,3
+ v_work(j,i)=v_work(j,i-1)+d_t(j,i-1)
+ enddo
+ enddo
+ do i=nnt,nct
+ if (iabs(itype(i)).ne.10 .and. iabs(itype(i)).ne.ntyp1) then
+ do j=1,3
+ v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres)
+ enddo
+ endif
+ enddo
+#ifdef DEBUG
+ write (iout,*) "v_work"
+ do i=1,2*nres
+ write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3)
+ enddo
+#endif
+ do j=1,3
+ ind=0
+ do ichain=1,nchain
+ n=dimen_chain(ichain)
+ iposc=iposd_chain(ichain)
+c write (iout,*) "friction_force j",j," ichain",ichain,
+c & " n",n," iposc",iposc,iposc+n-1
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ do i=innt,inct
+ vvec(ind+1)=v_work(j,i)
+ ind=ind+1
+ if (iabs(itype(i)).ne.10) then
+ vvec(ind+1)=v_work(j,i+nres)
+ ind=ind+1
+ endif
+ enddo
+#ifdef DEBUG
+ write (iout,*) "vvec ind",ind
+ write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
+#endif
+c write (iout,*) "chain",i," ind",ind," n",n
+ call fivediagmult(n,DMfric(iposc),DU1fric(iposc),
+ & DU2fric(iposc),vvec,rs)
+ do i=iposc,iposc+n-1
+ fric_work(3*(i-1)+j)=-rs(i)
+ enddo
+ enddo
+ enddo
+#ifdef DEBUG
+ write (iout,*) "Vector fric_work"
+ write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
+#endif
+#else
do i=0,MAXRES2
do j=1,3
friction(j,i)=0.0d0
enddo
enddo
endif
+#endif
return
end
c-----------------------------------------------------
subroutine stochastic_force(stochforcvec)
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
+ double precision time00
#endif
include 'COMMON.VAR'
include 'COMMON.CHAIN'
include 'COMMON.LOCAL'
include 'COMMON.INTERACT'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
include 'COMMON.TIME1'
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#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./
+ integer i,j,ind,ii,iti
+ double precision anorm_distr
+#ifdef FIVEDIAG
+ integer ichain,innt,inct,iposc
+#endif
+
do i=0,MAXRES2
do j=1,3
stochforc(j,i)=0.0d0
enddo
enddo
- x=0.0d0
+ x=0.0d0
#ifdef MPI
time00=MPI_Wtime()
force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
enddo
enddo
+#ifdef DEBUG
+ write (iout,*) "Stochastic forces on sites"
+ do i=1,nres
+ write (iout,'(i5,2(3f10.5,5x))') i,(force(j,i),j=1,3),
+ & (force(j,i+nres),j=1,3)
+ enddo
+#endif
#ifdef MPI
time_fsample=time_fsample+MPI_Wtime()-time00
#else
time_fsample=time_fsample+tcpu()-time00
#endif
+#ifdef FIVEDIAG
+ ind=0
+ do ichain=1,nchain
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ iposc=iposd_chain(ichain)
+c write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
+c & " inct",inct," iposc",iposc
+ do j=1,3
+ stochforcvec(ind+j)=0.5d0*force(j,innt)
+ enddo
+ if (iabs(itype(innt)).eq.10) then
+ do j=1,3
+ stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres)
+ enddo
+ ind=ind+3
+ else
+ ind=ind+3
+ do j=1,3
+ stochforcvec(ind+j)=force(j,innt+nres)
+ enddo
+ ind=ind+3
+ endif
+ do i=innt+1,inct-1
+ do j=1,3
+ stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1))
+ enddo
+ if (iabs(itype(i)).eq.10) then
+ do j=1,3
+ stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres)
+ enddo
+ ind=ind+3
+ else
+ ind=ind+3
+ do j=1,3
+ stochforcvec(ind+j)=force(j,i+nres)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ do j=1,3
+ stochforcvec(ind+j)=0.5d0*force(j,inct-1)
+ enddo
+ if (iabs(itype(inct)).eq.10) then
+ do j=1,3
+ stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres)
+ enddo
+ ind=ind+3
+ else
+ ind=ind+3
+ do j=1,3
+ stochforcvec(ind+j)=force(j,inct+nres)
+ enddo
+ ind=ind+3
+ endif
+c write (iout,*) "chain",ichain," ind",ind
+ enddo
+#ifdef DEBUG
+ write (iout,*) "stochforcvec"
+ write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind)
+#endif
+#else
c Compute the stochastic forces acting on virtual-bond vectors.
do j=1,3
ff(j)=0.0d0
enddo
endif
enddo
-
do j=1,3
stochforcvec(j)=stochforc(j,0)
enddo
ind=ind+3
endif
enddo
+#endif
if (lprn) then
write (iout,*) "stochforcvec"
do i=1,3*dimen
enddo
endif
-
return
end
c------------------------------------------------------------------
subroutine setup_fricmat
- implicit real*8 (a-h,o-z)
+ implicit none
#ifdef MPI
include 'mpif.h'
+ integer ierr
+ double precision time00
#endif
include 'DIMENSIONS'
include 'COMMON.VAR'
include 'COMMON.LOCAL'
include 'COMMON.INTERACT'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
include 'COMMON.SETUP'
include 'COMMON.TIME1'
c integer licznik /0/
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.IOUNITS'
integer IERROR
- integer i,j,ind,ind1,m
+ integer i,j,k,l,ind,ind1,m,ii,iti,it,nzero,innt,inct
+ integer ichain,nind
logical lprn /.false./
- double precision dtdi,gamvec(MAXRES2),
- & ginvfric(maxres2,maxres2),Ghalf(mmaxres2),fcopy(maxres2,maxres2)
+ double precision dtdi,gamvec(MAXRES2)
common /syfek/ gamvec
+#ifndef FIVEDIAG
+ double precision ginvfric(maxres2,maxres2),Ghalf(mmaxres2),
+ & fcopy(maxres2,maxres2)
double precision work(8*maxres2)
integer iwork(maxres2)
common /przechowalnia/ ginvfric,Ghalf,fcopy
+#endif
+
#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
enddo
c Load the friction coefficients corresponding to side chains
m=nct-nnt
+ if (lprn) write (iout,*) "m",m
ind=0
C gamsc(ntyp1)=1.0d0
do i=nnt,nct
gamvec(ii)=gamsc(iabs(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
-
+#ifdef FIVEDIAG
+ DMfric=0.0d0
+ DU1fric=0.0d0
+ DU2fric=0.0d0
+ ind=1
+ do ichain=1,nchain
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+c write (iout,*) "ichain",ichain," innt",innt," inct",inct
+c DMfric part
+ DMfric(ind)=gamvec(innt-nnt+1)/4
+ if (iabs(itype(innt)).eq.10) then
+ DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1)
+ ind=ind+1
+ else
+ DMfric(ind+1)=gamvec(m+innt-nnt+1)
+ ind=ind+2
+ endif
+c write (iout,*) "DMfric init ind",ind
+c DMfric
+ do i=innt+1,inct-1
+ DMfric(ind)=gamvec(i-nnt+1)/2
+ if (iabs(itype(i)).eq.10) then
+ DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1)
+ ind=ind+1
+ else
+ DMfric(ind+1)=gamvec(m+i-nnt+1)
+ ind=ind+2
+ endif
+ enddo
+c write (iout,*) "DMfric endloop ind",ind
+ if (inct.gt.innt) then
+ DMfric(ind)=gamvec(inct-1-nnt+1)/4
+ if (iabs(itype(inct)).eq.10) then
+ DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1)
+ ind=ind+1
+ else
+ DMfric(ind+1)=gamvec(inct+m-nnt+1)
+ ind=ind+2
+ endif
+ endif
+c write (iout,*) "DMfric end ind",ind
+ enddo
+c DU1fric part
+ do ichain=1,nchain
+ ind=iposd_chain(ichain)
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ do i=innt,inct
+ if (iabs(itype(i)).ne.10) then
+ ind=ind+2
+ else
+ DU1fric(ind)=gamvec(i-nnt+1)/4
+ ind=ind+1
+ endif
+ enddo
+ enddo
+c DU2fric part
+ do ichain=1,nchain
+ ind=iposd_chain(ichain)
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ do i=innt,inct-1
+ if (iabs(itype(i)).ne.10) then
+ DU2fric(ind)=gamvec(i-nnt+1)/4
+ DU2fric(ind+1)=0.0d0
+ ind=ind+2
+ else
+ DU2fric(ind)=0.0d0
+ ind=ind+1
+ endif
+ enddo
+ enddo
+ if (lprn) then
+ write(iout,*)"The upper part of the five-diagonal friction matrix"
+ do ichain=1,nchain
+ write (iout,'(a,i5)') 'Chain',ichain
+ innt=iposd_chain(ichain)
+ inct=iposd_chain(ichain)+dimen_chain(ichain)-1
+ do i=innt,inct
+ if (i.lt.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),
+ & DU2fric(i)
+ else if (i.eq.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i)
+ else
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i)
+ endif
+ enddo
+ enddo
+ endif
+ 10 continue
+#else
c The friction matrix
do k=1,dimen
do i=1,dimen
c call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
endif
#endif
+#endif
return
end
c-------------------------------------------------------------------------------
c Code adapted from TINKER
c AL 9/3/04
c
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
double precision radius(maxres2),gamvec(maxres2)
+ double precision twosix
parameter (twosix=1.122462048309372981d0)
logical lprn /.false./
+ integer i,j,iti,ind
+ double precision probe,area,ratio
c
c determine new friction coefficients every few SD steps
c
c Load side chain radii
do i=nnt,nct
iti=itype(i)
- radius(i+nres)=restok(iti)
+ if (iti.ne.ntyp1) radius(i+nres)=restok(iti)
enddo
c do i=1,2*nres
c write (iout,*) "i",i," radius",radius(i)
projects / unres.git / blobdiff