changes in wham

[unres.git] / source / wham / src-M / parmread.F

diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F
index d4fa521..097e1f7 100644 (file)
--- a/source/wham/src-M/parmread.F
+++ b/source/wham/src-M/parmread.F
@@ -21,6 +21,8 @@ C
       include 'COMMON.SCCOR'
       include 'COMMON.SCROT'
       include 'COMMON.FREE'
+      include 'COMMON.SHIELD'
+      include 'COMMON.CONTROL'
       character*1 t1,t2,t3
       character*1 onelett(4) /"G","A","P","D"/
       character*1 toronelet(-2:2) /"p","a","G","A","P"/
@@ -141,6 +143,9 @@ c
       wbond=ww(18)
       wsccor=ww(19)
       whpb=ww(15)
+      wstrain=ww(15)
+      wliptran=ww(22)
+      wshield=ww(25)
       endif
 
       call card_concat(controlcard,.false.)
@@ -212,7 +217,7 @@ c Read the virtual-bond parameters, masses, and moments of inertia
 c and Stokes' radii of the peptide group and side chains
 c
 #ifdef CRYST_BOND
-      read (ibond,*) vbldp0,akp
+      read (ibond,*) vbldp0,vbldpdum,akp
       do i=1,ntyp
         nbondterm(i)=1
         read (ibond,*) vbldsc0(1,i),aksc(1,i)
@@ -224,7 +229,7 @@ c
         endif
       enddo
 #else
-      read (ibond,*) ijunk,vbldp0,akp,rjunk
+      read (ibond,*) ijunk,vbldp0,vbldpdum,akp,rjunk
       do i=1,ntyp
         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
      &   j=1,nbondterm(i))
@@ -250,6 +255,11 @@ c
           enddo
         enddo
       endif
+       read(iliptranpar,*) pepliptran
+       do i=1,ntyp
+       read(iliptranpar,*) liptranene(i)
+       enddo
+       close(iliptranpar)
 #ifdef CRYST_THETA
 C
 C Read the parameters of the probability distribution/energy expression 
@@ -394,7 +404,7 @@ C
         enddo
       enddo
       enddo
-      write (iout,*) "KURWA1"
+C      write (iout,*) "KURWA1"
       do iblock=1,2
       do i=0,nthetyp
         do j=-nthetyp,nthetyp
@@ -418,7 +428,7 @@ C
           enddo
         enddo
       enddo
-       write(iout,*) "KURWA1.1"
+C       write(iout,*) "KURWA1.1"
 C
 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 C coefficients of theta-and-gamma-dependent terms are zero.
@@ -438,7 +448,7 @@ C
         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
       enddo
       enddo
-       write(iout,*) "KURWA1.5"
+C       write(iout,*) "KURWA1.5"
 C Substitution for D aminoacids from symmetry.
       do iblock=1,2
       do i=-nthetyp,0
@@ -517,7 +527,7 @@ C
       call flush(iout)
       endif
 #endif
-      write(iout,*) 'KURWA2'
+C      write(iout,*) 'KURWA2'
 #ifdef CRYST_SC
 C
 C Read the parameters of the probability distribution/energy expression
@@ -625,7 +635,7 @@ C
       enddo
 #endif
       close(irotam)
-      write (iout,*) 'KURWAKURWA'
+C      write (iout,*) 'KURWAKURWA'
 #ifdef CRYST_TOR
 C
 C Read torsional parameters in old format
@@ -1100,6 +1110,13 @@ C---------------------- GB or BP potential -----------------------------
       read (isidep,*)(sigii(i),i=1,ntyp)
       read (isidep,*)(chip(i),i=1,ntyp)
       read (isidep,*)(alp(i),i=1,ntyp)
+      do i=1,ntyp
+       read (isidep,*)(epslip(i,j),j=i,ntyp)
+C       print *,"WARNING!!"
+C       do j=1,ntyp
+C       epslip(i,j)=epslip(i,j)+0.05d0
+C       enddo
+      enddo
 C For the GB potential convert sigma'**2 into chi'
       if (ipot.eq.4) then
        do i=1,ntyp
@@ -1138,6 +1155,7 @@ C Calculate the "working" parameters of SC interactions.
       do i=2,ntyp
         do j=1,i-1
          eps(i,j)=eps(j,i)
+          epslip(i,j)=epslip(j,i)
         enddo
       enddo
       do i=1,ntyp
@@ -1155,6 +1173,7 @@ C Calculate the "working" parameters of SC interactions.
       do i=1,ntyp
        do j=i,ntyp
          epsij=eps(i,j)
+          epsijlip=epslip(i,j)
          if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
            rrij=sigma(i,j)
           else
@@ -1166,10 +1185,16 @@ C Calculate the "working" parameters of SC interactions.
          epsij=eps(i,j)
          sigeps=dsign(1.0D0,epsij)
          epsij=dabs(epsij)
-         aa(i,j)=epsij*rrij*rrij
-         bb(i,j)=-sigeps*epsij*rrij
-         aa(j,i)=aa(i,j)
-         bb(j,i)=bb(i,j)
+         aa_aq(i,j)=epsij*rrij*rrij
+         bb_aq(i,j)=-sigeps*epsij*rrij
+         aa_aq(j,i)=aa_aq(i,j)
+         bb_aq(j,i)=bb_aq(i,j)
+          sigeps=dsign(1.0D0,epsijlip)
+          epsijlip=dabs(epsijlip)
+          aa_lip(i,j)=epsijlip*rrij*rrij
+          bb_lip(i,j)=-sigeps*epsijlip*rrij
+          aa_lip(j,i)=aa_lip(i,j)
+          bb_lip(j,i)=bb_lip(i,j)
          if (ipot.gt.2) then
            sigt1sq=sigma0(i)**2
            sigt2sq=sigma0(j)**2
@@ -1279,7 +1304,7 @@ C      V3SS = 13.7d0
       write (iout,*) dyn_ss,'dyndyn'
       if (dyn_ss) then
         ss_depth=ebr/wsc-0.25*eps(1,1)
-        write(iout,*) akcm,whpb,wsc,'KURWA'
+C        write(iout,*) akcm,whpb,wsc,'KURWA'
         Ht=Ht/wsc-0.25*eps(1,1)
 
         akcm=akcm*whpb/wsc
@@ -1300,5 +1325,22 @@ C      if (lprint) then
       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
      & ' v3ss:',v3ss
 C      endif
+      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
+      buff_shield=1.0d0
+      endif 
       return
       end