C the shielding factor is set this factor is describing how each
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
- write (iout,*) "shield_mode",shield_mode
+C write (iout,*) "shield_mode",shield_mode
if (shield_mode.gt.0) then
call set_shield_fac
endif
& wstrain*ghpbc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
+ & +welec*gshieldc(j,i)
enddo
enddo
& wstrain*ghpbc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
+ & +welec*gshieldc(j,i)
enddo
enddo
do i=-1,nct
do j=1,3
#ifdef SPLITELE
+C print *,gradbufc(1,13)
+C print *,welec*gelc(1,13)
+C print *,wel_loc*gel_loc(1,13)
+C print *,0.5d0*(wscp*gvdwc_scpp(1,13))
+C print *,welec*gelc_long(1,13)+wvdwpp*gvdwpp(1,13)
+C print *,wel_loc*gel_loc_long(1,13)
+C print *,gradafm(1,13),"AFM"
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& wel_loc*gel_loc(j,i)+
& 0.5d0*(wscp*gvdwc_scpp(j,i)+
& +wscloc*gscloc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
+ & +welec*gshieldc(j,i)
+ & +welec*gshieldc_loc(j,i)
+
+
#else
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& wel_loc*gel_loc(j,i)+
& +wscloc*gscloc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
+ & +welec*gshieldc(j,i)
+ & +welec*gshieldc_loc(j,i)
+
#endif
gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
& wsccor*gsccorx(j,i)
& +wscloc*gsclocx(j,i)
& +wliptran*gliptranx(j,i)
+ & +welec*gshieldx(j,i)
enddo
enddo
#ifdef DEBUG
c
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
+CTU KURWA
do i=iatel_s,iatel_e
+C do i=75,75
if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
C changes suggested by Ana to avoid out of bounds
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
+C I TU KURWA
do j=ielstart(i),ielend(i)
+C do j=16,17
C write (iout,*) i,j
if (j.le.1) cycle
if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
el2=fac4*fac
C MARYSIA
- eesij=(el1+el2)
+C eesij=(el1+el2)
C 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
if (shield_mode.gt.0) then
- ees=ees+eesij*fac_shield(i)*fac_shield(j)
+C fac_shield(i)=0.4
+C fac_shield(j)=0.6
+ el1=el1*fac_shield(i)*fac_shield(j)
+ el2=el2*fac_shield(i)*fac_shield(j)
+ eesij=(el1+el2)
+ ees=ees+eesij
else
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+ eesij=(el1+el2)
ees=ees+eesij
endif
evdw1=evdw1+evdwij*sss
erij(1)=xj*rmij
erij(2)=yj*rmij
erij(3)=zj*rmij
+
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
ggg(1)=facel*xj
ggg(2)=facel*yj
ggg(3)=facel*zj
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (shield_mode.gt.0)) then
+C print *,i,j
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
+ gshieldx(k,iresshield)=gshieldx(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+ gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
+C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
+C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C if (iresshield.gt.i) then
+C do ishi=i+1,iresshield-1
+C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
+C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C enddo
+C else
+C do ishi=iresshield,i
+C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
+C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C enddo
+C endif
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
+ gshieldx(k,iresshield)=gshieldx(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+ gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
+
+C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
+C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+C if (iresshield.gt.j) then
+C do ishi=j+1,iresshield-1
+C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
+C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+C
+C enddo
+C else
+C do ishi=iresshield,j
+C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
+C & -grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+C enddo
+C endif
+ enddo
+ enddo
+
+ do k=1,3
+ gshieldc(k,i)=gshieldc(k,i)+
+ & grad_shield(k,i)*eesij/fac_shield(i)
+ gshieldc(k,j)=gshieldc(k,j)+
+ & grad_shield(k,j)*eesij/fac_shield(j)
+ gshieldc(k,i-1)=gshieldc(k,i-1)+
+ & grad_shield(k,i)*eesij/fac_shield(i)
+ gshieldc(k,j-1)=gshieldc(k,j-1)+
+ & grad_shield(k,j)*eesij/fac_shield(j)
+
+ enddo
+ endif
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gelc(k,i)=gelc(k,i)+ghalf
c gelc(k,j)=gelc(k,j)+ghalf
c enddo
c 9/28/08 AL Gradient compotents will be summed only at the end
+C print *,"before", gelc_long(1,i), gelc_long(1,j)
do k=1,3
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
+C & +grad_shield(k,j)*eesij/fac_shield(j)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
+C & +grad_shield(k,i)*eesij/fac_shield(i)
+C gelc_long(k,i-1)=gelc_long(k,i-1)
+C & +grad_shield(k,i)*eesij/fac_shield(i)
+C gelc_long(k,j-1)=gelc_long(k,j-1)
+C & +grad_shield(k,j)*eesij/fac_shield(j)
enddo
+C print *,"bafter", gelc_long(1,i), gelc_long(1,j)
+
*
* Loop over residues i+1 thru j-1.
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
ggg(1)=fac*xj
+C+eesij*grad_shield(1,i)+eesij*grad_shield(1,j)
ggg(2)=fac*yj
+C+eesij*grad_shield(2,i)+eesij*grad_shield(2,j)
ggg(3)=fac*zj
+C+eesij*grad_shield(3,i)+eesij*grad_shield(3,j)
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gelc(k,i)=gelc(k,i)+ghalf
cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
cd & (dcosg(k),k=1,3)
do k=1,3
- ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
+ & fac_shield(i)*fac_shield(j)
enddo
c do k=1,3
c ghalf=0.5D0*ggg(k)
cgrad gelc(l,k)=gelc(l,k)+ggg(l)
cgrad enddo
cgrad enddo
+C print *,"before22", gelc_long(1,i), gelc_long(1,j)
do k=1,3
gelc(k,i)=gelc(k,i)
- & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
- & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
+ & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
+ & *fac_shield(i)*fac_shield(j)
gelc(k,j)=gelc(k,j)
- & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
- & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
+ & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
+ & *fac_shield(i)*fac_shield(j)
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
+C print *,"before33", gelc_long(1,i), gelc_long(1,j)
+
C MARYSIA
c endif !sscale
IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
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/
+ &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)
+ &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
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(k+nres,j)-(c(i,j)+c(i+1,j))/2.0d0
+ 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(k+nres,j)-c(k,j)
- pept_group(j)=c(i,j)-c(i+1,j)
+ 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
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
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))=k
+ 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
else
scale_fac_dist=-sh_frac_dist*sh_frac_dist
& *(2.0*sh_frac_dist-3.0d0)
- fac_help_scale=6.0*(scale_fac_dist-scale_fac_dist**2)
+ fac_help_scale=6.0*(sh_frac_dist-sh_frac_dist**2)
& /dist_pep_side/buff_shield*0.5
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 print *,"jestem",scale_fac_dist,fac_help_scale,
+C & sh_frac_dist_grad(j)
enddo
endif
+C if ((i.eq.3).and.(k.eq.2)) then
+C print *,i,sh_frac_dist,dist_pep,fac_help_scale,scale_fac_dist
+C & ,"TU"
+C 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+short**2/dist_pep_side**2)
+ costhet=1.0d0/dsqrt(1.0+short**2/dist_pep_side**2)
C now costhet_grad
- costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**3
+C costhet=0.0d0
+ costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**4
+C costhet_fac=0.0d0
do j=1,3
costhet_grad(j)=costhet_fac*pep_side(j)
enddo
do j=1,3
pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
enddo
- fac_alfa_sin=1.0-(pep_side0pept_group/
- & (dist_pep_side*dist_side_calf))**2
+ cosalfa=(pep_side0pept_group/
+ & (dist_pep_side*dist_side_calf))
+ fac_alfa_sin=1.0-cosalfa**2
fac_alfa_sin=dsqrt(fac_alfa_sin)
rkprim=fac_alfa_sin*(long-short)+short
- cosphi=1.0d0/dsqrt(1+rkprim**2/dist_pep_side**2)
+C now costhet_grad
+ cosphi=1.0d0/dsqrt(1.0+rkprim**2/dist_pep_side**2)
+ cosphi_fac=cosphi**3*rkprim**2*(-0.5)/dist_pep_side**4
+
+ do j=1,3
+ cosphi_grad_long(j)=cosphi_fac*pep_side(j)
+ &+cosphi**3*0.5/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))
+
+ cosphi_grad_loc(j)=cosphi**3*0.5/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)
+ enddo
+
VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi)
& /VSolvSphere_div
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
+ & +(sh_frac_dist_grad(j)
C gradient po costhet
- &+scale_fac_dist*costhet_grad(j)
+ &-scale_fac_dist*costhet_grad(j)/(1.0-costhet)
+ &-scale_fac_dist*(cosphi_grad_long(j))
+ &/(1.0-cosphi) )*div77_81
+ &*VofOverlap
C grad_shield_side is Cbeta sidechain gradient
grad_shield_side(j,ishield_list(i),i)=
- & sh_frac_dist_grad(j)*VofOverlap*2.0d0
- & +scale_fac_dist*costhet_grad(j)*2.0d0
-C grad_shield_side_ca is Calfa sidechain gradient
- grad_shield_side_ca(j,ishield_list(i),i)=
+ & (sh_frac_dist_grad(j)*-2.0d0
+ & +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet)
+ & +scale_fac_dist*(cosphi_grad_long(j))
+ & *2.0d0/(1.0-cosphi))
+ & *div77_81*VofOverlap
+
+ grad_shield_loc(j,ishield_list(i),i)=
+ & scale_fac_dist*cosphi_grad_loc(j)
+ & *2.0d0/(1.0-cosphi)
+ & *div77_81*VofOverlap
enddo
VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
-C if ((cosphi.le.0.0).or.(costhet.le.0.0)) write(iout,*) "ERROR",
-C & cosphi,costhet
-C now should be fac_side_grad(k) which will be gradient of factor k which also
-C affect the gradient of peptide group i fac_pept_grad(i) and i+1
- write(2,*) "myvolume",VofOverlap,VSolvSphere_div,VolumeTotal
- enddo
-C write(2,*) "TOTAL VOLUME",i,VolumeTotal
-C the scaling factor of the shielding effect
+ enddo
fac_shield(i)=VolumeTotal*div77_81+div4_81
- write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
enddo
return
end
projects / unres.git / commitdiff