From 48f04f24e913a3e10867d2038b30efcd48a60a9f Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Tue, 8 Mar 2016 13:06:56 +0100 Subject: [PATCH] bug fix after Ana and cluster lipid (still in progress) --- source/cluster/wham/src-M/CMakeLists.txt | 4 +- source/cluster/wham/src-M/COMMON.CHAIN | 6 +- source/cluster/wham/src-M/COMMON.FFIELD | 4 +- source/cluster/wham/src-M/COMMON.IOUNITS | 10 +- source/cluster/wham/src-M/DIMENSIONS | 2 +- source/cluster/wham/src-M/energy_p_new.F | 369 ++++++++++++++++++-- .../cluster/wham/src-M/include_unres/COMMON.CALC | 4 +- .../wham/src-M/include_unres/COMMON.INTERACT | 16 +- source/cluster/wham/src-M/initialize_p.F | 7 +- source/cluster/wham/src-M/parmread.F | 36 +- source/cluster/wham/src-M/probabl.F | 4 +- source/cluster/wham/src-M/read_coords.F | 7 +- source/cluster/wham/src-M/readrtns.F | 20 ++ source/cluster/wham/src-M/ssMD.F | 128 ++++++- source/unres/src_MD-M/energy_p_new_barrier.F | 2 +- source/wham/src-M/energy_p_new.F | 4 +- 16 files changed, 547 insertions(+), 76 deletions(-) diff --git a/source/cluster/wham/src-M/CMakeLists.txt b/source/cluster/wham/src-M/CMakeLists.txt index 37d1ef8..30193dd 100644 --- a/source/cluster/wham/src-M/CMakeLists.txt +++ b/source/cluster/wham/src-M/CMakeLists.txt @@ -64,7 +64,7 @@ set(UNRES_CLUSTER_WHAM_M_PP_SRC # 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 () @@ -105,7 +105,7 @@ set(CPPFLAGS "${CPPFLAGS} -DUNRES -DISNAN") #========================================= 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") diff --git a/source/cluster/wham/src-M/COMMON.CHAIN b/source/cluster/wham/src-M/COMMON.CHAIN index 03b65c5..861df2f 100644 --- a/source/cluster/wham/src-M/COMMON.CHAIN +++ b/source/cluster/wham/src-M/COMMON.CHAIN @@ -11,6 +11,8 @@ & 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 diff --git a/source/cluster/wham/src-M/COMMON.FFIELD b/source/cluster/wham/src-M/COMMON.FFIELD index ccafd30..fa85436 100644 --- a/source/cluster/wham/src-M/COMMON.FFIELD +++ b/source/cluster/wham/src-M/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,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) diff --git a/source/cluster/wham/src-M/COMMON.IOUNITS b/source/cluster/wham/src-M/COMMON.IOUNITS index c97090d..d171ae0 100644 --- a/source/cluster/wham/src-M/COMMON.IOUNITS +++ b/source/cluster/wham/src-M/COMMON.IOUNITS @@ -10,10 +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,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 @@ -35,9 +37,9 @@ C CSA I/O units & files & 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 diff --git a/source/cluster/wham/src-M/DIMENSIONS b/source/cluster/wham/src-M/DIMENSIONS index 86f775a..a74fcfd 100644 --- a/source/cluster/wham/src-M/DIMENSIONS +++ b/source/cluster/wham/src-M/DIMENSIONS @@ -9,7 +9,7 @@ C Max. number of processors. 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 diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F index 6d5ec5a..27cfaac 100644 --- a/source/cluster/wham/src-M/energy_p_new.F +++ b/source/cluster/wham/src-M/energy_p_new.F @@ -96,6 +96,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 @@ -111,7 +116,7 @@ c print *,ecorr,ecorr5,ecorr6,eturn6 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 @@ -123,7 +128,7 @@ C write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t & +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 @@ -133,7 +138,7 @@ C & +wliptran*eliptran & +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 @@ -145,7 +150,7 @@ C & +wliptran*eliptran & +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) @@ -155,7 +160,7 @@ C & +wliptran*eliptran & +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 @@ -192,6 +197,7 @@ C & +wliptran*eliptran energia(20)=edihcnstr energia(21)=evdw_t energia(24)=ethetacnstr + energia(22)=eliptran c detecting NaNQ #ifdef ISNAN #ifdef AIX @@ -506,8 +512,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) c ROZNICA z WHAM @@ -521,7 +527,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 @@ -673,8 +679,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) cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) @@ -684,7 +690,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 @@ -810,23 +816,23 @@ 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 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, @@ -898,6 +904,28 @@ c if (icall.gt.0) lprn=.true. 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) @@ -972,6 +1000,34 @@ c alf12=0.0D0 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 @@ -1027,13 +1083,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 @@ -1044,8 +1100,9 @@ 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 +C#define DEBUG #ifdef DEBUG write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, @@ -1055,6 +1112,7 @@ c if (lprn) then & evdwij write (iout,*) "pratial sum", evdw,evdw_t #endif +C#undef DEBUG c endif if (calc_grad) then C Calculate gradient components. @@ -1063,6 +1121,12 @@ 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 @@ -1117,6 +1181,35 @@ c if (icall.gt.0) lprn=.true. 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 @@ -1147,9 +1240,76 @@ c chip12=0.0D0 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) @@ -1170,15 +1330,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 @@ -1284,10 +1444,10 @@ C---------------------------------------------------------------------------- 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 @@ -1296,9 +1456,12 @@ C Calculate the components of the gradient in DC and X 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------------------------------------------------------------------------------ @@ -8739,4 +8902,150 @@ C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) 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------------------------------------------------------------------------------------- diff --git a/source/cluster/wham/src-M/include_unres/COMMON.CALC b/source/cluster/wham/src-M/include_unres/COMMON.CALC index 67b4bb9..bf255c9 100644 --- a/source/cluster/wham/src-M/include_unres/COMMON.CALC +++ b/source/cluster/wham/src-M/include_unres/COMMON.CALC @@ -5,11 +5,11 @@ & 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 diff --git a/source/cluster/wham/src-M/include_unres/COMMON.INTERACT b/source/cluster/wham/src-M/include_unres/COMMON.INTERACT index 0db8895..1c0b8db 100644 --- a/source/cluster/wham/src-M/include_unres/COMMON.INTERACT +++ b/source/cluster/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, @@ -10,10 +12,12 @@ & 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), @@ -26,3 +30,7 @@ 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/cluster/wham/src-M/initialize_p.F b/source/cluster/wham/src-M/initialize_p.F index 7bf3f4d..409deb0 100644 --- a/source/cluster/wham/src-M/initialize_p.F +++ b/source/cluster/wham/src-M/initialize_p.F @@ -58,6 +58,7 @@ C ibond=28 isccor=29 jrms=30 + iliptran=60 C C Set default weights of the energy terms. C @@ -84,8 +85,10 @@ 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 diff --git a/source/cluster/wham/src-M/parmread.F b/source/cluster/wham/src-M/parmread.F index 4292bd6..66a7672 100644 --- a/source/cluster/wham/src-M/parmread.F +++ b/source/cluster/wham/src-M/parmread.F @@ -72,6 +72,11 @@ C Assign virtual-bond length 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 @@ -870,6 +875,14 @@ C---------------------- GB or BP potential ----------------------------- 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) @@ -907,6 +920,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 @@ -933,13 +947,19 @@ C Calculate the "working" parameters of SC interactions. 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) @@ -971,7 +991,7 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) 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 diff --git a/source/cluster/wham/src-M/probabl.F b/source/cluster/wham/src-M/probabl.F index 8d9ac48..7bc3f1a 100644 --- a/source/cluster/wham/src-M/probabl.F +++ b/source/cluster/wham/src-M/probabl.F @@ -150,13 +150,13 @@ c write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) 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 diff --git a/source/cluster/wham/src-M/read_coords.F b/source/cluster/wham/src-M/read_coords.F index 56962d9..2911e60 100644 --- a/source/cluster/wham/src-M/read_coords.F +++ b/source/cluster/wham/src-M/read_coords.F @@ -393,7 +393,8 @@ c------------------------------------------------------------------------------ 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 @@ -418,8 +419,8 @@ c------------------------------------------------------------------------------ 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:" diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F index 6f8f5e4..a8c8d1f 100644 --- a/source/cluster/wham/src-M/readrtns.F +++ b/source/cluster/wham/src-M/readrtns.F @@ -37,6 +37,23 @@ C Reading the dimensions of box in x,y,z coordinates 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 @@ -176,6 +193,7 @@ C call reada(weightcard,'WSC',wsc,1.0d0) 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) @@ -780,6 +798,8 @@ C Get parameter filenames and open the parameter files. 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 diff --git a/source/cluster/wham/src-M/ssMD.F b/source/cluster/wham/src-M/ssMD.F index 9fe7d17..42983f8 100644 --- a/source/cluster/wham/src-M/ssMD.F +++ b/source/cluster/wham/src-M/ssMD.F @@ -150,11 +150,117 @@ c-------END TESTING CODE 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) @@ -187,9 +293,9 @@ c om12=dxi*dxj+dyi*dyj+dzi*dzj 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 @@ -223,7 +329,7 @@ c-------TESTING CODE 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 @@ -248,8 +354,8 @@ c-------END TESTING CODE 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 @@ -316,8 +422,8 @@ C write(iout,*) eij,'TU?3' 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) 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 3ad8f52..aafba26 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -3463,7 +3463,7 @@ c end if 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) diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index 730500e..fb18913 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -224,11 +224,11 @@ c detecting NaNQ #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. -- 1.7.9.5