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
if (shield_mode.gt.0) then
call set_shield_fac
endif
include 'COMMON.FFIELD'
include 'COMMON.TIME1'
include 'COMMON.SPLITELE'
+ include 'COMMON.SHIELD'
dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
& erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
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)
+ else
ees=ees+eesij
+ endif
evdw1=evdw1+evdwij*sss
cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
C print *,'AFM',Eafmforce,totTafm*velAFMconst,dist
return
end
+C-----------------------------------------------------------
+C first for shielding is setting of function of side-chains
+ subroutine set_shield_fac
+ 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/
+
+C the vector between center of side_chain and peptide group
+ double precision pep_side(3),long,side_calf(3),
+ &pept_group(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
+ 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
+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)
+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)
+C now sscale fraction
+ sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
+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)=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.0*sh_frac_dist-3.0d0)
+ fac_help_scale=6.0*(scale_fac_dist-scale_fac_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)
+ 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+short**2/dist_pep_side**2)
+C now costhet_grad
+ costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side
+ 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.0
+ 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
+ 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)
+ VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi)
+ & /VSolvSphere_div
+ 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
+ fac_shield(i)=VolumeTotal*div77_81+div4_81
+ write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+ enddo
+ return
+ end
projects / unres.git / blobdiff