# Set comipiler flags for different sourcefiles
#================================================
if (Fortran_COMPILER_NAME STREQUAL "ifort")
- set(FFLAGS0 "-ip -w -I. -I${CMAKE_CURRENT_SOURCE_DIR}/include_unres" )
+ set(FFLAGS0 "-mcmodel=medium -shared-intel -ip -w -I. -I${CMAKE_CURRENT_SOURCE_DIR}/include_unres" )
elseif (Fortran_COMPILER_NAME STREQUAL "gfortran")
set(FFLAGS0 "-std=legacy -I. -I${CMAKE_CURRENT_SOURCE_DIR}/include_unres" )
else ()
#=========================================
if (Fortran_COMPILER_NAME STREQUAL "ifort")
# Add ifort preprocessor flags
- set(CPPFLAGS "${CPPFLAGS} -DPGI")
+ set(CPPFLAGS "${CPPFLAGS} -DPGI" )
elseif (Fortran_COMPILER_NAME STREQUAL "f95")
# Add new gfortran flags
set(CPPFLAGS "${CPPFLAGS} -DG77")
& nstart_seq,
& nend_sup, chain_length,tabperm(maxperm,maxsym)
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
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,scalscp,cutoff_corr,delt_corr,
- & r0_corr
+ & r0_corr,wliptran
integer ipot,n_ene_comp,rescale_mode
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),scalscp,
+ & wvdwpp,wbond,wliptran,weights(max_ene),scalscp,
& scal14,cutoff_corr,delt_corr,r0_corr,ipot,n_ene_comp,
& rescale_mode
common /potentials/ potname(5)
C General I/O units & files
integer inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,irotam,
& itorp,itordp,ifourier,ielep,isidep,iscpp,icbase,istat,
- & ientin,ientout,isidep1,ibond,isccor,jrms,jplot
+ & ientin,ientout,isidep1,ibond,isccor,jrms,jplot,
+ & iliptranpar
common /iounits/ inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,
& irotam,itorp,itordp,ifourier,ielep,isidep,iscpp,icbase,
- & istat,ientin,ientout,isidep1,ibond,isccor,jrms,jplot
+ & istat,ientin,ientout,isidep1,ibond,isccor,jrms,jplot,
+ & iliptranpar
character*256 outname,intname,pdbname,mol2name,statname,intinname,
& entname,restartname,prefix,scratchdir,sidepname,pdbfile,
& sccorname,rmsname,prefintin,prefout
& icsa_bank_reminimized,icsa_native_int,icsa_in
C Parameter files
character*256 bondname,thetname,rotname,torname,tordname,
- & fouriername,elename,sidename,scpname,patname
+ & fouriername,elename,sidename,scpname,patname,liptranname
common /parfiles/ thetname,rotname,torname,tordname,bondname,
- & fouriername,elename,sidename,scpname,patname
+ & fouriername,elename,sidename,scpname,patname,liptranname
character*3 pot
C-----------------------------------------------------------------------
C INP - main input file
parameter (maxprocs=16)
C Max. number of AA residues
integer maxres,maxres2
- parameter (maxres=100)
+ parameter (maxres=800)
C Appr. max. number of interaction sites
parameter (maxres2=2*maxres)
C Max. number of variables
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
if (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) then
call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
endif
-C write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
+ write (iout,*) "ft(6)",fact(6),wliptran,eliptran
#ifdef SPLITELE
if (shield_mode.gt.0) then
etot=fact(1)*wsc*(evdw+fact(6)*evdw_t)+fact(1)*wscp*evdw2
& +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+ethetacnstr
-C & +wliptran*eliptran
+ & +wliptran*eliptran
else
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees
& +wvdwpp*evdw1
& +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+ethetacnstr
-C & +wliptran*eliptran
+ & +wliptran*eliptran
endif
#else
if (shield_mode.gt.0) then
& +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+ethetacnstr
-C & +wliptran*eliptran
+ & +wliptran*eliptran
else
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
& +welec*fact(1)*(ees+evdw1)
& +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+ethetacnstr
-C & +wliptran*eliptran
+ & +wliptran*eliptran
endif
#endif
energia(20)=edihcnstr
energia(21)=evdw_t
energia(24)=ethetacnstr
+ energia(22)=eliptran
c detecting NaNQ
#ifdef ISNAN
#ifdef AIX
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)
c ROZNICA z WHAM
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
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)
cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
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
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
evdwij=evdwij*eps2rt*eps3rt
ij=icant(itypi,itypj)
aux=eps1*eps2rt**2*eps3rt**2
- 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
cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
cd & restyp(itypi),i,restyp(itypj),j,
cd & epsi,sigm,chi1,chi2,chip1,chip2,
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)
dzi=dc_norm(3,nres+i)
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,*) "czy jest 0", aa-aa_lip(itypi,itypj),
+ & aa-aa_aq(itypi,itypj)
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
yj_safe=yj
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
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
+C#define DEBUG
#ifdef DEBUG
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& evdwij
write (iout,*) "pratial sum", evdw,evdw_t
#endif
+C#undef DEBUG
c endif
if (calc_grad) then
C Calculate gradient components.
sigder=fac*sigder
fac=rij*fac
fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
+ gg_lipi(3)=eps1*(eps2rt*eps2rt)
+ &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
C Calculate the radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
dsci_inv=vbld_inv(i+nres)
+C returning the ith atom to box
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ 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
C
C Calculate SC interaction energy.
C
c alf1=0.0D0
c alf2=0.0D0
c alf12=0.0D0
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+C returning jth atom to box
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ 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
+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
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+C finding the closest
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
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
gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
enddo
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
& +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
& +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
- gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipi(k)
& +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
& +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
enddo
C
do k=i,j-1
do l=1,3
- gvdwc(l,k)=gvdwc(l,k)+gg(l)
+ gvdwc(l,k)=gvdwc(l,k)+gg(l)+gg_lipi(l)
enddo
enddo
+ do l=1,3
+ gvdwc(l,j)=gvdwc(l,j)+gg_lipj(l)
+ enddo
return
end
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-----------------------------------------------------------------------
+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
+ write(iout,*) "I am in?"
+ 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
+C-------------------------------------------------------------------------------------
& faceps1,faceps1_inv,eps1_om12,facsig,sigsq,sigsq_om1,sigsq_om2,
& sigsq_om12,facp,facp_inv,facp1,eps2rt,eps2rt_om1,eps2rt_om2,
& eps2rt_om12,eps3rt,eom1,eom2,eom12,evdwij,eps2der,eps3der,sigder,
- & dsci_inv,dscj_inv,gg
+ & dsci_inv,dscj_inv,gg,gg_lipi,gg_lipj
common /calc/ erij(3),rij,xj,yj,zj,dxi,dyi,dzi,dxj,dyj,dzj,
& chi1,chi2,chi12,chip1,chip2,chip12,alf1,alf2,alf12,om1,om2,om12,
& om1om2,chiom1,chiom2,chiom12,chipom1,chipom2,chipom12,eps1,
& faceps1,faceps1_inv,eps1_om12,facsig,sigsq,sigsq_om1,sigsq_om2,
& sigsq_om12,facp,facp_inv,facp1,eps2rt,eps2rt_om1,eps2rt_om2,
& eps2rt_om12,eps3rt,eom1,eom2,eom12,evdwij,eps2der,eps3der,sigder,
- & dsci_inv,dscj_inv,gg(3),i,j
+ & dsci_inv,dscj_inv,gg(3),gg_lipi(3),gg_lipj(3),i,j
- 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,
& iscpstart(maxres,maxint_gr),iscpend(maxres,maxint_gr),
& iatsc_s,iatsc_e,iatel_s,iatel_e,iatscp_s,iatscp_e,ispp,iscp
C 12/1/95 Array EPS included in the COMMON block.
- double precision eps,sigma,sigmaii,rs0,chi,chip,chip0,alp,signa0,
+ double precision eps,epslip,sigma,sigmaii,rs0,chi,chip,chip0,
+ & alp,signa0,
& sigii,sigma0,rr0,r0,r0e,r0d,rpp,epp,elpp6,elpp3,eps_scp,rscp,
& eps_orig
- common /body/eps(ntyp,ntyp),sigma(ntyp,ntyp),sigmaii(ntyp,ntyp),
+ common /body/eps(ntyp,ntyp),epslip(ntyp,ntyp),
+ & sigma(ntyp,ntyp),sigmaii(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),
& 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)
+
ibond=28
isccor=29
jrms=30
+ iliptran=60
C
C Set default weights of the energy terms.
C
enddo
do i=1,ntyp
do j=1,ntyp
- aa(i,j)=0.0D0
- bb(i,j)=0.0D0
+ aa_aq(i,j)=0.0D0
+ bb_aq(i,j)=0.0D0
+ aa_lip(i,j)=0.0D0
+ bb_lip(i,j)=0.0D0
augm(i,j)=0.0D0
sigma(i,j)=0.0D0
r0(i,j)=0.0D0
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
read (isidep,*)(chip0(i),i=1,ntyp)
read (isidep,*)(alp(i),i=1,ntyp)
C For the GB potential convert sigma'**2 into chi'
+ do i=1,ntyp
+ read (isidep,*)(epslip(i,j),j=i,ntyp)
+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
+ enddo
if (ipot.eq.4) then
do i=1,ntyp
chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
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
r0(j,i)=rrij
rrij=rrij**expon
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)
- if (ipot.gt.2) then
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ 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
sigii1=sigii(i)
endif
if (lprint) then
write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
- & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
+ & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
& sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
endif
enddo
c write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
c call enerprint(energia(0),fT)
c call pdbout(totfree(i),16,i)
-C#define DEBUG
+#define DEBUG
#ifdef DEBUG
write (iout,*) i," energia",(energia(j),j=0,max_ene)
write (iout,*) "etot", etot
write (iout,*) "ft(6)", ft(6)
#endif
-C#undef DEBUG
+#undef DEBUG
do k=1,max_ene
enetb(k,i)=energia(k)
enddo
chalen=int((nct-nnt+2)/symetr)
call int_from_cart1(.false.)
do j=nnt+1,nct
- if (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) then
+ if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)
+ & .and.(itype(j).ne.ntyp1)) then
if (j.gt.2) then
if (itel(j).ne.0 .and. itel(j-1).ne.0) then
write (iout,*) "Conformation",jjj,jj+1
enddo
do j=nnt,nct
itj=itype(j)
- if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0)
- & then
+ if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0
+ & .and. itype(j).ne.ntyp1) then
write (iout,*) "Conformation",jjj,jj+1
write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)
write (iout,*) "The Cartesian geometry is:"
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
C Shielding mode
call readi(controlcard,'SHIELD',shield_mode,0)
write (iout,*) "SHIELD MODE",shield_mode
call reada(weightcard,"EBR",ebr,-5.50D0)
call reada(weightcard,'WSHIELD',wshield,1.0d0)
write(iout,*) 'WSHIELD',wshield
+ call reada(weightcard,'WLT',wliptran,0.0D0)
call reada(weightcard,"ATRISS",atriss,0.301D0)
call reada(weightcard,"BTRISS",btriss,0.021D0)
call reada(weightcard,"CTRISS",ctriss,1.001D0)
open (isidep1,file=sidepname,status="old")
call getenv('SCCORPAR',sccorname)
open (isccor,file=sccorname,status="old")
+ call getenv('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old')
#ifndef OLDSCP
C
C 8/9/01 In the newest version SCp interaction constants are read from a file
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
dsci_inv=vbld_inv(i+nres)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ xi=dmod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=dmod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=dmod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+C define scaling factor for lipids
+
+C 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 ((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-
+ & ((positi-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-positi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
itypj=itype(j)
- xj=c(1,nres+j)-c(1,nres+i)
- yj=c(2,nres+j)-c(2,nres+i)
- zj=c(3,nres+j)-c(3,nres+i)
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(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-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/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
+
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
ljXs=sig-sig0ij
ljA=eps1*eps2rt**2*eps3rt**2
- ljB=ljA*bb(itypi,itypj)
- ljA=ljA*aa(itypi,itypj)
- ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+ ljB=ljA*bb
+ ljA=ljA*aa
+ ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
ssXs=d0cm
deltat1=1.0d0-om1
c Stop and plot energy and derivative as a function of distance
if (checkstop) then
ssm=ssC-0.25D0*ssB*ssB/ssA
- ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
+ ljm=-0.25D0*ljB*bb/aa
if (ssm.lt.ljm .and.
& dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
nicheck=1000
havebond=.false.
ljd=rij-ljXs
fac=(1.0D0/ljd)**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
eij=eps1*eps2rt*eps3rt*(e1+e2)
C write(iout,*) eij,'TU?1'
eps2der=eij*eps3rt
eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
else
havebond=.false.
- ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
- d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB
+ ljm=-0.25D0*ljB*bb/aa
+ d_ljm(1)=-0.5D0*bb/aa*ljB
d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
+ alf12/eps3rt)
num_cont_hb(i)=num_conti
enddo
do i=iturn4_start,iturn4_end
- if (i.le.1) cycle
+ if (i.lt.1) cycle
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)
#ifdef MPL
c endif
#endif
-C#define DEBUG
+#define DEBUG
#ifdef DEBUG
call enerprint(energia,fact)
#endif
-C#undef DEBUG
+#undef DEBUG
if (calc_grad) then
C
C Sum up the components of the Cartesian gradient.
projects / unres.git / commitdiff