X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Fstochfric.F;h=dc0b088a7cc66e98582d64a8cdd336051b8637f9;hb=a30bd29e64da2aa47b84963fdd0bf4192ead2738;hp=f38dfda2a71465827e092543c26b83d3cc0a4e41;hpb=020e579626d686ec20ecd9f0cc4c8313f474e152;p=unres.git diff --git a/source/unres/src-HCD-5D/stochfric.F b/source/unres/src-HCD-5D/stochfric.F index f38dfda..dc0b088 100644 --- a/source/unres/src-HCD-5D/stochfric.F +++ b/source/unres/src-HCD-5D/stochfric.F @@ -1,5 +1,5 @@ subroutine friction_force - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.VAR' include 'COMMON.CHAIN' @@ -8,20 +8,93 @@ 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 @@ -150,14 +223,16 @@ c enddo 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' @@ -166,24 +241,39 @@ c----------------------------------------------------- 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() @@ -207,11 +297,80 @@ c Compute the stochastic forces acting on bodies. Store in force. 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 @@ -240,7 +399,6 @@ c Compute the stochastic forces acting on virtual-bond vectors. enddo endif enddo - do j=1,3 stochforcvec(j)=stochforc(j,0) enddo @@ -259,6 +417,7 @@ c Compute the stochastic forces acting on virtual-bond vectors. ind=ind+3 endif enddo +#endif if (lprn) then write (iout,*) "stochforcvec" do i=1,3*dimen @@ -311,14 +470,15 @@ c Compute the stochastic forces acting on virtual-bond vectors. 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' @@ -328,6 +488,11 @@ c------------------------------------------------------------------ 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/ @@ -335,28 +500,30 @@ c save licznik #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 @@ -364,6 +531,7 @@ c Load the friction coefficients corresponding to peptide groups 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 @@ -373,23 +541,105 @@ C gamsc(ntyp1)=1.0d0 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 @@ -531,6 +781,7 @@ c write (iout,*) "My chunk of fricmat" c call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy) endif #endif +#endif return end c------------------------------------------------------------------------------- @@ -540,7 +791,7 @@ 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) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' @@ -548,8 +799,12 @@ c #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' @@ -558,8 +813,11 @@ c 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 @@ -578,7 +836,7 @@ c Load peptide group radii 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)