introduction of shielding effect (volume of overlap)
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Tue, 1 Sep 2015 14:29:59 +0000 (16:29 +0200)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Tue, 1 Sep 2015 14:29:59 +0000 (16:29 +0200)
CMakeCache.txt
source/unres/src_MD-M/COMMON.CONTROL
source/unres/src_MD-M/COMMON.SHIELD [new file with mode: 0644]
source/unres/src_MD-M/energy_p_new_barrier.F
source/unres/src_MD-M/parmread.F
source/unres/src_MD-M/readrtns_CSA.F
source/unres/src_MD-M/refsys.f

index a5f54f8..8bf73bf 100644 (file)
@@ -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
 
index 6c0980c..04cbfb7 100644 (file)
@@ -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 (file)
index 0000000..7974844
--- /dev/null
@@ -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)
+
+       
index 53d2eb4..49676d8 100644 (file)
@@ -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
 
index 1102c89..574552c 100644 (file)
@@ -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
index 91d7c07..4d70e89 100644 (file)
@@ -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
       
index 35a881e..4b7b763 100644 (file)
@@ -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.