From: Adam Sieradzan Date: Tue, 29 Sep 2015 11:02:05 +0000 (+0200) Subject: wham in lipid still diff X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=48f00408a2e488d2c34fd65d26fc69b214b4fbb2 wham in lipid still diff --- diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 6ab2a4e..889e5a5 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -233,7 +233,7 @@ c Cutoff range for interactions bordliptop=(boxzsize+lipthick)/2.0 bordlipbot=bordliptop-lipthick C endif - if ((bordliptop.gt.boxzsize).or.(borlipbot.lt.0.0)) + if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" buflipbot=bordlipbot+lipbufthick bufliptop=bordliptop-lipbufthick diff --git a/source/wham/src-M/COMMON.ALLPARM b/source/wham/src-M/COMMON.ALLPARM index 88b480f..5b6284a 100644 --- a/source/wham/src-M/COMMON.ALLPARM +++ b/source/wham/src-M/COMMON.ALLPARM @@ -66,8 +66,10 @@ & app_all(2,2,max_parm),bpp_all(2,2,max_parm), & ael6_all(2,2,max_parm),ael3_all(2,2,max_parm), & aad_all(ntyp,2,max_parm),bad_all(ntyp,2,max_parm), - & aa_all(ntyp,ntyp,max_parm),bb_all(ntyp,ntyp,max_parm), + & aa_aq_all(ntyp,ntyp,max_parm),bb_aq_all(ntyp,ntyp,max_parm), + & aa_lip_all(ntyp,ntyp,max_parm),bb_lip_all(ntyp,ntyp,max_parm), & augm_all(ntyp,ntyp,max_parm),eps_all(ntyp,ntyp,max_parm), + & epslip_all(ntyp,ntyp,max_parm), & sigma_all(ntyp,ntyp,max_parm),r0_all(ntyp,ntyp,max_parm), & chi_all(ntyp,ntyp,max_parm),chip_all(ntyp,max_parm), & alp_all(ntyp,max_parm),ebr_all(max_parm),d0cm_all(max_parm), @@ -98,7 +100,8 @@ & v0_all,v1_all,v2_all,vlor1_all,vlor2_all,vlor3_all,v1c_all, & v1s_all,v2c_all,v2s_all,b1_all,b2_all,cc_all,dd_all,ee_all, & ctilde_all,dtilde_all,b1tilde_all,app_all,bpp_all,ael6_all, - & ael3_all,aad_all,bad_all,aa_all,bb_all,augm_all, + & ael3_all,aad_all,bad_all,aa_aq_all,bb_aq_all,augm_all, + & aa_lip_all,bb_lip_all,epslip_all, & eps_all,sigma_all,r0_all,chi_all,chip_all,alp_all,ebr_all, & d0cm_all,akcm_all,akth_all,akct_all,v1ss_all,v2ss_all,v3ss_all, & v1sccor_all,v2sccor_all,nbondterm_all, diff --git a/source/wham/src-M/COMMON.CHAIN b/source/wham/src-M/COMMON.CHAIN index b147f08..24b8c56 100644 --- a/source/wham/src-M/COMMON.CHAIN +++ b/source/wham/src-M/COMMON.CHAIN @@ -13,6 +13,8 @@ &nend_sup,chain_length,tabperm(maxperm,maxsym),nperm, & nstart_seq,ishift_pdb double precision boxxsize,boxysize,boxzsize,enecut,sscut,sss, - &sssgrad - common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad + &sssgrad, + & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick + common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad, + & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick diff --git a/source/wham/src-M/COMMON.IOUNITS b/source/wham/src-M/COMMON.IOUNITS index 23783bb..188d55e 100644 --- a/source/wham/src-M/COMMON.IOUNITS +++ b/source/wham/src-M/COMMON.IOUNITS @@ -10,11 +10,12 @@ C----------------------------------------------------------------------- C General I/O units & files integer inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,irotam, & itorp,itordp,ifourier,ielep,isidep,iscpp,isccor,icbase, - & istat,ientin,ientout,isidep1,ibond,ihist,izsc,idistr + & istat,ientin,ientout,isidep1,ibond,ihist,izsc,idistr, + & iliptranpar common /iounits/ inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep, & irotam,itorp,itordp,ifourier,ielep,isidep,iscpp,isccor, & icbase,istat,ientin,ientout,isidep1,ibond,ihist,izsc, - & idistr + & idistr,iliptranpar character*256 outname,intname,pdbname,mol2name,statname,intinname, & entname,restartname,prefix,scratchdir,sidepname,pdbfile, & histname,zscname @@ -23,9 +24,11 @@ C General I/O units & files & sidepname,pdbfile,histname,zscname C Parameter files character*256 bondname,thetname,rotname,torname,tordname, - & fouriername,elename,sidename,scpname,sccorname,patname + & fouriername,elename,sidename,scpname,sccorname,patname, + & liptranname common /parfiles/ thetname,rotname,torname,tordname,bondname, - & fouriername,elename,sidename,scpname,sccorname,patname + & fouriername,elename,sidename,scpname,sccorname,patname, + & liptranname character*3 pot C----------------------------------------------------------------------- C INP - main input file diff --git a/source/wham/src-M/DIMENSIONS.ZSCOPT b/source/wham/src-M/DIMENSIONS.ZSCOPT index 2f5ee76..950f76d 100644 --- a/source/wham/src-M/DIMENSIONS.ZSCOPT +++ b/source/wham/src-M/DIMENSIONS.ZSCOPT @@ -3,7 +3,7 @@ c Maximum number of structures in the database, energy components, proteins, c and structural classes c#ifdef JUBL - parameter (maxstr=200000,max_ene=21,maxprot=7,maxclass=5000) + parameter (maxstr=200000,max_ene=22,maxprot=7,maxclass=5000) parameter (maxclass1=10) c Maximum number of structures to be dealt with by one processor parameter (maxstr_proc=10000) diff --git a/source/wham/src-M/enecalc1.F b/source/wham/src-M/enecalc1.F index c29d517..81668b5 100644 --- a/source/wham/src-M/enecalc1.F +++ b/source/wham/src-M/enecalc1.F @@ -228,7 +228,7 @@ c call intout endif C write (iout,*) "Czy tu dochodze" potE(iii+1,iparm)=energia(0) - do k=1,21 + do k=1,22 enetb(k,iii+1,iparm)=energia(k) enddo #ifdef DEBUG diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index cede380..c228c20 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -90,6 +90,11 @@ C C 21/5/07 Calculate local sicdechain correlation energy C call eback_sc_corr(esccor) + + if (wliptran.gt.0) then + call Eliptransfer(eliptran) + endif + C C 12/1/95 Multi-body terms C @@ -114,7 +119,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor + & +wbond*estr+wsccor*fact(1)*esccor+wliptran*eliptran #else etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2 & +welec*fact(1)*(ees+evdw1) @@ -123,7 +128,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor + & +wbond*estr+wsccor*fact(1)*esccor+wliptran*eliptran #endif energia(0)=etot energia(1)=evdw @@ -157,6 +162,8 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t energia(19)=esccor energia(20)=edihcnstr energia(21)=evdw_t + energia(22)=eliptran + c detecting NaNQ #ifdef ISNAN #ifdef AIX @@ -195,10 +202,12 @@ C & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(2)*gsccorx(j,i) + & +wliptran*gliptranx(j,i) enddo #else do i=1,nct @@ -214,10 +223,12 @@ C & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(1)*gsccorx(j,i) + & +wliptran*gliptranx(j,i) enddo #endif enddo @@ -270,6 +281,7 @@ C------------------------------------------------------------------------ esccor=energia(19) edihcnstr=energia(20) estr=energia(18) + eliptran=energia(22) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1, & wvdwpp, @@ -278,7 +290,7 @@ C------------------------------------------------------------------------ & ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5), & eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2), & eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5), - & esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot + & esccor,wsccor*fact(1),edihcnstr,ebr*nss,eliptran,wliptran,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -301,6 +313,7 @@ C------------------------------------------------------------------------ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ + & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ & 'ETOT= ',1pE16.6,' (total)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),estr,wbond, @@ -309,7 +322,7 @@ C------------------------------------------------------------------------ & ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2), & eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3), & eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor, - & edihcnstr,ebr*nss,etot + & edihcnstr,ebr*nss,eliptran,wliptran,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -331,6 +344,7 @@ C------------------------------------------------------------------------ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ + & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -396,8 +410,8 @@ C Change 12/1/95 to calculate four-body interactions c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj eps0ij=eps(itypi,itypj) fac=rrij**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=e1+e2 ij=icant(itypi,itypj) eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij) @@ -408,7 +422,7 @@ cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)') cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij else evdw_t=evdw_t+evdwij @@ -566,8 +580,8 @@ C rij=1.0D0/r_inv_ij r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=e_augm+e1+e2 ij=icant(itypi,itypj) eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm) @@ -580,7 +594,7 @@ cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm, cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij else evdw_t=evdw_t+evdwij @@ -712,8 +726,8 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives. C Calculate whole angle-dependent part of epsilon and contributions C to its derivatives fac=(rrij*sigsq)**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt @@ -723,15 +737,15 @@ C to its derivatives eneps_temp(1,ij)=eneps_temp(1,ij)+e1*aux & /dabs(eps(itypi,itypj)) eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj) - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij else evdw_t=evdw_t+evdwij endif if (calc_grad) then if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa write (iout,'(2(a3,i3,2x),15(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, & epsi,sigm,chi1,chi2,chip1,chip2, @@ -807,6 +821,28 @@ C returning the ith atom to box if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize + if ((zi.gt.bordlipbot) + &.and.(zi.lt.bordliptop)) then +C the energy transfer exist + if (zi.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0 + endif dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) @@ -851,6 +887,33 @@ C returning jth atom to box if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize + if ((zj.gt.bordlipbot) + &.and.(zj.lt.bordliptop)) then +C the energy transfer exist + if (zj.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zj-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zj.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) + sslipj=sscalelip(fracinbuf) + ssgradlipj=sscagradlip(fracinbuf)/lipbufthick + else + sslipj=1.0d0 + ssgradlipj=0.0 + endif + else + sslipj=0.0d0 + ssgradlipj=0.0 + endif + aa=aa_lip(itypi,itypj)*(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 + 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 xj_safe=xj @@ -896,6 +959,7 @@ c write (iout,*) i,j,xj,yj,zj if (sss.le.0.0) cycle C Calculate angle-dependent terms of energy and contributions to their C derivatives. + call sc_angular sigsq=1.0D0/sigsq sig=sig0ij*dsqrt(sigsq) @@ -909,13 +973,13 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt evdwij=evdwij*eps2rt*eps3rt - if (bb(itypi,itypj).gt.0) then + if (bb.gt.0) then evdw=evdw+evdwij*sss else evdw_t=evdw_t+evdwij*sss @@ -929,8 +993,8 @@ c write (iout,*) "i",i," j",j," itypi",itypi," itypj",itypj, c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)), c & aux*e2/eps(itypi,itypj) c if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa #ifdef DEBUG write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, @@ -1060,15 +1124,15 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt fac_augm=rrij**expon e_augm=augm(itypi,itypj)*fac_augm evdwij=evdwij*eps2rt*eps3rt - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij+e_augm else evdw_t=evdw_t+evdwij+e_augm @@ -7596,6 +7660,125 @@ cd write (2,*) 'eel_turn6',ekont*eel_turn6 return end crc------------------------------------------------- +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + subroutine Eliptransfer(eliptran) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' + include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' + include 'COMMON.SBRIDGE' +C this is done by Adasko +C print *,"wchodze" +C structure of box: +C water +C--bordliptop-- buffore starts +C--bufliptop--- here true lipid starts +C lipid +C--buflipbot--- lipid ends buffore starts +C--bordlipbot--buffore ends + eliptran=0.0 + do i=1,nres +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 +C print *,i +C first for peptide groups +c for each residue check if it is in lipid or lipid water border area + if ((positi.gt.bordlipbot) + &.and.(positi.lt.bordliptop)) then +C the energy transfer exist + if (positi.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((positi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslip=sscalelip(fracinbuf) + ssgradlip=-sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*pepliptran + gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0 + gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0 +C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran + elseif (positi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + sslip=sscalelip(fracinbuf) + ssgradlip=sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*pepliptran + gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0 + gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0 +C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran +C print *, "doing sscalefor top part" +C print *,i,sslip,fracinbuf,ssgradlip + else + eliptran=eliptran+pepliptran +C print *,"I am in true lipid" + endif +C else +C eliptran=elpitran+0.0 ! I am in water + endif + enddo +C print *, "nic nie bylo w lipidzie?" +C now multiply all by the peptide group transfer factor +C eliptran=eliptran*pepliptran +C now the same for side chains +CV do i=1,1 + do i=1,nres + if (itype(i).eq.ntyp1) cycle + positi=(mod(c(3,i+nres),boxzsize)) + if (positi.le.0) positi=positi+boxzsize +C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop +c for each residue check if it is in lipid or lipid water border area +C respos=mod(c(3,i+nres),boxzsize) +C print *,positi,bordlipbot,buflipbot + if ((positi.gt.bordlipbot) + & .and.(positi.lt.bordliptop)) then +C the energy transfer exist + if (positi.lt.buflipbot) then + fracinbuf=1.0d0- + & ((positi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslip=sscalelip(fracinbuf) + ssgradlip=-sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*liptranene(itype(i)) + gliptranx(3,i)=gliptranx(3,i) + &+ssgradlip*liptranene(itype(i)) + gliptranc(3,i-1)= gliptranc(3,i-1) + &+ssgradlip*liptranene(itype(i)) +C print *,"doing sccale for lower part" + elseif (positi.gt.bufliptop) then + fracinbuf=1.0d0- + &((bordliptop-positi)/lipbufthick) + sslip=sscalelip(fracinbuf) + ssgradlip=sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*liptranene(itype(i)) + gliptranx(3,i)=gliptranx(3,i) + &+ssgradlip*liptranene(itype(i)) + gliptranc(3,i-1)= gliptranc(3,i-1) + &+ssgradlip*liptranene(itype(i)) +C print *, "doing sscalefor top part",sslip,fracinbuf + else + eliptran=eliptran+liptranene(itype(i)) +C print *,"I am in true lipid" + endif + endif ! if in lipid or buffor +C else +C eliptran=elpitran+0.0 ! I am in water + enddo + return + end + + +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + SUBROUTINE MATVEC2(A1,V1,V2) implicit real*8 (a-h,o-z) include 'DIMENSIONS' @@ -7759,4 +7942,34 @@ C----------------------------------------------------------------------- return end C----------------------------------------------------------------------- +C----------------------------------------------------------------------- + double precision function sscalelip(r) + double precision r,gamm + include "COMMON.SPLITELE" +C if(r.lt.r_cut-rlamb) then +C sscale=1.0d0 +C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then +C gamm=(r-(r_cut-rlamb))/rlamb + sscalelip=1.0d0+r*r*(2*r-3.0d0) +C else +C sscale=0d0 +C endif + return + end +C----------------------------------------------------------------------- + double precision function sscagradlip(r) + double precision r,gamm + include "COMMON.SPLITELE" +C if(r.lt.r_cut-rlamb) then +C sscagrad=0.0d0 +C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then +C gamm=(r-(r_cut-rlamb))/rlamb + sscagradlip=r*(6*r-6.0d0) +C else +C sscagrad=0.0d0 +C endif + return + end + +C----------------------------------------------------------------------- diff --git a/source/wham/src-M/include_unres/COMMON.DERIV b/source/wham/src-M/include_unres/COMMON.DERIV index 79f8630..80cbed6 100644 --- a/source/wham/src-M/include_unres/COMMON.DERIV +++ b/source/wham/src-M/include_unres/COMMON.DERIV @@ -1,5 +1,6 @@ double precision dcdv,dxdv,dxds,gradx,gradc,gvdwc,gelc,gvdwpp, & gradx_scp,gvdwc_scp,ghpbx,ghpbc,gloc,gvdwx,gradcorr,gradxorr, + & gliptranc,gliptranx, & gradcorr5,gradcorr6,gel_loc,gcorr3_turn,gcorr4_turn,gcorr6_turn, & gel_loc_loc,gel_loc_turn3,gel_loc_turn4,gel_loc_turn6,gcorr_loc, & g_corr5_loc,g_corr6_loc,gradb,gradbx,gsccorc,gsccorx,gsccor_loc, @@ -10,6 +11,8 @@ & gradx(3,maxres,2),gradc(3,maxres,2),gvdwx(3,maxres), & gvdwc(3,maxres),gelc(3,maxres),gvdwpp(3,maxres), & gradx_scp(3,maxres), + & gliptranc(3,-1:maxres), + & gliptranx(3,-1:maxres), & gvdwc_scp(3,maxres),ghpbx(3,maxres),ghpbc(3,maxres), & gloc(maxvar,2),gradcorr(3,maxres),gradxorr(3,maxres), & gradcorr5(3,maxres),gradcorr6(3,maxres), diff --git a/source/wham/src-M/include_unres/COMMON.FFIELD b/source/wham/src-M/include_unres/COMMON.FFIELD index 0c169f7..6c432a9 100644 --- a/source/wham/src-M/include_unres/COMMON.FFIELD +++ b/source/wham/src-M/include_unres/COMMON.FFIELD @@ -6,11 +6,11 @@ C----------------------------------------------------------------------- double precision wsc,wscp,welec,wstrain,wtor,wtor_d,wang,wscloc, & wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4, & wturn6,wvdwpp,wbond,weights,scal14,cutoff_corr,delt_corr, - & r0_corr + & r0_corr,wliptran integer ipot,n_ene_comp common /ffield/ wsc,wscp,welec,wstrain,wtor,wtor_d,wang,wscloc, & wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4, - & wturn6,wvdwpp,wbond,weights(max_ene), + & wturn6,wvdwpp,wbond,wliptran,weights(max_ene), & scal14,cutoff_corr,delt_corr,r0_corr,ipot,n_ene_comp common /potentials/ potname(5) character*3 potname diff --git a/source/wham/src-M/include_unres/COMMON.INTERACT b/source/wham/src-M/include_unres/COMMON.INTERACT index d38c470..7d6b59f 100644 --- a/source/wham/src-M/include_unres/COMMON.INTERACT +++ b/source/wham/src-M/include_unres/COMMON.INTERACT @@ -1,8 +1,10 @@ - double precision aa,bb,augm,aad,bad,app,bpp,ael6,ael3 + double precision aa_aq,bb_aq,augm,aad,bad,app,bpp,ael6,ael3, + & aa_lip,bb_lip integer nnt,nct,nint_gr,istart,iend,itype,itel,itypro,ielstart, & ielend,nscp_gr,iscpstart,iscpend,iatsc_s,iatsc_e,iatel_s, & iatel_e,iatscp_s,iatscp_e,ispp,iscp,expon,expon2 - common /interact/aa(ntyp,ntyp),bb(ntyp,ntyp),augm(ntyp,ntyp), + common /interact/aa_aq(ntyp,ntyp),bb_aq(ntyp,ntyp), + & augm(ntyp,ntyp),aa_lip(ntyp,ntyp),bb_lip(ntyp,ntyp), & aad(ntyp,2),bad(ntyp,2),app(2,2),bpp(2,2),ael6(2,2),ael3(2,2), & expon,expon2,nnt,nct,nint_gr(maxres),istart(maxres,maxint_gr), & iend(maxres,maxint_gr),itype(maxres),itel(maxres),itypro, @@ -12,8 +14,9 @@ C 12/1/95 Array EPS included in the COMMON block. double precision eps,sigma,sigmaii,rs0,chi,chip,chip0,alp,signa0, & sigii,sigma0,rr0,r0,r0e,r0d,rpp,epp,elpp6,elpp3,eps_scp,rscp, - & eps_orig + & eps_orig,epslip common /body/eps(ntyp,ntyp),sigma(ntyp,ntyp),sigmaii(ntyp,ntyp), + &epslip(ntyp,ntyp), & rs0(ntyp,ntyp),chi(ntyp,ntyp),chip(ntyp),chip0(ntyp),alp(ntyp), & sigma0(ntyp),sigii(ntyp),rr0(ntyp),r0(ntyp,ntyp),r0e(ntyp,ntyp), & r0d(ntyp,2),rpp(2,2),epp(2,2),elpp6(2,2),elpp3(2,2), @@ -26,4 +29,8 @@ c 12/5/03 modified 09/18/03 Bond stretching parameters. & aksc(maxbondterm,ntyp),abond0(maxbondterm,ntyp), & distchainmax,nbondterm(ntyp) &,vbldpDUM +C 01/29/15 Lipidic parameters + double precision pepliptran,liptranene + common /lipid/ pepliptran,liptranene(ntyp) + diff --git a/source/wham/src-M/initialize_p.F b/source/wham/src-M/initialize_p.F index d1cc94b..f47a951 100644 --- a/source/wham/src-M/initialize_p.F +++ b/source/wham/src-M/initialize_p.F @@ -62,6 +62,8 @@ C ihist=30 iweight=31 izsc=32 +C Lipidic input file for parameters range 60-79 + iliptranpar=60 C C Set default weights of the energy terms. C @@ -88,8 +90,10 @@ C enddo do i=1,ntyp do j=1,ntyp - aa(i,j)=0.0D0 - bb(i,j)=0.0D0 + aa_lip(i,j)=0.0D0 + bb_lip(i,j)=0.0D0 + aa_aq(i,j)=0.0D0 + bb_aq(i,j)=0.0D0 augm(i,j)=0.0D0 sigma(i,j)=0.0D0 r0(i,j)=0.0D0 @@ -261,17 +265,17 @@ c------------------------------------------------------------------------- & "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ", & "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ", & "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP", - & "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T"/ + & "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T","ELT"/ data wname / & "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC", & "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD", - & "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC"/ + & "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC","WLT"/ data ww0 /1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0, & 1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,0.4d0,1.0d0,1.0d0, - & 0.0d0,0.0/ - data nprint_ene /21/ + & 0.0d0,0.0,0.0/ + data nprint_ene /22/ data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19, - & 16,15,17,20,21/ + & 16,15,17,20,21,22/ end c--------------------------------------------------------------------------- subroutine init_int_table diff --git a/source/wham/src-M/make_ensemble1.F b/source/wham/src-M/make_ensemble1.F index 5d7b750..699995e 100644 --- a/source/wham/src-M/make_ensemble1.F +++ b/source/wham/src-M/make_ensemble1.F @@ -23,7 +23,7 @@ double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl, & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/ double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors, - & escloc, + & escloc,eliptran, & ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3, & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist @@ -162,6 +162,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft estr=enetb(18,i,iparm) esccor=enetb(19,i,iparm) edihcnstr=enetb(20,i,iparm) + eliptran=enetb(22,i,iparm) #ifdef SPLITELE etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees & +wvdwpp*evdw1 @@ -171,7 +172,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor - & +wbond*estr + & +wbond*estr+wliptran*eliptran #else etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 & +ft(1)*welec*(ees+evdw1) @@ -181,7 +182,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor - & +wbond*estr + & +wbond*estr+wliptran*eliptran #endif #ifdef MPI Fdimless(i)= diff --git a/source/wham/src-M/openunits.F b/source/wham/src-M/openunits.F index b9f54b7..5200b3e 100644 --- a/source/wham/src-M/openunits.F +++ b/source/wham/src-M/openunits.F @@ -48,6 +48,8 @@ C Get parameter filenames and open the parameter files. open (isidep,file=sidename,status='old') call mygetenv('SIDEP',sidepname) open (isidep1,file=sidepname,status="old") + call mygetenv('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old',action='read') #ifndef OLDSCP C C 8/9/01 In the newest version SCp interaction constants are read from a file diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F index ea6cb14..914c090 100644 --- a/source/wham/src-M/parmread.F +++ b/source/wham/src-M/parmread.F @@ -78,7 +78,7 @@ c wvdwpp=ww(16) wbond=ww(18) wsccor=ww(19) - + wliptran=ww(22) endif call card_concat(controlcard,.false.) @@ -188,6 +188,11 @@ c 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 @@ -1033,6 +1038,13 @@ C---------------------- GB or BP potential ----------------------------- read (isidep,*)(sigii(i),i=1,ntyp) read (isidep,*)(chip(i),i=1,ntyp) read (isidep,*)(alp(i),i=1,ntyp) + do i=1,ntyp + read (isidep,*)(epslip(i,j),j=i,ntyp) +C print *,"WARNING!!" +C do j=1,ntyp +C epslip(i,j)=epslip(i,j)+0.05d0 +C enddo + enddo C For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp @@ -1071,6 +1083,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 @@ -1088,6 +1101,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 @@ -1099,10 +1113,16 @@ C Calculate the "working" parameters of SC interactions. 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) + 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 diff --git a/source/wham/src-M/readrtns.F b/source/wham/src-M/readrtns.F index 9023170..276e7a6 100644 --- a/source/wham/src-M/readrtns.F +++ b/source/wham/src-M/readrtns.F @@ -81,6 +81,23 @@ c Cutoff range for interactions call reada(controlcard,"R_CUT",r_cut,15.0d0) call reada(controlcard,"LAMBDA",rlamb,0.3d0) + call reada(controlcard,"LIPTHICK",lipthick,0.0d0) + call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) + if (lipthick.gt.0.0d0) then + bordliptop=(boxzsize+lipthick)/2.0 + bordlipbot=bordliptop-lipthick +C endif + if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) + & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" + buflipbot=bordlipbot+lipbufthick + bufliptop=bordliptop-lipbufthick + if ((lipbufthick*2.0d0).gt.lipthick) + &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF" + endif + write(iout,*) "bordliptop=",bordliptop + write(iout,*) "bordlipbot=",bordlipbot + write(iout,*) "bufliptop=",bufliptop + write(iout,*) "buflipbot=",buflipbot call readi(controlcard,'SYM',symetr,1) write (iout,*) "DISTCHAINMAX",distchainmax write (iout,*) "delta",delta diff --git a/source/wham/src-M/store_parm.F b/source/wham/src-M/store_parm.F index 8c44422..f44b6f4 100644 --- a/source/wham/src-M/store_parm.F +++ b/source/wham/src-M/store_parm.F @@ -40,6 +40,7 @@ c Store weights ww_all(16,iparm)=wvdwpp ww_all(17,iparm)=wbond ww_all(19,iparm)=wsccor + ww_all(22,iparm)=wliptran c Store bond parameters vbldp0_all(iparm)=vbldp0 akp_all(iparm)=akp @@ -222,13 +223,16 @@ c Store the parameters of electrostatic interactions c Store sidechain parameters do i=1,ntyp do j=1,ntyp - aa_all(j,i,iparm)=aa(j,i) - bb_all(j,i,iparm)=bb(j,i) + aa_aq_all(j,i,iparm)=aa_aq(j,i) + bb_aq_all(j,i,iparm)=bb_aq(j,i) + aa_lip_all(j,i,iparm)=aa_lip(j,i) + bb_lip_all(j,i,iparm)=bb_lip(j,i) r0_all(j,i,iparm)=r0(j,i) sigma_all(j,i,iparm)=sigma(j,i) chi_all(j,i,iparm)=chi(j,i) augm_all(j,i,iparm)=augm(j,i) eps_all(j,i,iparm)=eps(j,i) + epslip_all(j,i,iparm)=epslip(j,i) enddo enddo do i=1,ntyp @@ -311,6 +315,7 @@ c Restore weights wvdwpp=ww_all(16,iparm) wbond=ww_all(17,iparm) wsccor=ww_all(19,iparm) + wliptran=ww_all(22,iparm) c Restore bond parameters vbldp0=vbldp0_all(iparm) akp=akp_all(iparm) @@ -492,13 +497,16 @@ c Restore the parameters of electrostatic interactions c Restore sidechain parameters do i=1,ntyp do j=1,ntyp - aa(j,i)=aa_all(j,i,iparm) - bb(j,i)=bb_all(j,i,iparm) + aa_aq(j,i)=aa_aq_all(j,i,iparm) + bb_aq(j,i)=bb_aq_all(j,i,iparm) + aa_lip(j,i)=aa_lip_all(j,i,iparm) + bb_lip(j,i)=bb_lip_all(j,i,iparm) r0(j,i)=r0_all(j,i,iparm) sigma(j,i)=sigma_all(j,i,iparm) chi(j,i)=chi_all(j,i,iparm) augm(j,i)=augm_all(j,i,iparm) eps(j,i)=eps_all(j,i,iparm) + epslip(j,i)=epslip_all(j,i,iparm) enddo enddo do i=1,ntyp diff --git a/source/wham/src-M/wham_calc1.F b/source/wham/src-M/wham_calc1.F index 15d6716..45573d1 100644 --- a/source/wham/src-M/wham_calc1.F +++ b/source/wham/src-M/wham_calc1.F @@ -84,7 +84,8 @@ c parameter (MaxHdim=200000) & eplus,eminus,logfac,tanhT,tt double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors, & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3, - & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor + & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor, + & eliptran integer ind_point(maxpoint),upindE,indE character*16 plik @@ -312,6 +313,8 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft estr=enetb(18,i,iparm) esccor=enetb(19,i,iparm) edihcnstr=enetb(20,i,iparm) + eliptran=enetb(22,i,iparm) + #ifdef DEBUG write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6), & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc, @@ -327,7 +330,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor - & +wbond*estr + & +wbond*estr+wliptran*eliptran #else etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 & +ft(1)*welec*(ees+evdw1) @@ -337,7 +340,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor - & +wbond*estr + & +wbond*estr+wliptran*eliptran #endif #ifdef DEBUG write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3), @@ -776,7 +779,7 @@ c write (iout,*) "ftbis",ftbis & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor - & +wbond*estr + & +wbond*estr+wliptran*eliptran eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees & +ftprim(1)*wtor*etors+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ @@ -799,7 +802,7 @@ c write (iout,*) "ftbis",ftbis & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor - & +wbond*estr + & +wbond*estr+wliptran*eliptran eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) & +ftprim(1)*wtor*etors+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+