correction in wham and UNRES for lipid and correlation

diff --git a/PARAM/sc_GB_opt_ext.1gab_3S_qclass5no310-shan2-sc-16-10-8k_lip b/PARAM/sc_GB_opt_ext.1gab_3S_qclass5no310-shan2-sc-16-10-8k_lip
index 3adbfc8..9708c27 100644 (file)
--- a/PARAM/sc_GB_opt_ext.1gab_3S_qclass5no310-shan2-sc-16-10-8k_lip
+++ b/PARAM/sc_GB_opt_ext.1gab_3S_qclass5no310-shan2-sc-16-10-8k_lip
@@ -214,21 +214,21 @@
    3.992146       3.451319       2.676042       3.353784
    4.172640       2.402115       2.402115
 
    3.992146       3.451319       2.676042       3.353784
    4.172640       2.402115       2.402115
 

-  -61.60000      -35.26000       3.800125       48.12500
-   35.00000       2.916757       3.655463       4.547977
+  -7.2938        -4.0239         3.800125       5.9101
+   3.8743         2.916757       3.655463       4.547977

    2.618191       2.618191
 
    2.618191       2.618191
 

-  -23.62500       3.505321       39.90000       28.44000
+  -3.0057         3.505321       6.0636         3.5303

    2.690483       3.371882       4.195157       2.415078
    2.415078
 
    3.824000       4.378585       3.785405       2.935082
    3.678429       4.576550       2.634640       2.634640
 
    2.690483       3.371882       4.195157       2.415078
    2.415078
 
    3.824000       4.378585       3.785405       2.935082
    3.678429       4.576550       2.634640       2.634640
 

-  -1.447000      -28.07000       3.360749       4.211903
+  -0.1779        −3.4493         3.360749       4.211903

    5.240276       3.016734       3.016734
 
    5.240276       3.016734       3.016734
 

-  -76.82000       2.905459       3.641304       4.530360
+  -6.5557         2.905459       3.641304       4.530360

    2.608049       2.608049
 
    2.252800       2.823350       3.512697       2.022197
    2.608049       2.608049
 
    2.252800       2.823350       3.512697       2.022197

diff --git a/PARAM/scinter_GB_ext_lip.parm b/PARAM/scinter_GB_ext_lip.parm
index c508b3e..e600c30 100644 (file)
--- a/PARAM/scinter_GB_ext_lip.parm
+++ b/PARAM/scinter_GB_ext_lip.parm
@@ -134,6 +134,7 @@
          .0572167103        -.0468608825         .0151048455         .0084963678
          .0278930397         .0076922911         .1033536738        -.0098256036
          .0611385674         .0448303346         .0861379100         .0861379100
          .0572167103        -.0468608825         .0151048455         .0084963678
          .0278930397         .0076922911         .1033536738        -.0098256036
          .0611385674         .0448303346         .0861379100         .0861379100

+

    2.252000       2.758289       2.828747       2.573807    
    2.573807       2.314337       3.329045       3.123838    
    1.679138       1.245463       2.417517       2.138542    
    2.252000       2.758289       2.828747       2.573807    
    2.573807       2.314337       3.329045       3.123838    
    1.679138       1.245463       2.417517       2.138542    

@@ -213,32 +214,32 @@
    3.178800       3.464740       3.195954       3.486507    
    3.992146       3.451319       2.676042       3.353784    
    4.172640       2.402115       2.402115    
    3.178800       3.464740       3.195954       3.486507    
    3.992146       3.451319       2.676042       3.353784    
    4.172640       2.402115       2.402115    

-  
-  -61.60000      -35.26000       3.800125       48.12500    
-   35.00000       2.916757       3.655463       4.547977    
-   2.618191       2.618191    
-  
-  -23.62500       3.505321       39.90000       28.44000    
-   2.690483       3.371882       4.195157       2.415078    
-   2.415078    
-  
-   3.824000       4.378585       3.785405       2.935082    
-   3.678429       4.576550       2.634640       2.634640    
-  
-  -1.447000      -28.07000       3.360749       4.211903    
-   5.240276       3.016734       3.016734    
-  
-  -76.82000       2.905459       3.641304       4.530360    
-   2.608049       2.608049    
-  
-   2.252800       2.823350       3.512697       2.022197    
-   2.022197    
-  
-   3.538400       4.402332       2.534345       2.534345    
-  
-   5.477200       3.153127       3.153127    
-  
-   1.815200       1.815200    
-  
-   1.815200    
-  
+
+   -7.2938        -4.0239         3.800125       5.9101
+   3.8743         2.916757       3.655463       4.547977
+   2.618191       2.618191
+
+  -3.0057         3.505321       6.0636         3.5303
+   2.690483       3.371882       4.195157       2.415078
+   2.415078
+
+   3.824000       4.378585       3.785405       2.935082
+   3.678429       4.576550       2.634640       2.634640
+
+  -0.1779        -3.4493         3.360749       4.211903
+   5.240276       3.016734       3.016734
+
+  -6.5557         2.905459       3.641304       4.530360
+   2.608049       2.608049
+
+   2.252800       2.823350       3.512697       2.022197
+   2.022197
+
+   3.538400       4.402332       2.534345       2.534345
+
+   5.477200       3.153127       3.153127
+
+   1.815200       1.815200
+
+   1.815200
+ 

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index c90fc3a..64288d6 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -1780,6 +1780,7 @@ C lipbufthick is thickenes of lipid buffore
      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0

+C      write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj)

 C      if (aa.ne.aa_aq(itypi,itypj)) write(63,'(2e10.5)')
 C     &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
 C      if (ssgradlipj.gt.0.0d0) print *,"??WTF??"
 C      if (aa.ne.aa_aq(itypi,itypj)) write(63,'(2e10.5)')
 C     &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
 C      if (ssgradlipj.gt.0.0d0) print *,"??WTF??"

@@ -2056,6 +2057,7 @@ C lipbufthick is thickenes of lipid buffore
      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
 C      if (aa.ne.aa_aq(itypi,itypj)) write(63,'2e10.5') 
 C     &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
 C      if (aa.ne.aa_aq(itypi,itypj)) write(63,'2e10.5') 
 C     &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))

+C      write(iout,*) "tu,", i,j,aa,bb,aa_lip(itypi,itypj),sslipi,sslipj

       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
       xj_safe=xj
       yj_safe=yj
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
       xj_safe=xj
       yj_safe=yj

@@ -5848,6 +5850,7 @@ C       Checking if it involves dummy (NH3+ or COO-) group
          if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
 C YES   vbldpDUM is the equlibrium length of spring for Dummy atom
         diff = vbld(i)-vbldpDUM
          if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
 C YES   vbldpDUM is the equlibrium length of spring for Dummy atom
         diff = vbld(i)-vbldpDUM

+        if (energy_dec) write(iout,*) "dum_bond",i,diff 

          else
 C NO    vbldp0 is the equlibrium lenght of spring for peptide group
         diff = vbld(i)-vbldp0
          else
 C NO    vbldp0 is the equlibrium lenght of spring for peptide group
         diff = vbld(i)-vbldp0

@@ -5861,6 +5864,7 @@ C NO    vbldp0 is the equlibrium lenght of spring for peptide group
 c        write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
 c        endif
       enddo
 c        write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
 c        endif
       enddo

+      

       estr=0.5d0*AKP*estr+estr1
 c
 c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
       estr=0.5d0*AKP*estr+estr1
 c
 c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included

@@ -11120,7 +11124,7 @@ C       do i=1,1
         if (itype(i).eq.ntyp1) cycle
 
         positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
         if (itype(i).eq.ntyp1) cycle
 
         positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))

-        if (positi.le.0) positi=positi+boxzsize
+        if (positi.le.0.0) positi=positi+boxzsize

 C        print *,i
 C first for peptide groups
 c for each residue check if it is in lipid or lipid water border area
 C        print *,i
 C first for peptide groups
 c for each residue check if it is in lipid or lipid water border area

diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F
index e96dfc3..2d2c23b 100644 (file)
--- a/source/unres/src_MD-M/parmread.F
+++ b/source/unres/src_MD-M/parmread.F
@@ -1236,6 +1236,8 @@ C---------------------- GB or BP potential -----------------------------
 C now we start reading lipid
       do i=1,ntyp
        read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,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       print *,"WARNING!!"
 C       do j=1,ntyp
 C       epslip(i,j)=epslip(i,j)+0.05d0

diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index eb6146a..2b99394 100644 (file)
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -895,7 +895,7 @@ C
       include 'COMMON.SBRIDGE'
       logical lprn
       common /srutu/icall
       include 'COMMON.SBRIDGE'
       logical lprn
       common /srutu/icall

-      integer icant
+      integer icant,xshift,yshift,zshift

       external icant
       do i=1,210
         do j=1,2
       external icant
       do i=1,210
         do j=1,2

@@ -1034,6 +1034,8 @@ C lipbufthick is thickenes of lipid buffore
      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0

+C      write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj)
+

 C        write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj)
 C checking the distance
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 C        write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj)
 C checking the distance
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2

@@ -2098,15 +2100,15 @@ C      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         gcorr_loc(i)=0.0d0
       enddo
       do i=iatel_s,iatel_e
         gcorr_loc(i)=0.0d0
       enddo
       do i=iatel_s,iatel_e

-          if (i.eq.1) then 
+C          if (i.eq.1) then 

            if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
            if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1

-     &  .or. itype(i+2).eq.ntyp1) cycle
-          else
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
-     &  .or. itype(i+2).eq.ntyp1
-     &  .or. itype(i-1).eq.ntyp1
+C     &  .or. itype(i+2).eq.ntyp1) cycle
+C          else
+C        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C     &  .or. itype(i+2).eq.ntyp1
+C     &  .or. itype(i-1).eq.ntyp1

      &) cycle
      &) cycle

-         endif
+C         endif

         if (itel(i).eq.0) goto 1215
         dxi=dc(1,i)
         dyi=dc(2,i)
         if (itel(i).eq.0) goto 1215
         dxi=dc(1,i)
         dyi=dc(2,i)

@@ -2126,16 +2128,16 @@ C      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         num_conti=0
 C        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
         num_conti=0
 C        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)

-          if (j.eq.1) then
-           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
-     & .or.itype(j+2).eq.ntyp1
-     &) cycle  
-          else     
+          if (j.le.1) cycle
+C           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
+C     & .or.itype(j+2).eq.ntyp1
+C     &) cycle  
+C          else     

           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1

-     & .or.itype(j+2).eq.ntyp1
-     & .or.itype(j-1).eq.ntyp1
+C     & .or.itype(j+2).eq.ntyp1
+C     & .or.itype(j-1).eq.ntyp1

      &) cycle
      &) cycle

-         endif
+C         endif

 C
 C) cycle
           if (itel(j).eq.0) goto 1216
 C
 C) cycle
           if (itel(j).eq.0) goto 1216

@@ -3028,6 +3030,18 @@ C Third- and fourth-order contributions from turns
      &    aggj(3,4),aggj1(3,4),a_temp(2,2)
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2
       if (j.eq.i+2) then
      &    aggj(3,4),aggj1(3,4),a_temp(2,2)
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2
       if (j.eq.i+2) then

+      if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+C     & .or.((i+5).gt.nres)
+C     & .or.((i-1).le.0)
+C end of changes suggested by Ana
+     &    .or. itype(i+2).eq.ntyp1
+     &    .or. itype(i+3).eq.ntyp1
+C     &    .or. itype(i+5).eq.ntyp1
+C     &    .or. itype(i).eq.ntyp1
+C     &    .or. itype(i-1).eq.ntyp1
+     &    ) goto 178
+

 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C               Third-order contributions
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C               Third-order contributions

@@ -3158,14 +3172,15 @@ C Cartesian derivatives
       else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then
       if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds
       else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then
       if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds

-     & .or.((i+5).gt.nres)
-     & .or.((i-1).le.0)
+C     & .or.((i+5).gt.nres)
+C     & .or.((i-1).le.0)

 C end of changes suggested by Ana
      &    .or. itype(i+3).eq.ntyp1
      &    .or. itype(i+4).eq.ntyp1
 C end of changes suggested by Ana
      &    .or. itype(i+3).eq.ntyp1
      &    .or. itype(i+4).eq.ntyp1

-     &    .or. itype(i+5).eq.ntyp1
+C     &    .or. itype(i+5).eq.ntyp1

      &    .or. itype(i).eq.ntyp1
      &    .or. itype(i).eq.ntyp1

-     &    .or. itype(i-1).eq.ntyp1) goto 178
+C     &    .or. itype(i-1).eq.ntyp1
+     &    ) goto 178

 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C               Fourth-order contributions
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C               Fourth-order contributions

@@ -3851,6 +3866,7 @@ C     &       gnmr1(vbld(i),-1.0d0,distchainmax)
 C        else
          if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
         diff = vbld(i)-vbldpDUM
 C        else
          if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
         diff = vbld(i)-vbldpDUM

+C         write(iout,*) i,diff

          else
           diff = vbld(i)-vbldp0
 c          write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
          else
           diff = vbld(i)-vbldp0
 c          write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff

diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F
index 097e1f7..8db78b9 100644 (file)
--- a/source/wham/src-M/parmread.F
+++ b/source/wham/src-M/parmread.F
@@ -1112,7 +1112,8 @@ C---------------------- GB or BP potential -----------------------------
       read (isidep,*)(alp(i),i=1,ntyp)
       do i=1,ntyp
        read (isidep,*)(epslip(i,j),j=i,ntyp)
       read (isidep,*)(alp(i),i=1,ntyp)
       do i=1,ntyp
        read (isidep,*)(epslip(i,j),j=i,ntyp)

-C       print *,"WARNING!!"
+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
 C       do j=1,ntyp
 C       epslip(i,j)=epslip(i,j)+0.05d0
 C       enddo

diff --git a/source/wham/src-M/wham_calc1.F b/source/wham/src-M/wham_calc1.F
index 77fc439..eff98da 100644 (file)
--- a/source/wham/src-M/wham_calc1.F
+++ b/source/wham/src-M/wham_calc1.F
@@ -228,7 +228,7 @@ c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
         do iparm=1,nParmSet
 #ifdef DEBUG
           write (iout,'(2i5,21f8.2)') i,iparm,
         do iparm=1,nParmSet
 #ifdef DEBUG
           write (iout,'(2i5,21f8.2)') i,iparm,

-     &     (enetb(k,i,iparm),k=1,21)
+     &     (enetb(k,i,iparm),k=1,22)

 #endif
           call restore_parm(iparm)
 #ifdef DEBUG
 #endif
           call restore_parm(iparm)
 #ifdef DEBUG