X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;fp=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=49676d8b767988768b851732bbc2c35e1c47c13d;hb=e910def1e1c1c0f0b35e8ce1cbf2a9ee59628771;hp=53d2eb45f6612702be760cf5553fbc3add875445;hpb=d101c97dea752458d76055fdbae49c26fff03c1f;p=unres.git diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 53d2eb4..49676d8 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -141,6 +141,7 @@ C Introduction of shielding effect first for each peptide group 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 @@ -3535,6 +3536,7 @@ C------------------------------------------------------------------------------- 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), @@ -3670,7 +3672,11 @@ C MARYSIA 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, @@ -10512,4 +10518,109 @@ C Eafmforce=-forceAFMconst*(dist-distafminit) 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