From d305127baa0f7e25ecb5422e6e0e34aa15db4776 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Wed, 3 Feb 2016 12:55:22 +0100 Subject: [PATCH] correction in wham and UNRES for lipid and correlation --- ..._ext.1gab_3S_qclass5no310-shan2-sc-16-10-8k_lip | 10 ++-- PARAM/scinter_GB_ext_lip.parm | 59 ++++++++++---------- source/unres/src_MD-M/energy_p_new_barrier.F | 6 +- source/unres/src_MD-M/parmread.F | 2 + source/wham/src-M/energy_p_new.F | 56 ++++++++++++------- source/wham/src-M/parmread.F | 3 +- source/wham/src-M/wham_calc1.F | 2 +- 7 files changed, 81 insertions(+), 57 deletions(-) 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 --- 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 - -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 - -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 - -1.447000 -28.07000 3.360749 4.211903 + -0.1779 −3.4493 3.360749 4.211903 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 diff --git a/PARAM/scinter_GB_ext_lip.parm b/PARAM/scinter_GB_ext_lip.parm index c508b3e..e600c30 100644 --- 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 + 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 - - -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 --- 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 +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??" @@ -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)) +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 @@ -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 (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 @@ -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 + 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 (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 diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index e96dfc3..2d2c23b 100644 --- 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) + 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 diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index eb6146a..2b99394 100644 --- 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 - integer icant + integer icant,xshift,yshift,zshift 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 +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 @@ -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 - if (i.eq.1) then +C if (i.eq.1) then 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 - endif +C endif 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) - 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 - & .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 - endif +C endif 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 + 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 @@ -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 - & .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 - & .or. itype(i+5).eq.ntyp1 +C & .or. itype(i+5).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 @@ -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 write(iout,*) i,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 --- 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) -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 diff --git a/source/wham/src-M/wham_calc1.F b/source/wham/src-M/wham_calc1.F index 77fc439..eff98da 100644 --- 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, - & (enetb(k,i,iparm),k=1,21) + & (enetb(k,i,iparm),k=1,22) #endif call restore_parm(iparm) #ifdef DEBUG -- 1.7.9.5