Corrected WHAM and UNRES reading in SS dyn...

diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F
index 3296e45..c841df2 100644 (file)
--- a/source/unres/src_MD-M/parmread.F
+++ b/source/unres/src_MD-M/parmread.F
@@ -1256,7 +1256,7 @@ c      lprint=.false.
 C
 C Define the constants of the disulfide bridge
 C
-      ebr=-12.00D0
+C      ebr=-12.00D0
 c
 c Old arbitrary potential - commented out.
 c

diff --git a/source/wham/src-M/geomout.F b/source/wham/src-M/geomout.F
index c9110ae..6747e76 100644 (file)
--- a/source/wham/src-M/geomout.F
+++ b/source/wham/src-M/geomout.F
@@ -54,11 +54,11 @@
         write (ipdb,30) ica(nct),ica(nct)+1
       endif
       do i=1,nss
-       if (dyn_ss) then
-        write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
-       else
-        write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
-       endif
+C       if (dyn_ss) then
+C        write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
+C       else
+C        write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
+C       endif
       enddo
       write (ipdb,'(a6)') 'ENDMDL'
   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)

diff --git a/source/wham/src-M/include_unres/COMMON.SBRIDGE b/source/wham/src-M/include_unres/COMMON.SBRIDGE
index 144a864..1380548 100644 (file)
--- a/source/wham/src-M/include_unres/COMMON.SBRIDGE
+++ b/source/wham/src-M/include_unres/COMMON.SBRIDGE
@@ -12,7 +12,7 @@
       common /links_split/ link_start,link_end
       double precision Ht,dyn_ssbond_ij,dtriss,atriss,btriss,ctriss
       logical dyn_ss,dyn_ss_mask
-      common /dyn_ssbond/ dtriss,atriss,btriss,ctriss,
+      common /dyn_ssbond/ dtriss,atriss,btriss,ctriss,Ht,
      &  dyn_ssbond_ij(maxres,maxres),
      &  idssb(maxdim),jdssb(maxdim),
-     &  Ht,dyn_ss,dyn_ss_mask(maxres)
+     &  dyn_ss,dyn_ss_mask(maxres)

diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F
index d50cd03..750f03e 100644 (file)
--- a/source/wham/src-M/parmread.F
+++ b/source/wham/src-M/parmread.F
@@ -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,14 @@ 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)
@@ -63,24 +72,45 @@ c If reading not own parameters, skip assignment
 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)
-      if (dyn_ss) then
-        ss_depth=ebr/wsc-0.25*eps(1,1)
-        Ht=Ht/wsc-0.25*eps(1,1)
-        akcm=akcm*wstrain/wsc
-        akth=akth*wstrain/wsc
-        akct=akct*wstrain/wsc
-        v1ss=v1ss*wstrain/wsc
-        v2ss=v2ss*wstrain/wsc
-        v3ss=v3ss*wstrain/wsc
-      else
-        ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
-      endif
+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
 
@@ -104,7 +134,7 @@ c
       wvdwpp=ww(16)
       wbond=ww(18)
       wsccor=ww(19)
-
+      whpb=ww(15)
       endif
 
       call card_concat(controlcard,.false.)
@@ -1217,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
@@ -1228,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

diff --git a/source/wham/src-M/ssMD.F b/source/wham/src-M/ssMD.F
index 283adf3..5080b18 100644 (file)
--- a/source/wham/src-M/ssMD.F
+++ b/source/wham/src-M/ssMD.F
@@ -251,6 +251,7 @@ c-------END TESTING CODE
         e1=fac*fac*aa(itypi,itypj)
         e2=fac*bb(itypi,itypj)
         eij=eps1*eps2rt*eps3rt*(e1+e2)
+C        write(iout,*) eij,'TU?1'
         eps2der=eij*eps3rt
         eps3der=eij*eps2rt
         eij=eij*eps2rt*eps3rt
@@ -267,7 +268,7 @@ c-------END TESTING CODE
         havebond=.true.
         ssd=rij-ssXs
         eij=ssA*ssd*ssd+ssB*ssd+ssC
-
+C        write(iout,*) 'TU?2',ssc,ssd
         ed=2*akcm*ssd+akct*deltat12
         pom1=akct*ssd
         pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
@@ -303,6 +304,7 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
           h1=h_base(f1,hd1)
           h2=h_base(f2,hd2)
           eij=ssm*h1+Ht*h2
+C         write(iout,*) eij,'TU?3'
           delta_inv=1.0d0/(xm-ssxm)
           deltasq_inv=delta_inv*delta_inv
           fac=ssm*hd1-Ht*hd2
@@ -325,6 +327,7 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
           h1=h_base(f1,hd1)
           h2=h_base(f2,hd2)
           eij=Ht*h1+ljm*h2
+C        write(iout,*) 'TU?4',ssA
           delta_inv=1.0d0/(ljxm-xm)
           deltasq_inv=delta_inv*delta_inv
           fac=Ht*hd1-ljm*hd2
@@ -396,7 +399,7 @@ c$$$        if (ed.gt.0.0d0) havebond=.true.
 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
 
       endif
-
+      write(iout,*) 'havebond',havebond
       if (havebond) then
 #ifndef CLUST
 #ifndef WHAM