introduction of shielding effect (volume of overlap)

[unres.git] / source / unres / src_MD-M / readrtns_CSA.F

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index 91d7c07..4d70e89 100644 (file)
--- 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