From e910def1e1c1c0f0b35e8ce1cbf2a9ee59628771 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Tue, 1 Sep 2015 16:29:59 +0200 Subject: [PATCH] introduction of shielding effect (volume of overlap) --- CMakeCache.txt | 16 ++-- source/unres/src_MD-M/COMMON.CONTROL | 5 +- source/unres/src_MD-M/COMMON.SHIELD | 9 +++ source/unres/src_MD-M/energy_p_new_barrier.F | 111 ++++++++++++++++++++++++++ source/unres/src_MD-M/parmread.F | 55 ++++++++++--- source/unres/src_MD-M/readrtns_CSA.F | 29 ++++++- source/unres/src_MD-M/refsys.f | 3 +- 7 files changed, 206 insertions(+), 22 deletions(-) create mode 100644 source/unres/src_MD-M/COMMON.SHIELD diff --git a/CMakeCache.txt b/CMakeCache.txt index a5f54f8..8bf73bf 100644 --- a/CMakeCache.txt +++ b/CMakeCache.txt @@ -1,13 +1,13 @@ # This is the CMakeCache file. # For build in directory: /users2/adasko/unres -# It was generated by CMake: /users2/adasko/apps/cmake/bin/cmake +# It was generated by CMake: /usr/bin/cmake # You can edit this file to change values found and used by cmake. # If you do not want to change any of the values, simply exit the editor. # If you do want to change a value, simply edit, save, and exit the editor. # The syntax for the file is as follows: # KEY:TYPE=VALUE # KEY is the name of a variable in the cache. -# TYPE is a hint to GUIs for the type of VALUE, DO NOT EDIT TYPE!. +# TYPE is a hint to GUI's for the type of VALUE, DO NOT EDIT TYPE!. # VALUE is the current value for the KEY. ######################## @@ -26,15 +26,15 @@ CMAKE_CACHE_MAJOR_VERSION:INTERNAL=2 //Minor version of cmake used to create the current loaded cache CMAKE_CACHE_MINOR_VERSION:INTERNAL=8 //Patch version of cmake used to create the current loaded cache -CMAKE_CACHE_PATCH_VERSION:INTERNAL=12 +CMAKE_CACHE_PATCH_VERSION:INTERNAL=7 //Path to CMake executable. -CMAKE_COMMAND:INTERNAL=/users2/adasko/apps/cmake/bin/cmake +CMAKE_COMMAND:INTERNAL=/usr/bin/cmake //Path to cpack program executable. -CMAKE_CPACK_COMMAND:INTERNAL=/users2/adasko/apps/cmake/bin/cpack +CMAKE_CPACK_COMMAND:INTERNAL=/usr/bin/cpack //Path to ctest program executable. -CMAKE_CTEST_COMMAND:INTERNAL=/users2/adasko/apps/cmake/bin/ctest +CMAKE_CTEST_COMMAND:INTERNAL=/usr/bin/ctest //Path to cache edit program executable. -CMAKE_EDIT_COMMAND:INTERNAL=/users2/adasko/apps/cmake/bin/ccmake +CMAKE_EDIT_COMMAND:INTERNAL=/usr/bin/ccmake //Path to CMake installation. -CMAKE_ROOT:INTERNAL=/users2/adasko/apps/cmake/share/cmake-2.8 +CMAKE_ROOT:INTERNAL=/usr/share/cmake-2.8 diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL index 6c0980c..04cbfb7 100644 --- a/source/unres/src_MD-M/COMMON.CONTROL +++ b/source/unres/src_MD-M/COMMON.CONTROL @@ -1,5 +1,6 @@ integer modecalc,iscode,indpdb,indback,indphi,iranconf,icheckgrad, - & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide + & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide, + & shield_mode logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec, & sideadd,lsecondary,read_cart,unres_pdb, & vdisulf,searchsc,lmuca,dccart,extconf,out1file, @@ -8,7 +9,7 @@ & icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,iprint, & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file, - & selfguide,AFMlog, + & selfguide,AFMlog,shield_mode, & constr_dist,gnorm_check,gradout,split_ene,with_theta_constr, & symetr C... minim = .true. means DO minimization. diff --git a/source/unres/src_MD-M/COMMON.SHIELD b/source/unres/src_MD-M/COMMON.SHIELD new file mode 100644 index 0000000..7974844 --- /dev/null +++ b/source/unres/src_MD-M/COMMON.SHIELD @@ -0,0 +1,9 @@ + double precision VSolvSphere,VSolvSphere_div,long_r_sidechain, + & short_r_sidechain,fac_shield + integer ishield_list,shield_list + common /shield/ VSolvSphere,VSolvSphere_div, + & long_r_sidechain(ntyp), + & short_r_sidechain(ntyp),fac_shield(maxres), + & ishield_list(maxres),shield_list(maxres) + + 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 diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 1102c89..574552c 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -26,6 +26,8 @@ C include 'COMMON.SBRIDGE' include 'COMMON.MD' include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.SHIELD' character*1 t1,t2,t3 character*1 onelett(4) /"G","A","P","D"/ character*1 toronelet(-2:2) /"p","a","G","A","P"/ @@ -945,16 +947,16 @@ c B2(1,i) = b(2) c B2(2,i) = b(4) c B2(1,-i) =b(2) c B2(2,-i) =-b(4) - B1tilde(1,i) = b(3) - B1tilde(2,i) =-b(5) - B1tilde(1,-i) =-b(3) - B1tilde(2,-i) =b(5) + B1tilde(1,i) = b(3,i) + B1tilde(2,i) =-b(5,i) + B1tilde(1,-i) =-b(3,i) + B1tilde(2,-i) =b(5,i) b1tilde(1,i)=0.0d0 b1tilde(2,i)=0.0d0 - B2(1,i) = b(2) - B2(2,i) = b(4) - B2(1,-i) =b(2) - B2(2,-i) =-b(4) + B2(1,i) = b(2,i) + B2(2,i) = b(4,i) + B2(1,-i) =b(2,i) + B2(2,-i) =-b(4,i) c b2(1,i)=0.0d0 c b2(2,i)=0.0d0 @@ -1333,6 +1335,25 @@ C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, C & ' v3ss:',v3ss C endif +C set the variables used for shielding effect +C write (iout,*) "SHIELD MODE",shield_mode +C if (shield_mode.gt.0) then +C VSolvSphere the volume of solving sphere +C print *,pi,"pi" +C rpp(1,1) is the energy r0 for peptide group contact and will be used for it +C there will be no distinction between proline peptide group and normal peptide +C group in case of shielding parameters +C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 +C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 +C write (iout,*) VSolvSphere,VSolvSphere_div +C long axis of side chain +C do i=1,ntyp +C long_r_sidechain(i)=vbldsc0(1,i) +C short_r_sidechain(i)=sigma0(i) +C enddo +C lets set the buffor value +C buff_shield=1.0d0 +C endif return 111 write (iout,*) "Error reading bending energy parameters." goto 999 @@ -1404,6 +1425,22 @@ c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--' #else call getenv(var,val) #endif - +C set the variables used for shielding effect +C if (shield_mode.gt.0) then +C VSolvSphere the volume of solving sphere +C print *,pi,"pi" +C rpp(1,1) is the energy r0 for peptide group contact and will be used for it +C there will be no distinction between proline peptide group and normal peptide +C group in case of shielding parameters +C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 +C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 +C long axis of side chain +C do i=1,ntyp +C long_r_sidechain(i)=vbldsc0(1,i) +C short_r_sidechain(i)=sigma0(i) +C enddo +C lets set the buffor value +C buff_shield=1.0d0 +C endif return end diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 91d7c07..4d70e89 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -81,6 +81,7 @@ C include 'COMMON.INTERACT' include 'COMMON.SETUP' include 'COMMON.SPLITELE' + include 'COMMON.SHIELD' COMMON /MACHSW/ KDIAG,ICORFL,IXDR character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/ character*80 ucase @@ -145,6 +146,15 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword selfguide=(index(controlcard,'SELFGUIDE')) print *,'AFMlog',AFMlog,selfguide,"KUPA" call readi(controlcard,'IPRINT',iprint,0) +C SHIELD keyword sets if the shielding effect of side-chains is used +C 0 denots no shielding is used all peptide are equally despite the +C solvent accesible area +C 1 the newly introduced function +C 2 reseved for further possible developement + call readi(controlcard,'SHIELD',shield_mode,0) +C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then + write(iout,*) "shield_mode",shield_mode +C endif call readi(controlcard,'MAXGEN',maxgen,10000) call readi(controlcard,'MAXOVERLAP',maxoverlap,1000) call readi(controlcard,"KDIAG",kdiag,0) @@ -248,8 +258,23 @@ C endif write(iout,*) "bordlipbot=",bordlipbot write(iout,*) "bufliptop=",bufliptop write(iout,*) "buflipbot=",buflipbot - - + write (iout,*) "SHIELD MODE",shield_mode + if (shield_mode.gt.0) then + pi=3.141592d0 +C VSolvSphere the volume of solving sphere +C print *,pi,"pi" +C rpp(1,1) is the energy r0 for peptide group contact and will be used for it +C there will be no distinction between proline peptide group and normal peptide +C group in case of shielding parameters + VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 + VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 + write (iout,*) VSolvSphere,VSolvSphere_div +C long axis of side chain + do i=1,ntyp + long_r_sidechain(i)=vbldsc0(1,i) + short_r_sidechain(i)=sigma0(i) + enddo + endif if (me.eq.king .or. .not.out1file ) & write (iout,*) "DISTCHAINMAX",distchainmax diff --git a/source/unres/src_MD-M/refsys.f b/source/unres/src_MD-M/refsys.f index 35a881e..4b7b763 100644 --- a/source/unres/src_MD-M/refsys.f +++ b/source/unres/src_MD-M/refsys.f @@ -3,6 +3,8 @@ c This subroutine calculates unit vectors of a local reference system c defined by atoms (i2), (i3), and (i4). The x axis is the axis from implicit real*8 (a-h,o-z) include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include "COMMON.CHAIN" c this subroutine calculates unity vectors of a local reference system c defined by atoms (i2), (i3), and (i4). the x axis is the axis from c atom (i3) to atom (i2), and the xy plane is the plane defined by atoms @@ -13,7 +15,6 @@ c form a linear fragment. returns vectors e1, e2, and e3. logical fail double precision e1(3),e2(3),e3(3) double precision u(3),z(3) - include 'COMMON.IOUNITS' double precision coinc/1.0D-13/,align /1.0D-13/ c print *,'just initialize' fail=.false. -- 1.7.9.5