From 54453b9e3cc18886558071b1bce406b21e6566cf Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Wed, 2 Dec 2015 09:53:22 +0100 Subject: [PATCH] correction in shield --- source/unres/src_MD-M/COMMON.SHIELD | 4 +- source/unres/src_MD-M/checkder_p.F | 6 +- source/unres/src_MD-M/energy_p_new_barrier.F | 179 +++++++++++++++++++++++++- source/unres/src_MD-M/readrtns_CSA.F | 2 + 4 files changed, 184 insertions(+), 7 deletions(-) diff --git a/source/unres/src_MD-M/COMMON.SHIELD b/source/unres/src_MD-M/COMMON.SHIELD index 56a1c71..a8110d5 100644 --- a/source/unres/src_MD-M/COMMON.SHIELD +++ b/source/unres/src_MD-M/COMMON.SHIELD @@ -1,10 +1,10 @@ double precision VSolvSphere,VSolvSphere_div,long_r_sidechain, & short_r_sidechain,fac_shield,grad_shield_side,grad_shield, - & buff_shield + & buff_shield,wshield integer ishield_list,shield_list,ees0plist common /shield/ VSolvSphere,VSolvSphere_div,buff_shield, & long_r_sidechain(ntyp), - & short_r_sidechain(ntyp),fac_shield(maxres), + & short_r_sidechain(ntyp),fac_shield(maxres),wshield, & grad_shield_side(3,maxcont,-1:maxres),grad_shield(3,-1:maxres), & grad_shield_loc(3,maxcont,-1:maxres), & ishield_list(maxres),shield_list(maxcont,maxres), diff --git a/source/unres/src_MD-M/checkder_p.F b/source/unres/src_MD-M/checkder_p.F index 135aab7..270f4cc 100644 --- a/source/unres/src_MD-M/checkder_p.F +++ b/source/unres/src_MD-M/checkder_p.F @@ -14,8 +14,8 @@ C Check the gradient of Cartesian coordinates in internal coordinates. * Check the gradient of the virtual-bond and SC vectors in the internal * coordinates. * - aincr=1.0d-7 - aincr2=5.0d-8 + aincr=1.5d-7 + aincr2=2.5d-8 call cartder write (iout,'(a)') '**************** dx/dalpha' write (iout,'(a)') @@ -282,7 +282,7 @@ c rlambd=0.3d0 c call intcartderiv c call checkintcartgrad call zerograd - aincr=8.0D-6 + aincr=8.0D-7 write(iout,*) 'Calling CHECK_ECARTINT.' nf=0 icall=0 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 215fb6c..393ff9a 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -142,8 +142,10 @@ 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 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 @@ -11369,7 +11371,8 @@ C now costhet_grad enddo VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi) - & /VSolvSphere_div*4.0d0 + & /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, @@ -11448,4 +11451,176 @@ C print *, x(i+1),yy(i),i 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 + & *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))) + & )*div77_81 +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))) + & )*div77_81 + + 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) + & )) + & *div77_81 + enddo + VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist + enddo + fac_shield(i)=VolumeTotal*div77_81+div4_81 +C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) + enddo + return + end diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 79840a4..68fb297 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -600,6 +600,7 @@ C include 'COMMON.BOUNDS' include 'COMMON.MD' include 'COMMON.SETUP' + include 'COMMON.SHIELD' character*4 sequence(maxres) integer rescode double precision x(maxvar) @@ -641,6 +642,7 @@ C Read weights of the subsequent energy terms. call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) call reada(weightcard,'TEMP0',temp0,300.0d0) + call reada(weightcard,'WSHIELD',wshield,1.0d0) call reada(weightcard,'WLT',wliptran,0.0D0) if (index(weightcard,'SOFT').gt.0) ipot=6 C 12/1/95 Added weight for the multi-body term WCORR -- 1.7.9.5