Corrected WHAM and UNRES reading in SS dyn...
[unres.git] / source / wham / src-M / parmread.F
index 84f725c..750f03e 100644 (file)
@@ -35,6 +35,7 @@ C
       character*16 key
       integer iparm
       double precision ip,mp
+C      write (iout,*) "KURWA"
 C
 C Body
 C
@@ -55,6 +56,61 @@ C Assign virtual-bond length
 
       write (iout,*) "iparm",iparm," myparm",myparm
 c If reading not own parameters, skip assignment
+      call reada(weightcard,"D0CM",d0cm,3.78d0)
+      call reada(weightcard,"AKCM",akcm,15.1d0)
+      call reada(weightcard,"AKTH",akth,11.0d0)
+      call reada(weightcard,"AKCT",akct,12.0d0)
+      call reada(weightcard,"V1SS",v1ss,-1.08d0)
+      call reada(weightcard,"V2SS",v2ss,7.61d0)
+      call reada(weightcard,"V3SS",v3ss,13.7d0)
+      call reada(weightcard,"EBR",ebr,-5.50D0)
+      call reada(controlcard,"DTRISS",dtriss,1.0D0)
+      call reada(controlcard,"ATRISS",atriss,0.3D0)
+      call reada(controlcard,"BTRISS",btriss,0.02D0)
+      call reada(controlcard,"CTRISS",ctriss,1.0D0)
+      dyn_ss=(index(controlcard,'DYN_SS').gt.0)
+C      do i=1,maxres
+C        dyn_ss_mask(i)=.false.
+C      enddo
+C      ebr=-12.0D0
+c
+c Old arbitrary potential - commented out.
+c
+c      dbr= 4.20D0
+c      fbr= 3.30D0
+c
+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
+
+      do i=1,maxres-1
+        do j=i+1,maxres
+          dyn_ssbond_ij(i,j)=1.0d300
+        enddo
+      enddo
+      call reada(controlcard,"HT",Ht,0.0D0)
+C      if (dyn_ss) then
+C        ss_depth=ebr/wsc-0.25*eps(1,1)
+C        write(iout,*) HT,wsc,eps(1,1),'KURWA'
+C        Ht=Ht/wsc-0.25*eps(1,1)
+       
+C        akcm=akcm*whpb/wsc
+C        akth=akth*whpb/wsc
+C        akct=akct*whpb/wsc
+C        v1ss=v1ss*whpb/wsc
+C        v2ss=v2ss*whpb/wsc
+C        v3ss=v3ss*whpb/wsc
+C      else
+C        ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
+C      endif
 
       if (iparm.eq.myparm .or. .not.separate_parset) then
 
@@ -78,7 +134,7 @@ c
       wvdwpp=ww(16)
       wbond=ww(18)
       wsccor=ww(19)
-
+      whpb=ww(15)
       endif
 
       call card_concat(controlcard,.false.)
@@ -295,7 +351,7 @@ C
 C Read the parameters of Utheta determined from ab initio surfaces
 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 C
-      write (iout,*) "tu dochodze"
+c      write (iout,*) "tu dochodze"
       read (ithep,*) nthetyp,ntheterm,ntheterm2,
      &  ntheterm3,nsingle,ndouble
       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
@@ -303,7 +359,7 @@ C
       do i=-ntyp1,-1
         ithetyp(i)=-ithetyp(-i)
       enddo
-      write (iout,*) "tu dochodze"
+c      write (iout,*) "tu dochodze"
       do iblock=1,2
       do i=-maxthetyp,maxthetyp
         do j=-maxthetyp,maxthetyp
@@ -982,8 +1038,10 @@ C
         bpp (i,j)=-2.0D0*epp(i,j)*rri
         ael6(i,j)=elpp6(i,j)*4.2D0**6
         ael3(i,j)=elpp3(i,j)*4.2D0**3
+        lprint=.true.
         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
      &                    ael6(i,j),ael3(i,j)
+         lprint=.false.
         enddo
       enddo
 C
@@ -1024,13 +1082,17 @@ C----------------------- LJK potential --------------------------------
       endif
       goto 50
 C---------------------- GB or BP potential -----------------------------
-   30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
-     &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
-     &  (alp(i),i=1,ntyp)
+   30 do i=1,ntyp
+       read (isidep,*)(eps(i,j),j=i,ntyp)
+      enddo
+      read (isidep,*)(sigma0(i),i=1,ntyp)
+      read (isidep,*)(sigii(i),i=1,ntyp)
+      read (isidep,*)(chip(i),i=1,ntyp)
+      read (isidep,*)(alp(i),i=1,ntyp)
 C For the GB potential convert sigma'**2 into chi'
       if (ipot.eq.4) then
        do i=1,ntyp
-         chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
+         chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
         enddo
       endif
       if (lprint) then
@@ -1185,7 +1247,7 @@ C
 C
 C Define the constants of the disulfide bridge
 C
-      ebr=-5.50D0
+C      ebr=-12.0D0
 c
 c Old arbitrary potential - commented out.
 c
@@ -1196,21 +1258,36 @@ 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
 
-      if (lprint) then
+      if (dyn_ss) then
+        ss_depth=ebr/wsc-0.25*eps(1,1)
+C        write(iout,*) akcm,whpb,wsc,'KURWA'
+        Ht=Ht/wsc-0.25*eps(1,1)
+
+        akcm=akcm*whpb/wsc
+        akth=akth*whpb/wsc
+        akct=akct*whpb/wsc
+        v1ss=v1ss*whpb/wsc
+        v2ss=v2ss*whpb/wsc
+        v3ss=v3ss*whpb/wsc
+      else
+        ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
+      endif
+
+C      if (lprint) 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      endif
       return
       end