C peptide group is shielded by side-chains
C the matrix - shield_fac(i) the i index describe the ith between i and i+1
C write (iout,*) "shield_mode",shield_mode
- if (shield_mode.gt.0) then
+ if (shield_mode.eq.1) then
call set_shield_fac
+ else if (shield_mode.eq.2) then
+ call set_shield_fac2
endif
c print *,"Processor",myrank," left VEC_AND_DERIV"
if (ipot.lt.6) then
C Calculate the virtual-bond torsional energy.
C
cd print *,'nterm=',nterm
+C print *,"tor",tor_mode
if (wtor.gt.0) then
if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then
call etor(etors,edihcnstr)
C etor kcc is Kubo cumulant clustered rigorous attemp to derive the
C energy function
if ((tor_mode.eq.1).or.(tor_mode.eq.2)) then
- call etor(etors,edihcnstr)
+ call etor_kcc(etors,edihcnstr)
endif
else
etors=0
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
include 'COMMON.SBRIDGE'
+ include 'COMMON.CONTROL'
double precision kfac /2.4d0/
double precision x,x2,x3,x4,x5,licznik /1.12692801104297249644/
c facT=temp0/t_bath
#endif
stop 555
endif
+ if (shield_mode.gt.0) then
+ wscp=weights(2)*fact
+ wsc=weights(1)*fact
+ wvdwpp=weights(16)*fact
+ endif
welec=weights(3)*fact
wcorr=weights(4)*fact3
wcorr5=weights(5)*fact4
C Set lprn=.true. for debugging
lprn=.false.
c lprn=.true.
+C print *,"wchodze kcc"
+ if (tor_mode.ne.2) then
etors=0.0D0
+ endif
do i=iphi_start,iphi_end
C ANY TWO ARE DUMMY ATOMS in row CYCLE
c if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
cosphi=dcos(j*phii)
sinphi=dsin(j*phii)
sint1t2n=sint1t2**j
- sumnonchebyshev=sumnonchebyshev+
+ sumnonchebyshev=
& sint1t2n*(v1ij*cosphi+v2ij*sinphi)
actval=sint1t2n*(v1ij*cosphi+v2ij*sinphi)
C etors=etors+sint1t2n*(v1ij*cosphi+v2ij*sinphi)
C if (energy_dec) etors_ii=etors_ii+
C & v1ij*cosphi+v2ij*sinphi
C glocig is the gradient local i site in gamma
- glocig=glocig+j*(v2ij*cosphi-v1ij*sinphi)*sint1t2n
+ glocig=j*(v2ij*cosphi-v1ij*sinphi)*sint1t2n
C now gradient over theta_1
- glocit1=glocit1+actval/sinthet1*j*costhet1
- glocit2=glocit2+actval/sinthet2*j*costhet2
- enddo
+ glocit1=actval/sinthet1*j*costhet1
+ glocit2=actval/sinthet2*j*costhet2
C now the Czebyshev polinominal sum
- do j=1,nterm_kcc_Tb(itori,itori1)
- thybt1(j)=v1_chyb(j,itori,itori1)
- thybt2(j)=v2_chyb(j,itori,itori1)
+ do k=1,nterm_kcc_Tb(itori,itori1)
+ thybt1(k)=v1_chyb(k,j,itori,itori1)
+ thybt2(k)=v2_chyb(k,j,itori,itori1)
+C thybt1(k)=0.0
+C thybt2(k)=0.0
enddo
sumth1thyb=tschebyshev
- & (1,nterm_kcc_Tb(itori,itori1),thybt1(1),dcos(theta12))
+ & (1,nterm_kcc_Tb(itori,itori1),thybt1(1),dcos(theti12)**2)
gradthybt1=gradtschebyshev
- & (0,nterm_kcc_Tb(itori,itori1)-1,thybt1(1),dcos(theta12))
- & *(nterm_kcc_Tb(itori,itori1))*0.5*dsin(theta12)
+ & (0,nterm_kcc_Tb(itori,itori1)-1,thybt1(1),
+ & dcos(theti12)**2)
+ & *dcos(theti12)*(-dsin(theti12))
sumth2thyb=tschebyshev
- & (1,nterm_kcc_Tb(itori,itori1),thybt2(1),dcos(theta22))
+ & (1,nterm_kcc_Tb(itori,itori1),thybt2(1),dcos(theti22)**2)
gradthybt2=gradtschebyshev
- & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1),dcos(theta22))
- & *(nterm_kcc_Tb(itori,itori1))*0.5*dsin(theta22)
+ & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1),
+ & dcos(theti22)**2)
+ & *dcos(theti22)*(-dsin(theti22))
+C print *, sumth1thyb, gradthybt1, sumth2thyb, gradthybt2,
+C & gradtschebyshev
+C & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1),
+C & dcos(theti22)**2),
+C & dsin(theti22)
C now overal sumation
etors=etors+(1.0d0+sumth1thyb+sumth2thyb)*sumnonchebyshev
+C print *,"sumnon", sumnonchebyshev,sumth1thyb+sumth2thyb
C derivative over gamma
gloc(i-3,icg)=gloc(i-3,icg)+wtor*glocig
& *(1.0d0+sumth1thyb+sumth2thyb)
C derivative over theta1
- gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wang*
+ gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wtor*
& (glocit1*(1.0d0+sumth1thyb+sumth2thyb)+
& sumnonchebyshev*gradthybt1)
C now derivative over theta2
- gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wtor*
& (glocit2*(1.0d0+sumth1thyb+sumth2thyb)+
& sumnonchebyshev*gradthybt2)
enddo
-
+ enddo
C gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
! 6/20/98 - dihedral angle constraints
return
end
+C The rigorous attempt to derive energy function
+ subroutine ebend_kcc(etheta,ethetacnstr)
-
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.TORSION'
+ include 'COMMON.INTERACT'
+ include 'COMMON.DERIV'
+ include 'COMMON.CHAIN'
+ include 'COMMON.NAMES'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.FFIELD'
+ include 'COMMON.TORCNSTR'
+ include 'COMMON.CONTROL'
+ logical lprn
+ double precision thybt1(maxtermkcc)
+C Set lprn=.true. for debugging
+ lprn=.false.
+c lprn=.true.
+C print *,"wchodze kcc"
+ if (tormode.ne.2) etheta=0.0D0
+ do i=ithet_start,ithet_end
+c print *,i,itype(i-1),itype(i),itype(i-2)
+ if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
+ & .or.itype(i).eq.ntyp1) cycle
+ iti=itortyp_kcc(itype(i-1))
+ sinthet=dsin(theta(i)/2.0d0)
+ costhet=dcos(theta(i)/2.0d0)
+ do j=1,nbend_kcc_Tb(iti)
+ thybt1(j)=v1bend_chyb(j,iti)
+ enddo
+ sumth1thyb=tschebyshev
+ & (1,nbend_kcc_Tb(iti),thybt1(1),costhet)
+ ihelp=nbend_kcc_Tb(iti)-1
+ gradthybt1=gradtschebyshev
+ & (0,ihelp,thybt1(1),costhet)
+ etheta=etheta+sumth1thyb
+C print *,sumth1thyb,gradthybt1,sinthet*(-0.5d0)
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*
+ & gradthybt1*sinthet*(-0.5d0)
+ enddo
+ if (tormode.ne.2) then
+ ethetacnstr=0.0d0
+C print *,ithetaconstr_start,ithetaconstr_end,"TU"
+ do i=ithetaconstr_start,ithetaconstr_end
+ itheta=itheta_constr(i)
+ thetiii=theta(itheta)
+ difi=pinorm(thetiii-theta_constr0(i))
+ if (difi.gt.theta_drange(i)) then
+ difi=difi-theta_drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else if (difi.lt.-drange(i)) then
+ difi=difi+drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else
+ difi=0.0
+ endif
+ if (energy_dec) then
+ write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
+ & i,itheta,rad2deg*thetiii,
+ & rad2deg*theta_constr0(i), rad2deg*theta_drange(i),
+ & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
+ & gloc(itheta+nphi-2,icg)
+ endif
+ enddo
+ endif
+ return
+ end
c------------------------------------------------------------------------------
subroutine eback_sc_corr(esccor)
c 7/21/2007 Correlations between the backbone-local and side-chain-local
VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi)
& /VSolvSphere_div
+ & *wshield
C now the gradient...
C grad_shield is gradient of Calfa for peptide groups
C write(iout,*) "shield_compon",i,k,VSolvSphere,scale_fac_dist,
implicit none
include "DIMENSIONS"
integer i,m,n
- double precision x(n),y,yy(0:maxvar),aux
+ double precision x(n+1),y,yy(0:maxvar),aux
c Tschebyshev polynomial. Note that the first term is omitted
c m=0: the constant term is included
c m=1: the constant term is not included
yy(0)=1.0d0
yy(1)=2.0d0*y
do i=2,n
- yy(i)=2*yy(1)*yy(i-1)-yy(i-2)
+ yy(i)=2*y*yy(i-1)-yy(i-2)
enddo
aux=0.0d0
do i=m,n
- aux=aux+x(i)*yy(i)
+ aux=aux+x(i+1)*yy(i)*(i+1)
+C print *, x(i+1),yy(i),i
enddo
gradtschebyshev=aux
return
end
+C------------------------------------------------------------------------
+C first for shielding is setting of function of side-chains
+ subroutine set_shield_fac2
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SHIELD'
+ include 'COMMON.INTERACT'
+C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
+ double precision div77_81/0.974996043d0/,
+ &div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
+
+C the vector between center of side_chain and peptide group
+ double precision pep_side(3),long,side_calf(3),
+ &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
+ &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
+C the line belowe needs to be changed for FGPROC>1
+ do i=1,nres-1
+ if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+ ishield_list(i)=0
+Cif there two consequtive dummy atoms there is no peptide group between them
+C the line below has to be changed for FGPROC>1
+ VolumeTotal=0.0
+ do k=1,nres
+ if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
+ dist_pep_side=0.0
+ dist_side_calf=0.0
+ do j=1,3
+C first lets set vector conecting the ithe side-chain with kth side-chain
+ pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0
+C pep_side(j)=2.0d0
+C and vector conecting the side-chain with its proper calfa
+ side_calf(j)=c(j,k+nres)-c(j,k)
+C side_calf(j)=2.0d0
+ pept_group(j)=c(j,i)-c(j,i+1)
+C lets have their lenght
+ dist_pep_side=pep_side(j)**2+dist_pep_side
+ dist_side_calf=dist_side_calf+side_calf(j)**2
+ dist_pept_group=dist_pept_group+pept_group(j)**2
+ enddo
+ dist_pep_side=dsqrt(dist_pep_side)
+ dist_pept_group=dsqrt(dist_pept_group)
+ dist_side_calf=dsqrt(dist_side_calf)
+ do j=1,3
+ pep_side_norm(j)=pep_side(j)/dist_pep_side
+ side_calf_norm(j)=dist_side_calf
+ enddo
+C now sscale fraction
+ sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
+C print *,buff_shield,"buff"
+C now sscale
+ if (sh_frac_dist.le.0.0) cycle
+C If we reach here it means that this side chain reaches the shielding sphere
+C Lets add him to the list for gradient
+ ishield_list(i)=ishield_list(i)+1
+C ishield_list is a list of non 0 side-chain that contribute to factor gradient
+C this list is essential otherwise problem would be O3
+ shield_list(ishield_list(i),i)=k
+C Lets have the sscale value
+ if (sh_frac_dist.gt.1.0) then
+ scale_fac_dist=1.0d0
+ do j=1,3
+ sh_frac_dist_grad(j)=0.0d0
+ enddo
+ else
+ scale_fac_dist=-sh_frac_dist*sh_frac_dist
+ & *(2.0d0*sh_frac_dist-3.0d0)
+ fac_help_scale=6.0d0*(sh_frac_dist-sh_frac_dist**2)
+ & /dist_pep_side/buff_shield*0.5d0
+C remember for the final gradient multiply sh_frac_dist_grad(j)
+C for side_chain by factor -2 !
+ do j=1,3
+ sh_frac_dist_grad(j)=fac_help_scale*pep_side(j)
+C sh_frac_dist_grad(j)=0.0d0
+C scale_fac_dist=1.0d0
+C print *,"jestem",scale_fac_dist,fac_help_scale,
+C & sh_frac_dist_grad(j)
+ enddo
+ endif
+C this is what is now we have the distance scaling now volume...
+ short=short_r_sidechain(itype(k))
+ long=long_r_sidechain(itype(k))
+ costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2)
+ sinthet=short/dist_pep_side*costhet
+C now costhet_grad
+C costhet=0.6d0
+C sinthet=0.8
+ costhet_fac=costhet**3*short**2*(-0.5d0)/dist_pep_side**4
+C sinthet_fac=costhet**2*0.5d0*(short**3/dist_pep_side**4*costhet
+C & -short/dist_pep_side**2/costhet)
+C costhet_fac=0.0d0
+ do j=1,3
+ costhet_grad(j)=costhet_fac*pep_side(j)
+ enddo
+C remember for the final gradient multiply costhet_grad(j)
+C for side_chain by factor -2 !
+C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1
+C pep_side0pept_group is vector multiplication
+ pep_side0pept_group=0.0d0
+ do j=1,3
+ pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
+ enddo
+ cosalfa=(pep_side0pept_group/
+ & (dist_pep_side*dist_side_calf))
+ fac_alfa_sin=1.0d0-cosalfa**2
+ fac_alfa_sin=dsqrt(fac_alfa_sin)
+ rkprim=fac_alfa_sin*(long-short)+short
+C rkprim=short
+
+C now costhet_grad
+ cosphi=1.0d0/dsqrt(1.0d0+rkprim**2/dist_pep_side**2)
+C cosphi=0.6
+ cosphi_fac=cosphi**3*rkprim**2*(-0.5d0)/dist_pep_side**4
+ sinphi=rkprim/dist_pep_side/dsqrt(1.0d0+rkprim**2/
+ & dist_pep_side**2)
+C sinphi=0.8
+ do j=1,3
+ cosphi_grad_long(j)=cosphi_fac*pep_side(j)
+ &+cosphi**3*0.5d0/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa/
+ &((dist_pep_side*dist_side_calf))*
+ &((side_calf(j))-cosalfa*
+ &((pep_side(j)/dist_pep_side)*dist_side_calf))
+C cosphi_grad_long(j)=0.0d0
+ cosphi_grad_loc(j)=cosphi**3*0.5d0/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa
+ &/((dist_pep_side*dist_side_calf))*
+ &(pep_side(j)-
+ &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
+C cosphi_grad_loc(j)=0.0d0
+ enddo
+C print *,sinphi,sinthet
+ VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet))
+ & /VSolvSphere_div
+C & *wshield
+C now the gradient...
+ do j=1,3
+ grad_shield(j,i)=grad_shield(j,i)
+C gradient po skalowaniu
+ & +(sh_frac_dist_grad(j)*VofOverlap
+C gradient po costhet
+ & +scale_fac_dist*VSolvSphere/VSolvSphere_div/4.0d0*
+ &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(
+ & sinphi/sinthet*costhet*costhet_grad(j)
+ & +sinthet/sinphi*cosphi*cosphi_grad_long(j)))
+ & )*wshield
+C grad_shield_side is Cbeta sidechain gradient
+ grad_shield_side(j,ishield_list(i),i)=
+ & (sh_frac_dist_grad(j)*-2.0d0
+ & *VofOverlap
+ & -scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*
+ &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(
+ & sinphi/sinthet*costhet*costhet_grad(j)
+ & +sinthet/sinphi*cosphi*cosphi_grad_long(j)))
+ & )*wshield
+
+ grad_shield_loc(j,ishield_list(i),i)=
+ & scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*
+ &(1.0d0/(dsqrt(1.0d0-sinphi*sinthet))*(
+ & sinthet/sinphi*cosphi*cosphi_grad_loc(j)
+ & ))
+ & *wshield
+ enddo
+ VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
+ enddo
+ fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
+C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+ enddo
+ return
+ end
projects / unres.git / blobdiff