working cluster for nano parameters

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

diff --git a/source/cluster/wham/src-M/parmread.F b/source/cluster/wham/src-M/parmread.F
index 4292bd6..26502e7 100644 (file)
--- a/source/cluster/wham/src-M/parmread.F
+++ b/source/cluster/wham/src-M/parmread.F
@@ -72,6 +72,11 @@ C Assign virtual-bond length
           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 
@@ -870,6 +875,14 @@ C---------------------- GB or BP potential -----------------------------
       read (isidep,*)(chip0(i),i=1,ntyp)
       read (isidep,*)(alp(i),i=1,ntyp)
 C For the GB potential convert sigma'**2 into chi'
+      do i=1,ntyp
+       read (isidep,*)(epslip(i,j),j=i,ntyp)
+C       write(iout,*) "WARNING!!",i,ntyp
+       write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
+C       do j=1,ntyp
+C       epslip(i,j)=epslip(i,j)+0.05d0
+C       enddo
+      enddo
       if (ipot.eq.4) then
        do i=1,ntyp
          chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
@@ -907,6 +920,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
@@ -924,6 +938,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
@@ -933,13 +948,19 @@ C Calculate the "working" parameters of SC interactions.
          r0(j,i)=rrij
          rrij=rrij**expon
          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)
-         if (ipot.gt.2) then
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          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
            sigii1=sigii(i)
@@ -971,7 +992,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
@@ -979,6 +1000,26 @@ c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 C
 C Define the SC-p interaction constants
 C
+      write(iout,*) "tube param"
+      read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep,
+     & ccavtubpep,dcavtubpep,tubetranenepep
+      sigmapeptube=sigmapeptube**6
+      sigeps=dsign(1.0D0,epspeptube)
+      epspeptube=dabs(epspeptube)
+      pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
+      pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
+      write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
+      do i=1,ntyp
+       read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i),
+     & ccavtub(i),dcavtub(i),tubetranene(i)
+       sigmasctube=sigmasctube**6
+       sigeps=dsign(1.0D0,epssctube)
+       epssctube=dabs(epssctube)
+       sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
+       sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
+      write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
+      enddo
+
       do i=1,20
         do j=1,2
           eps_scp(i,j)=-1.5d0