correction in wham and UNRES for lipid and correlation
[unres.git] / source / unres / src_MD-M / parmread.F
index 154fcf4..2d2c23b 100644 (file)
@@ -26,12 +26,14 @@ 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"/
       logical lprint,LaTeX
       dimension blower(3,3,maxlob)
-      dimension b(13)
+C      dimension b(13)
       character*3 lancuch,ucase
 C
 C For printing parameters after they are read set the following in the UNRES
@@ -59,7 +61,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,mp,ip,pstok
+      read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
       do i=1,ntyp
         nbondterm(i)=1
         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
@@ -71,8 +73,9 @@ c
         endif
       enddo
 #else
-      read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
+      read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
       do i=1,ntyp
+      print *,i
         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
      &   j=1,nbondterm(i)),msc(i),isc(i),restok(i)
         dsc(i) = vbldsc0(1,i)
@@ -97,6 +100,14 @@ c
           enddo
         enddo
       endif
+C reading lipid parameters
+      write (iout,*) "iliptranpar",iliptranpar
+      call flush(iout)
+       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 
@@ -210,6 +221,8 @@ C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 C
       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
      &  ntheterm3,nsingle,ndouble
+      write (iout,*) "ithep",ithep
+      call flush(iout)
       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
       do i=-ntyp1,-1
@@ -352,11 +365,12 @@ C Control printout of the coefficients of virtual-bond-angle potentials
 C
       if (lprint) then
         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
-        do i=1,nthetyp+1
-          do j=1,nthetyp+1
-            do k=1,nthetyp+1
+        do iblock=1,2
+        do i=0,nthetyp
+          do j=-nthetyp,nthetyp
+            do k=-nthetyp,nthetyp
               write (iout,'(//4a)') 
-     &         'Type ',onelett(i),onelett(j),onelett(k) 
+     &         'Type ',toronelet(i),toronelet(j),toronelet(k) 
               write (iout,'(//a,10x,a)') " l","a[l]"
               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
               write (iout,'(i2,1pe15.5)')
@@ -385,16 +399,19 @@ C
             enddo
           enddo
         enddo
+        enddo
       enddo
       call flush(iout)
       endif
-      write (2,*) "Start reading THETA_PDB"
+      write (2,*) "Start reading THETA_PDB",ithep_pdb
       do i=1,ntyp
-        read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
+c      write (2,*) 'i=',i
+        read (ithep_pdb,*,err=111,end=111)
+     &     a0thet(i),(athet(j,i,1,1),j=1,2),
      &    (bthet(j,i,1,1),j=1,2)
-        read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
-        read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
-        read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
+        read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
+        read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
+        read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
         sigc0(i)=sigc0(i)**2
       enddo
       do i=1,ntyp
@@ -539,6 +556,7 @@ C
 C Read the parameters of the probability distribution/energy expression
 C of the side chains.
 C
+      write (2,*) "Start reading ROTAM_PDB"
       do i=1,ntyp
         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
         if (i.eq.10) then
@@ -577,6 +595,7 @@ C
         endif
       enddo
       close (irotam_pdb)
+      write (2,*) "End reading ROTAM_PDB"
 #endif
       close(irotam)
 
@@ -648,8 +667,9 @@ c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
       close (itorp)
       if (lprint) then
         write (iout,'(/a/)') 'Torsional constants:'
-        do i=1,ntortyp
-          do j=1,ntortyp
+        do iblock=1,2
+        do i=0,ntortyp-1
+          do j=-ntortyp+1,ntortyp-1
             write (iout,*) 'ityp',i,' jtyp',j
             write (iout,*) 'Fourier constants'
             do k=1,nterm(i,j,iblock)
@@ -663,6 +683,7 @@ c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
             enddo
           enddo
         enddo
+        enddo
       endif
 
 C
@@ -762,6 +783,80 @@ C Martix of D parameters for two dimesional fourier series
       enddo
       endif
 #endif
+C read Czybyshev torsional parameters
+      read (itorkcc,*,end=121,err=121) nkcctyp
+      read (itorkcc,*,end=121,err=121) (itortyp_kcc(i),i=1,ntyp)
+      do i=-ntyp,-1
+        itortyp_kcc(i)=-itortyp_kcc(-i)
+      enddo
+      do i=0,nkcctyp
+        do j=0,nkcctyp
+C first we read the cos and sin gamma parameters
+          read (itorkcc,*,end=121,err=121) 
+     &    nterm_kcc(j,i),nterm_kcc_Tb(j,i)
+C           read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
+          do k=1,nterm_kcc(j,i)
+            do l=1,nterm_kcc_Tb(j,i)
+              read (itorkcc,*,end=121,err=121) v11_chyb(l,k,j,i)
+            enddo
+            do l=1,nterm_kcc_Tb(j,i)
+              read (itorkcc,*,end=121,err=121) v21_chyb(l,k,j,i)
+            enddo
+            do l=1,nterm_kcc_Tb(j,i)
+              read (itorkcc,*,end=121,err=121) v12_chyb(l,k,j,i)
+            enddo
+            do l=1,nterm_kcc_Tb(j,i)
+              read (itorkcc,*,end=121,err=121) v22_chyb(l,k,j,i)
+            enddo
+            read (itorkcc,*,end=121,err=121) v1_kcc(k,j,i)
+            read (itorkcc,*,end=121,err=121) v2_kcc(k,j,i)
+          enddo
+        enddo
+      enddo
+      if (lprint) then
+c Print valence-torsional parameters
+        write (iout,'(a)') 
+     &    "Parameters of the valence-torsional potentials"
+        do i=0,nkcctyp
+        do j=0,nkcctyp
+        write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
+        write (iout,'(2a20,a15)') "v_kcc","v1_chyb","v2_chyb"
+        do k=1,nterm_kcc(j,i)
+          write (iout,'(i5,f15.10,i5,2f15.10)') 
+     &      k,v1_kcc(k,j,i),1,v11_chyb(1,k,j,i),v21_chyb(1,k,j,i)
+          do l=2,nterm_kcc_Tb(j,i)
+            write (iout,'(20x,i5,2f15.10)') 
+     &        l,v11_chyb(l,k,j,i),v21_chyb(l,k,j,i)
+          enddo
+          write (iout,'(i5,f15.10,i5,2f15.10)') 
+     &      k,v2_kcc(k,j,i),1,v12_chyb(1,k,j,i),v22_chyb(1,k,j,i)
+          do l=2,nterm_kcc_Tb(j,i)
+            write (iout,'(20x,i5,2f15.10)') 
+     &        l,v12_chyb(l,k,j,i),v22_chyb(l,k,j,i)
+          enddo
+          write (iout,'(a)')
+        enddo
+        enddo
+        enddo
+      endif
+C here will be the apropriate recalibrating for D-aminoacid
+C        read (ithetkcc,*,end=121,err=121) nkcctyp
+      do i=0,nkcctyp
+        read (ithetkcc,*,end=121,err=121) nbend_kcc_Tb(i)
+        do j=1,nbend_kcc_Tb(i)
+          read (ithetkcc,*,end=121,err=121) v1bend_chyb(j,i)
+        enddo
+      enddo
+      if (lprint) then
+        write (iout,'(a)') 
+     &    "Parameters of the valence-only potentials"
+        do i=0,nkcctyp
+        write (iout,'(2a)') "Type ",toronelet(i)
+        do k=1,nbend_kcc_Tb(i)
+          write(iout,'(i5,f15.10)') k,v1bend_chyb(k,i)
+        enddo
+        enddo
+      endif
 C Read of Side-chain backbone correlation parameters
 C Modified 11 May 2012 by Adasko
 CCC
@@ -869,7 +964,7 @@ cc maxinter is maximum interaction sites
             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
      &(1+vlor3sccor(k,i,j)**2)
           enddo
-          v0sccor(i,j,iblock)=v0ijsccor
+          v0sccor(l,i,j)=v0ijsccor
         enddo
       enddo
       enddo
@@ -878,6 +973,7 @@ cc maxinter is maximum interaction sites
 #endif      
       if (lprint) then
         write (iout,'(/a/)') 'Torsional constants:'
+        do l=1,maxinter
         do i=1,nsccortyp
           do j=1,nsccortyp
             write (iout,*) 'ityp',i,' jtyp',j
@@ -892,6 +988,7 @@ cc maxinter is maximum interaction sites
             enddo
           enddo
         enddo
+        enddo
       endif
 
 C
         write (iout,*) "Coefficients of the cumulants"
       endif
       read (ifourier,*) nloctyp
+#ifdef NEWCORR
+      read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
+      read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
+      itype2loc(ntyp1)=nloctyp
+      iloctyp(nloctyp)=ntyp1
+#else
+      do i=1,ntyp1
+        itype2loc(i)=itortyp(i)
+      enddo
+      iloctyp(0)=10
+      iloctyp(1)=9
+      iloctyp(2)=20
+      iloctyp(3)=ntyp1
+#endif
+      do i=1,ntyp1
+        itype2loc(-i)=-itype2loc(i)
+      enddo
+      do i=1,nloctyp
+        iloctyp(-i)=-iloctyp(i)
+      enddo
+      write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
+      write (iout,*) "nloctyp",nloctyp,
+     &  " iloctyp",(iloctyp(i),i=0,nloctyp)
       do i=0,nloctyp-1
         read (ifourier,*,end=115,err=115)
-        read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
+        read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
+#ifdef NEWCORR
+        read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
+        read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
+        read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
+        read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
+        read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
+#endif 
         if (lprint) then
         write (iout,*) 'Type',i
-        write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
+        write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
         endif
-        B1(1,i)  = b(3)
-        B1(2,i)  = b(5)
-        B1(1,-i) = b(3)
-        B1(2,-i) = -b(5)
+c        B1(1,i)  = b(3)
+c        B1(2,i)  = b(5)
+c        B1(1,-i) = b(3)
+c        B1(2,-i) = -b(5)
 c        b1(1,i)=0.0d0
 c        b1(2,i)=0.0d0
-        B1tilde(1,i) = b(3)
-        B1tilde(2,i) =-b(5)
-        B1tilde(1,-i) =-b(3)
-        B1tilde(2,-i) =b(5)
+c        B1tilde(1,i) = b(3)
+c        B1tilde(2,i) =-b(5)
+c        B1tilde(1,-i) =-b(3)
+c        B1tilde(2,-i) =b(5)
 c        b1tilde(1,i)=0.0d0
 c        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)
+c        B2(1,i)  = b(2)
+c        B2(2,i)  = b(4)
+c        B2(1,-i)  =b(2)
+c        B2(2,-i)  =-b(4)
+cc        B1tilde(1,i) = b(3,i)
+cc        B1tilde(2,i) =-b(5,i)
+C        B1tilde(1,-i) =-b(3,i)
+C        B1tilde(2,-i) =b(5,i)
+cc        b1tilde(1,i)=0.0d0
+cc        b1tilde(2,i)=0.0d0
+cc        B2(1,i)  = b(2,i)
+cc        B2(2,i)  = b(4,i)
+C        B2(1,-i)  =b(2,i)
+C        B2(2,-i)  =-b(4,i)
 
 c        b2(1,i)=0.0d0
 c        b2(2,i)=0.0d0
-        CC(1,1,i)= b(7)
-        CC(2,2,i)=-b(7)
-        CC(2,1,i)= b(9)
-        CC(1,2,i)= b(9)
-        CC(1,1,-i)= b(7)
-        CC(2,2,-i)=-b(7)
-        CC(2,1,-i)=-b(9)
-        CC(1,2,-i)=-b(9)
+        CC(1,1,i)= b(7,i)
+        CC(2,2,i)=-b(7,i)
+        CC(2,1,i)= b(9,i)
+        CC(1,2,i)= b(9,i)
+        CC(1,1,-i)= b(7,i)
+        CC(2,2,-i)=-b(7,i)
+        CC(2,1,-i)=-b(9,i)
+        CC(1,2,-i)=-b(9,i)
 c        CC(1,1,i)=0.0d0
 c        CC(2,2,i)=0.0d0
 c        CC(2,1,i)=0.0d0
 c        CC(1,2,i)=0.0d0
-        Ctilde(1,1,i)=b(7)
-        Ctilde(1,2,i)=b(9)
-        Ctilde(2,1,i)=-b(9)
-        Ctilde(2,2,i)=b(7)
-        Ctilde(1,1,-i)=b(7)
-        Ctilde(1,2,-i)=-b(9)
-        Ctilde(2,1,-i)=b(9)
-        Ctilde(2,2,-i)=b(7)
+        Ctilde(1,1,i)=b(7,i)
+        Ctilde(1,2,i)=b(9,i)
+        Ctilde(2,1,i)=-b(9,i)
+        Ctilde(2,2,i)=b(7,i)
+        Ctilde(1,1,-i)=b(7,i)
+        Ctilde(1,2,-i)=-b(9,i)
+        Ctilde(2,1,-i)=b(9,i)
+        Ctilde(2,2,-i)=b(7,i)
 
 c        Ctilde(1,1,i)=0.0d0
 c        Ctilde(1,2,i)=0.0d0
 c        Ctilde(2,1,i)=0.0d0
 c        Ctilde(2,2,i)=0.0d0
-        DD(1,1,i)= b(6)
-        DD(2,2,i)=-b(6)
-        DD(2,1,i)= b(8)
-        DD(1,2,i)= b(8)
-        DD(1,1,-i)= b(6)
-        DD(2,2,-i)=-b(6)
-        DD(2,1,-i)=-b(8)
-        DD(1,2,-i)=-b(8)
+        DD(1,1,i)= b(6,i)
+        DD(2,2,i)=-b(6,i)
+        DD(2,1,i)= b(8,i)
+        DD(1,2,i)= b(8,i)
+        DD(1,1,-i)= b(6,i)
+        DD(2,2,-i)=-b(6,i)
+        DD(2,1,-i)=-b(8,i)
+        DD(1,2,-i)=-b(8,i)
 c        DD(1,1,i)=0.0d0
 c        DD(2,2,i)=0.0d0
 c        DD(2,1,i)=0.0d0
 c        DD(1,2,i)=0.0d0
-        Dtilde(1,1,i)=b(6)
-        Dtilde(1,2,i)=b(8)
-        Dtilde(2,1,i)=-b(8)
-        Dtilde(2,2,i)=b(6)
-        Dtilde(1,1,-i)=b(6)
-        Dtilde(1,2,-i)=-b(8)
-        Dtilde(2,1,-i)=b(8)
-        Dtilde(2,2,-i)=b(6)
+        Dtilde(1,1,i)=b(6,i)
+        Dtilde(1,2,i)=b(8,i)
+        Dtilde(2,1,i)=-b(8,i)
+        Dtilde(2,2,i)=b(6,i)
+        Dtilde(1,1,-i)=b(6,i)
+        Dtilde(1,2,-i)=-b(8,i)
+        Dtilde(2,1,-i)=b(8,i)
+        Dtilde(2,2,-i)=b(6,i)
 
 c        Dtilde(1,1,i)=0.0d0
 c        Dtilde(1,2,i)=0.0d0
 c        Dtilde(2,1,i)=0.0d0
 c        Dtilde(2,2,i)=0.0d0
-        EE(1,1,i)= b(10)+b(11)
-        EE(2,2,i)=-b(10)+b(11)
-        EE(2,1,i)= b(12)-b(13)
-        EE(1,2,i)= b(12)+b(13)
-        EE(1,1,-i)= b(10)+b(11)
-        EE(2,2,-i)=-b(10)+b(11)
-        EE(2,1,-i)=-b(12)+b(13)
-        EE(1,2,-i)=-b(12)-b(13)
-
+        EEold(1,1,i)= b(10,i)+b(11,i)
+        EEold(2,2,i)=-b(10,i)+b(11,i)
+        EEold(2,1,i)= b(12,i)-b(13,i)
+        EEold(1,2,i)= b(12,i)+b(13,i)
+        EEold(1,1,-i)= b(10,i)+b(11,i)
+        EEold(2,2,-i)=-b(10,i)+b(11,i)
+        EEold(2,1,-i)=-b(12,i)+b(13,i)
+        EEold(1,2,-i)=-b(12,i)-b(13,i)
+        write(iout,*) "TU DOCHODZE"
+        print *,"JESTEM"
 c        ee(1,1,i)=1.0d0
 c        ee(2,2,i)=1.0d0
 c        ee(2,1,i)=0.0d0
 c        ee(1,2,i)=0.0d0
 c        ee(2,1,i)=ee(1,2,i)
       enddo
+c      lprint=.true.
       if (lprint) then
       do i=1,nloctyp
         write (iout,*) 'Type',i
@@ -1011,10 +1150,11 @@ c        ee(2,1,i)=ee(1,2,i)
         enddo
         write(iout,*) 'EE'
         do j=1,2
-          write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
+          write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
         enddo
       enddo
       endif
+c      lprint=.false.
 
 C 
 C Read electrostatic-interaction parameters
@@ -1061,7 +1201,7 @@ C
      & ', exponents are ',expon,2*expon 
       goto (10,20,30,30,40) ipot
 C----------------------- LJ potential ---------------------------------
-   10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
+   10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
      &   (sigma0(i),i=1,ntyp)
       if (lprint) then
        write (iout,'(/a/)') 'Parameters of the LJ potential:'
@@ -1073,7 +1213,7 @@ C----------------------- LJ potential ---------------------------------
       endif
       goto 50
 C----------------------- LJK potential --------------------------------
-   20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
+   20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
       if (lprint) then
        write (iout,'(/a/)') 'Parameters of the LJK potential:'
@@ -1087,12 +1227,23 @@ C----------------------- LJK potential --------------------------------
       goto 50
 C---------------------- GB or BP potential -----------------------------
    30 do i=1,ntyp
-       read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
+       read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
       enddo
       read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
       read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
       read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
       read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
+C now we start reading lipid
+      do i=1,ntyp
+       read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
+       write(iout,*) "epslip", i, (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
+      write(iout,*) epslip(1,1),"OK?"
 C For the GB potential convert sigma'**2 into chi'
       if (ipot.eq.4) then
        do i=1,ntyp
@@ -1111,7 +1262,7 @@ C For the GB potential convert sigma'**2 into chi'
       endif
       goto 50
 C--------------------- GBV potential -----------------------------------
-   40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
+   40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
       if (lprint) then
@@ -1131,6 +1282,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
@@ -1159,10 +1311,17 @@ 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)
+          epsijlip=epslip(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
@@ -1195,7 +1354,7 @@ c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
           endif
          if (lprint) then
             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
-     &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
+     &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
          endif
         enddo
@@ -1252,7 +1411,7 @@ c      lprint=.false.
 C
 C Define the constants of the disulfide bridge
 C
-      ebr=-5.50D0
+C      ebr=-12.00D0
 c
 c Old arbitrary potential - commented out.
 c
@@ -1263,13 +1422,13 @@ c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 c energy surface of diethyl disulfide.
 c A. Liwo and U. Kozlowska, 11/24/03
 c
-      D0CM = 3.78d0
-      AKCM = 15.1d0
-      AKTH = 11.0d0
-      AKCT = 12.0d0
-      V1SS =-1.08d0
-      V2SS = 7.61d0
-      V3SS = 13.7d0
+C      D0CM = 3.78d0
+C      AKCM = 15.1d0
+C      AKTH = 11.0d0
+C      AKCT = 12.0d0
+C      V1SS =-1.08d0
+C      V2SS = 7.61d0
+C      V3SS = 13.7d0
 c      akcm=0.0d0
 c      akth=0.0d0
 c      akct=0.0d0
@@ -1277,14 +1436,33 @@ c      v1ss=0.0d0
 c      v2ss=0.0d0
 c      v3ss=0.0d0
       
-      if(me.eq.king) then
-      write (iout,'(/a)') "Disulfide bridge parameters:"
-      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
-      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
-      write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
-      write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
-     &  ' v3ss:',v3ss
-      endif
+C      if(me.eq.king) then
+C      write (iout,'(/a)') "Disulfide bridge parameters:"
+C      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
+C      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
+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
@@ -1299,11 +1477,15 @@ c      v3ss=0.0d0
       goto 999
   116 write (iout,*) "Error reading electrostatic energy parameters."
       goto 999
+ 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
+      goto 999
   117 write (iout,*) "Error reading side chain interaction parameters."
       goto 999
   118 write (iout,*) "Error reading SCp interaction parameters."
       goto 999
   119 write (iout,*) "Error reading SCCOR parameters"
+      go to 999
+  121 write (iout,*) "Error in Czybyshev parameters"
   999 continue
 #ifdef MPI
       call MPI_Finalize(Ierror)
@@ -1354,6 +1536,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