X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=dfdaba61f0936aee0e58058c28c9e7626406a17f;hb=22ea58d3be4c7c625b9c5fb0e9f6d1b2bf7428f6;hp=f640679cb349f97b1caa615629310829630ff7e9;hpb=d101c97dea752458d76055fdbae49c26fff03c1f;p=unres.git diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F index f640679..dfdaba6 100644 --- a/source/cluster/wham/src-M/energy_p_new.F +++ b/source/cluster/wham/src-M/energy_p_new.F @@ -22,6 +22,8 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.INTERACT' include 'COMMON.SBRIDGE' include 'COMMON.CHAIN' + include 'COMMON.SHIELD' + include 'COMMON.CONTROL' double precision fact(6) cd write(iout, '(a,i2)')'Calling etotal ipot=',ipot cd print *,'nnt=',nnt,' nct=',nct @@ -47,7 +49,14 @@ C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence). C C Calculate electrostatic (H-bonding) energy of the main chain. C - 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) + 106 continue +C write(iout,*) "shield_mode",shield_mode,ethetacnstr + if (shield_mode.eq.1) then + call set_shield_fac + else if (shield_mode.eq.2) then + call set_shield_fac2 + endif + call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) C C Calculate excluded-volume interaction energy between peptide groups C and side chains. @@ -87,6 +96,23 @@ 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 + + if (TUBElog.eq.1) then + print *,"just before call" + call calctube(Etube) + print *,"just after call",etube + elseif (TUBElog.eq.2) then + call calctube2(Etube) + elseif (TUBElog.eq.3) then + call calcnano(Etube) + else + Etube=0.0d0 + endif + write(iout,*), "Etube",etube C C 12/1/95 Multi-body terms C @@ -98,12 +124,29 @@ c print *,"calling multibody_eello" call multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,n_corr1) c write (*,*) 'n_corr=',n_corr,' n_corr1=',n_corr1 c print *,ecorr,ecorr5,ecorr6,eturn6 + else + ecorr=0.0d0 + ecorr5=0.0d0 + ecorr6=0.0d0 + eturn6=0.0d0 endif 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 + & +welec*fact(1)*ees + & +fact(1)*wvdwpp*evdw1 + & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +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+ethetacnstr + & +wliptran*eliptran+wtube*Etube + else etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees & +wvdwpp*evdw1 & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc @@ -112,7 +155,20 @@ 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 + & +wliptran*eliptran+wtube*Etube + endif #else + if (shield_mode.gt.0) then + etot=fact(1)wsc*(evdw+fact(6)*evdw_t)+fact(1)*wscp*evdw2 + & +welec*fact(1)*(ees+evdw1) + & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +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+ethetacnstr + & +wliptran*eliptran+wtube*Etube + else etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2 & +welec*fact(1)*(ees+evdw1) & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc @@ -121,7 +177,10 @@ 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 + & +wliptran*eliptran+wtube*Etube + endif #endif + energia(0)=etot energia(1)=evdw #ifdef SCP14 @@ -155,6 +214,8 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t energia(20)=edihcnstr energia(21)=evdw_t energia(24)=ethetacnstr + energia(22)=eliptran + energia(25)=Etube c detecting NaNQ #ifdef ISNAN #ifdef AIX @@ -181,6 +242,7 @@ C #ifdef SPLITELE do i=1,nct do j=1,3 + if (shield_mode.eq.0) then gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+ & welec*fact(1)*gelc(j,i)+wvdwpp*gvdwpp(j,i)+ & wbond*gradb(j,i)+ @@ -193,14 +255,63 @@ C & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) + & +wtube*gg_tube(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) + & +wtube*gg_tube_SC(j,i) + + else + gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i) + & +fact(1)*wscp*gvdwc_scp(j,i)+ + & welec*fact(1)*gelc(j,i)+fact(1)*wvdwpp*gvdwpp(j,i)+ + & wbond*gradb(j,i)+ + & wstrain*ghpbc(j,i)+ + & wcorr*fact(3)*gradcorr(j,i)+ + & wel_loc*fact(2)*gel_loc(j,i)+ + & wturn3*fact(2)*gcorr3_turn(j,i)+ + & wturn4*fact(3)*gcorr4_turn(j,i)+ + & wcorr5*fact(4)*gradcorr5(j,i)+ + & wcorr6*fact(5)*gradcorr6(j,i)+ + & wturn6*fact(5)*gcorr6_turn(j,i)+ + & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + & +wcorr*gshieldc_ec(j,i) + & +wcorr*gshieldc_loc_ec(j,i) + & +wturn3*gshieldc_t3(j,i) + & +wturn3*gshieldc_loc_t3(j,i) + & +wturn4*gshieldc_t4(j,i) + & +wturn4*gshieldc_loc_t4(j,i) + & +wel_loc*gshieldc_ll(j,i) + & +wel_loc*gshieldc_loc_ll(j,i) + & +wtube*gg_tube(j,i) + + gradx(j,i,icg)=fact(1)*wsc*gvdwx(j,i) + & +fact(1)*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) + & +welec*gshieldx(j,i) + & +wcorr*gshieldx_ec(j,i) + & +wturn3*gshieldx_t3(j,i) + & +wturn4*gshieldx_t4(j,i) + & +wel_loc*gshieldx_ll(j,i) + & +wtube*gg_tube_SC(j,i) + + + endif enddo #else - do i=1,nct + do i=1,nct do j=1,3 + if (shield_mode.eq.0) then gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+ & welec*fact(1)*gelc(j,i)+wstrain*ghpbc(j,i)+ & wbond*gradb(j,i)+ @@ -212,11 +323,40 @@ C & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) + & +wtube*gg_tube(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) - enddo + & +wliptran*gliptranx(j,i) + & +wtube*gg_tube_SC(j,i) + else + gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i)+ + & fact(1)*wscp*gvdwc_scp(j,i)+ + & welec*fact(1)*gelc(j,i)+wstrain*ghpbc(j,i)+ + & wbond*gradb(j,i)+ + & wcorr*fact(3)*gradcorr(j,i)+ + & wel_loc*fact(2)*gel_loc(j,i)+ + & wturn3*fact(2)*gcorr3_turn(j,i)+ + & wturn4*fact(3)*gcorr4_turn(j,i)+ + & wcorr5*fact(4)*gradcorr5(j,i)+ + & wcorr6*fact(5)*gradcorr6(j,i)+ + & wturn6*fact(5)*gcorr6_turn(j,i)+ + & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) + & +wtube*gg_tube(j,i) + + gradx(j,i,icg)=fact(1)*wsc*gvdwx(j,i)+ + & fact(1)*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) + & +wtube*gg_tube_SC(j,i) + endif + enddo #endif enddo @@ -272,6 +412,7 @@ C------------------------------------------------------------------------ edihcnstr=energia(20) estr=energia(18) ethetacnstr=energia(24) + etube=energia(25) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1, & wvdwpp, @@ -280,7 +421,8 @@ 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,ethetacnstr,ebr*nss,etot + & esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss,etube,wtube, + & etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -304,6 +446,7 @@ C------------------------------------------------------------------------ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ + & 'ETUBE=',1pE16.6,' WEIGHT=',1pD16.6,' (energy with nano)'/ & 'ETOT= ',1pE16.6,' (total)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),estr,wbond, @@ -312,7 +455,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,ethetacnstr,ebr*nss,etot + & edihcnstr,ethetacnstr,ebr*nss,etube,wtube,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -335,6 +478,7 @@ C------------------------------------------------------------------------ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ + & 'ETUBE=',1pE16.6,' WEIGHT=',1pD16.6,' (energy with nano)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -402,8 +546,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 @@ -417,7 +561,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 @@ -569,8 +713,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) @@ -580,7 +724,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 @@ -706,23 +850,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, @@ -794,6 +938,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) @@ -868,6 +1034,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 +C write(iout,*) "czy jest 0", bb-bb_lip(itypi,itypj), +C & bb-bb_aq(itypi,itypj) dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj yj_safe=yj @@ -923,13 +1117,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 @@ -940,15 +1134,19 @@ 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) -c write (iout,'(2(a3,i3,2x),17(0pf7.3))') -c & restyp(itypi),i,restyp(itypj),j, -c & epsi,sigm,chi1,chi2,chip1,chip2, -c & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, -c & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, -c & evdwij -c write (iout,*) "pratial sum", evdw,evdw_t + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa +C#define DEBUG +#ifdef DEBUG +C write (iout,'(2(a3,i3,2x),17(0pf7.3))') +C & restyp(itypi),i,restyp(itypj),j, +C & epsi,sigm,chi1,chi2,chip1,chip2, +C & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, +C & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, +C & evdwij + write (iout,*) "pratial sum", evdw,evdw_t,e1,e2,fac,aa +#endif +C#undef DEBUG c endif if (calc_grad) then C Calculate gradient components. @@ -957,6 +1155,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 @@ -1011,6 +1215,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 @@ -1041,9 +1274,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) @@ -1064,15 +1364,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 @@ -1178,10 +1478,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 @@ -1190,9 +1490,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------------------------------------------------------------------------------ @@ -1837,6 +2140,8 @@ C include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' + include 'COMMON.SHIELD' + dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), @@ -1907,15 +2212,15 @@ cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e gcorr_loc(i)=0.0d0 enddo do i=iatel_s,iatel_e - if (i.eq.1) then +C if (i.eq.1) then if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 - & .or. itype(i+2).eq.ntyp1) cycle - else - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 - & .or. itype(i+2).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +C & .or. itype(i+2).eq.ntyp1) cycle +C else +C if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 +C & .or. itype(i+2).eq.ntyp1 +C & .or. itype(i-1).eq.ntyp1 &) cycle - endif +C endif if (itel(i).eq.0) goto 1215 dxi=dc(1,i) dyi=dc(2,i) @@ -1935,16 +2240,16 @@ cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) - if (j.eq.1) then - if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 - & .or.itype(j+2).eq.ntyp1 - &) cycle - else +C if (j.le.1) cycle +C if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 +C & .or.itype(j+2).eq.ntyp1 +C &) cycle +C else if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 - & .or.itype(j+2).eq.ntyp1 - & .or.itype(j-1).eq.ntyp1 +C & .or.itype(j+2).eq.ntyp1 +C & .or.itype(j-1).eq.ntyp1 &) cycle - endif +C endif if (itel(j).eq.0) goto 1216 ind=ind+1 iteli=itel(i) @@ -1975,7 +2280,7 @@ C End diagnostics if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 xj_safe=xj yj_safe=yj zj_safe=zj @@ -1986,7 +2291,7 @@ C End diagnostics 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 + dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 if(dist_temp.lt.dist_init) then dist_init=dist_temp xj_temp=xj @@ -2032,7 +2337,26 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions c write (iout,*) "i",i,iteli," j",j,itelj," eesij",eesij C 12/26/95 - for the evaluation of multi-body H-bonding interactions ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) + if (shield_mode.gt.0) then +C fac_shield(i)=0.4 +C fac_shield(j)=0.6 +C#define DEBUG +#ifdef DEBUG + write(iout,*) "ees_compon",i,j,el1,el2, + & fac_shield(i),fac_shield(j) +#endif +C#undef DEBUG + el1=el1*fac_shield(i)**2*fac_shield(j)**2 + el2=el2*fac_shield(i)**2*fac_shield(j)**2 + eesij=(el1+el2) + ees=ees+eesij + else + fac_shield(i)=1.0 + fac_shield(j)=1.0 + eesij=(el1+el2) ees=ees+eesij + endif +C ees=ees+eesij evdw1=evdw1+evdwij*sss cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, @@ -2055,6 +2379,63 @@ C ggg(1)=facel*xj ggg(2)=facel*yj ggg(3)=facel*zj + + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. + & (shield_mode.gt.0)) then +C print *,i,j + do ilist=1,ishield_list(i) + iresshield=shield_list(ilist,i) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i) + & *2.0 + gshieldx(k,iresshield)=gshieldx(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0 + gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield +C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield) +C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) +C if (iresshield.gt.i) then +C do ishi=i+1,iresshield-1 +C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield +C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) +C +C enddo +C else +C do ishi=iresshield,i +C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield +C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) +C +C enddo +C endif +C enddo +C enddo + enddo + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j) + & *2.0 + gshieldx(k,iresshield)=gshieldx(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0 + gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield + enddo + enddo + + do k=1,3 + gshieldc(k,i)=gshieldc(k,i)+ + & grad_shield(k,i)*eesij/fac_shield(i)*2.0 + gshieldc(k,j)=gshieldc(k,j)+ + & grad_shield(k,j)*eesij/fac_shield(j)*2.0 + gshieldc(k,i-1)=gshieldc(k,i-1)+ + & grad_shield(k,i)*eesij/fac_shield(i)*2.0 + gshieldc(k,j-1)=gshieldc(k,j-1)+ + & grad_shield(k,j)*eesij/fac_shield(j)*2.0 + + enddo + endif + do k=1,3 ghalf=0.5D0*ggg(k) gelc(k,i)=gelc(k,i)+ghalf @@ -2138,15 +2519,19 @@ cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), cd & (dcosg(k),k=1,3) do k=1,3 ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k) + & *fac_shield(i)**2*fac_shield(j)**2 enddo do k=1,3 ghalf=0.5D0*ggg(k) gelc(k,i)=gelc(k,i)+ghalf & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) + & *fac_shield(i)**2*fac_shield(j)**2 + gelc(k,j)=gelc(k,j)+ghalf & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) + & *fac_shield(i)**2*fac_shield(j)**2 enddo do k=i+1,j-1 do l=1,3 @@ -2394,16 +2779,70 @@ C Contribution to the local-electrostatic energy coming from the i-j pair & +a33*muij(4) cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij cd write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) + if (shield_mode.eq.0) then + fac_shield(i)=1.0 + fac_shield(j)=1.0 +C else +C fac_shield(i)=0.4 +C fac_shield(j)=0.6 + endif + eel_loc_ij=eel_loc_ij + & *fac_shield(i)*fac_shield(j) eel_loc=eel_loc+eel_loc_ij C Partial derivatives in virtual-bond dihedral angles gamma if (calc_grad) then + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. + & (shield_mode.gt.0)) then +C print *,i,j + + do ilist=1,ishield_list(i) + iresshield=shield_list(ilist,i) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij + & /fac_shield(i) +C & *2.0 + gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i) + gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1) + & +rlocshield + enddo + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij + & /fac_shield(j) +C & *2.0 + gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j) + gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1) + & +rlocshield + + enddo + enddo + do k=1,3 + gshieldc_ll(k,i)=gshieldc_ll(k,i)+ + & grad_shield(k,i)*eel_loc_ij/fac_shield(i) + gshieldc_ll(k,j)=gshieldc_ll(k,j)+ + & grad_shield(k,j)*eel_loc_ij/fac_shield(j) + gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+ + & grad_shield(k,i)*eel_loc_ij/fac_shield(i) + gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+ + & grad_shield(k,j)*eel_loc_ij/fac_shield(j) + enddo + endif if (i.gt.1) & gel_loc_loc(i-1)=gel_loc_loc(i-1)+ & a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j) + & *fac_shield(i)*fac_shield(j) gel_loc_loc(j-1)=gel_loc_loc(j-1)+ & a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j) + & *fac_shield(i)*fac_shield(j) + cd call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij) cd write(iout,*) 'agg ',agg cd write(iout,*) 'aggi ',aggi @@ -2415,6 +2854,8 @@ C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) do l=1,3 ggg(l)=agg(l,1)*muij(1)+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4) + & *fac_shield(i)*fac_shield(j) + enddo do k=i+2,j2 do l=1,3 @@ -2425,12 +2866,20 @@ C Remaining derivatives of eello do l=1,3 gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+ & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4) + & *fac_shield(i)*fac_shield(j) + gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+ & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4) + & *fac_shield(i)*fac_shield(j) + gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+ & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4) + & *fac_shield(i)*fac_shield(j) + gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+ & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4) + & *fac_shield(i)*fac_shield(j) + enddo endif ENDIF @@ -2536,9 +2985,21 @@ c fac3=dsqrt(-ael6i)/r0ij**3 fac3=dsqrt(-ael6i)*r3ij ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1) ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2) + if (shield_mode.eq.0) then + fac_shield(i)=1.0d0 + fac_shield(j)=1.0d0 + else + ees0plist(num_conti,i)=j +C fac_shield(i)=0.4d0 +C fac_shield(j)=0.6d0 + endif c ees0mij=0.0D0 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) + & *fac_shield(i)*fac_shield(j) + ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) + & *fac_shield(i)*fac_shield(j) + C Diagnostics. Comment out or remove after debugging! c ees0p(num_conti,i)=0.5D0*fac3*ees0pij c ees0m(num_conti,i)=0.5D0*fac3*ees0mij @@ -2603,17 +3064,29 @@ C Derivatives due to the contact function gacontp_hb1(k,num_conti,i)=ghalfp & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) + & *fac_shield(i)*fac_shield(j) + gacontp_hb2(k,num_conti,i)=ghalfp & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) + & *fac_shield(i)*fac_shield(j) + gacontp_hb3(k,num_conti,i)=gggp(k) + & *fac_shield(i)*fac_shield(j) + gacontm_hb1(k,num_conti,i)=ghalfm & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) + & *fac_shield(i)*fac_shield(j) + gacontm_hb2(k,num_conti,i)=ghalfm & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) + & *fac_shield(i)*fac_shield(j) + gacontm_hb3(k,num_conti,i)=gggm(k) + & *fac_shield(i)*fac_shield(j) + enddo endif C Diagnostics. Comment out or remove after debugging! @@ -2659,6 +3132,9 @@ C Third- and fourth-order contributions from turns include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' + include 'COMMON.SHIELD' + include 'COMMON.CONTROL' + dimension ggg(3) double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2), & e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2), @@ -2667,6 +3143,18 @@ C Third- and fourth-order contributions from turns & aggj(3,4),aggj1(3,4),a_temp(2,2) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2 if (j.eq.i+2) then + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 +C changes suggested by Ana to avoid out of bounds +C & .or.((i+5).gt.nres) +C & .or.((i-1).le.0) +C end of changes suggested by Ana + & .or. itype(i+2).eq.ntyp1 + & .or. itype(i+3).eq.ntyp1 +C & .or. itype(i+5).eq.ntyp1 +C & .or. itype(i).eq.ntyp1 +C & .or. itype(i-1).eq.ntyp1 + & ) goto 179 + CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Third-order contributions @@ -2681,22 +3169,80 @@ cd call checkint_turn3(i,a_temp,eello_turn3_num) call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1)) call transpose2(auxmat(1,1),auxmat1(1,1)) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) + if (shield_mode.eq.0) then + fac_shield(i)=1.0 + fac_shield(j)=1.0 +C else +C fac_shield(i)=0.4 +C fac_shield(j)=0.6 + endif eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + cd write (2,*) 'i,',i,' j',j,'eello_turn3', cd & 0.5d0*(pizda(1,1)+pizda(2,2)), cd & ' eello_turn3_num',4*eello_turn3_num if (calc_grad) then +C Derivatives in shield mode + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. + & (shield_mode.gt.0)) then +C print *,i,j + + do ilist=1,ishield_list(i) + iresshield=shield_list(ilist,i) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,i)*eello_t3/fac_shield(i) +C & *2.0 + gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,i)*eello_t3/fac_shield(i) + gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1) + & +rlocshield + enddo + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,j)*eello_t3/fac_shield(j) +C & *2.0 + gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,j)*eello_t3/fac_shield(j) + gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1) + & +rlocshield + + enddo + enddo + + do k=1,3 + gshieldc_t3(k,i)=gshieldc_t3(k,i)+ + & grad_shield(k,i)*eello_t3/fac_shield(i) + gshieldc_t3(k,j)=gshieldc_t3(k,j)+ + & grad_shield(k,j)*eello_t3/fac_shield(j) + gshieldc_t3(k,i-1)=gshieldc_t3(k,i-1)+ + & grad_shield(k,i)*eello_t3/fac_shield(i) + gshieldc_t3(k,j-1)=gshieldc_t3(k,j-1)+ + & grad_shield(k,j)*eello_t3/fac_shield(j) + enddo + endif + C Derivatives in gamma(i) call matmat2(EUgder(1,1,i+1),EUg(1,1,i+2),auxmat2(1,1)) call transpose2(auxmat2(1,1),pizda(1,1)) call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1)) gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + C Derivatives in gamma(i+1) call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1)) call transpose2(auxmat2(1,1),pizda(1,1)) call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1)) gel_loc_turn3(i+1)=gel_loc_turn3(i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + C Cartesian derivatives do l=1,3 a_temp(1,1)=aggi(l,1) @@ -2706,6 +3252,8 @@ C Cartesian derivatives call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i)=gcorr3_turn(l,i) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) a_temp(2,1)=aggi1(l,3) @@ -2713,6 +3261,8 @@ C Cartesian derivatives call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) a_temp(2,1)=aggj(l,3) @@ -2720,6 +3270,8 @@ C Cartesian derivatives call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j)=gcorr3_turn(l,j) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -2727,9 +3279,24 @@ C Cartesian derivatives call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j1)=gcorr3_turn(l,j1) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + enddo endif + 179 continue else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 +C changes suggested by Ana to avoid out of bounds +C & .or.((i+5).gt.nres) +C & .or.((i-1).le.0) +C end of changes suggested by Ana + & .or. itype(i+3).eq.ntyp1 + & .or. itype(i+4).eq.ntyp1 +C & .or. itype(i+5).eq.ntyp1 + & .or. itype(i).eq.ntyp1 +C & .or. itype(i-1).eq.ntyp1 + & ) goto 178 + CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Fourth-order contributions @@ -2757,11 +3324,64 @@ cd call checkint_turn4(i,a_temp,eello_turn4_num) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) + if (shield_mode.eq.0) then + fac_shield(i)=1.0 + fac_shield(j)=1.0 +C else +C fac_shield(i)=0.4 +C fac_shield(j)=0.6 + endif eello_turn4=eello_turn4-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + eello_t4=-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), cd & ' eello_turn4_num',8*eello_turn4_num C Derivatives in gamma(i) if (calc_grad) then + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. + & (shield_mode.gt.0)) then +C print *,i,j + + do ilist=1,ishield_list(i) + iresshield=shield_list(ilist,i) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,i)*eello_t4/fac_shield(i) +C & *2.0 + gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,i)*eello_t4/fac_shield(i) + gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1) + & +rlocshield + enddo + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,j)*eello_t4/fac_shield(j) +C & *2.0 + gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,j)*eello_t4/fac_shield(j) + gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1) + & +rlocshield + + enddo + enddo + + do k=1,3 + gshieldc_t4(k,i)=gshieldc_t4(k,i)+ + & grad_shield(k,i)*eello_t4/fac_shield(i) + gshieldc_t4(k,j)=gshieldc_t4(k,j)+ + & grad_shield(k,j)*eello_t4/fac_shield(j) + gshieldc_t4(k,i-1)=gshieldc_t4(k,i-1)+ + & grad_shield(k,i)*eello_t4/fac_shield(i) + gshieldc_t4(k,j-1)=gshieldc_t4(k,j-1)+ + & grad_shield(k,j)*eello_t4/fac_shield(j) + enddo + endif + call transpose2(EUgder(1,1,i+1),e1tder(1,1)) call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1)) call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1)) @@ -2769,6 +3389,8 @@ C Derivatives in gamma(i) call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) + & *fac_shield(i)*fac_shield(j) + C Derivatives in gamma(i+1) call transpose2(EUgder(1,1,i+2),e2tder(1,1)) call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) @@ -2777,6 +3399,8 @@ C Derivatives in gamma(i+1) call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3) + & *fac_shield(i)*fac_shield(j) + C Derivatives in gamma(i+2) call transpose2(EUgder(1,1,i+3),e3tder(1,1)) call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1)) @@ -2788,6 +3412,8 @@ C Derivatives in gamma(i+2) call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + C Cartesian derivatives C Derivatives of this turn contributions in DC(i+2) if (j.lt.nres-1) then @@ -2807,6 +3433,8 @@ C Derivatives of this turn contributions in DC(i+2) s3=0.5d0*(pizda(1,1)+pizda(2,2)) ggg(l)=-(s1+s2+s3) gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + enddo endif C Remaining derivatives of this turn contribution @@ -2825,6 +3453,8 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) a_temp(2,1)=aggi1(l,3) @@ -2839,6 +3469,8 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) a_temp(2,1)=aggj(l,3) @@ -2853,6 +3485,8 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -2867,8 +3501,11 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) + enddo endif + 178 continue endif return end @@ -5549,6 +6186,8 @@ c------------------------------------------------------------------------------ include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.SHIELD' + double precision gx(3),gx1(3) logical lprn lprn=.false. @@ -5606,7 +6245,85 @@ C Calculate multi-body contributions to the gradient. & ekont*(coeffp*ees0pij*gacontp_hb3(ll,kk,k)+ & coeffm*ees0mij*gacontm_hb3(ll,kk,k)) enddo - enddo + enddo + if (shield_mode.gt.0) then + j=ees0plist(jj,i) + l=ees0plist(kk,k) +C print *,i,j,fac_shield(i),fac_shield(j), +C &fac_shield(k),fac_shield(l) + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. + & (fac_shield(k).gt.0).and.(fac_shield(l).gt.0)) then + do ilist=1,ishield_list(i) + iresshield=shield_list(ilist,i) + do m=1,3 + rlocshield=grad_shield_side(m,ilist,i)*ehbcorr/fac_shield(i) +C & *2.0 + gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ + & rlocshield + & +grad_shield_loc(m,ilist,i)*ehbcorr/fac_shield(i) + gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) + &+rlocshield + enddo + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do m=1,3 + rlocshield=grad_shield_side(m,ilist,j)*ehbcorr/fac_shield(j) +C & *2.0 + gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ + & rlocshield + & +grad_shield_loc(m,ilist,j)*ehbcorr/fac_shield(j) + gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) + & +rlocshield + enddo + enddo + do ilist=1,ishield_list(k) + iresshield=shield_list(ilist,k) + do m=1,3 + rlocshield=grad_shield_side(m,ilist,k)*ehbcorr/fac_shield(k) +C & *2.0 + gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ + & rlocshield + & +grad_shield_loc(m,ilist,k)*ehbcorr/fac_shield(k) + gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) + & +rlocshield + enddo + enddo + do ilist=1,ishield_list(l) + iresshield=shield_list(ilist,l) + do m=1,3 + rlocshield=grad_shield_side(m,ilist,l)*ehbcorr/fac_shield(l) +C & *2.0 + gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ + & rlocshield + & +grad_shield_loc(m,ilist,l)*ehbcorr/fac_shield(l) + gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) + & +rlocshield + enddo + enddo +C print *,gshieldx(m,iresshield) + do m=1,3 + gshieldc_ec(m,i)=gshieldc_ec(m,i)+ + & grad_shield(m,i)*ehbcorr/fac_shield(i) + gshieldc_ec(m,j)=gshieldc_ec(m,j)+ + & grad_shield(m,j)*ehbcorr/fac_shield(j) + gshieldc_ec(m,i-1)=gshieldc_ec(m,i-1)+ + & grad_shield(m,i)*ehbcorr/fac_shield(i) + gshieldc_ec(m,j-1)=gshieldc_ec(m,j-1)+ + & grad_shield(m,j)*ehbcorr/fac_shield(j) + + gshieldc_ec(m,k)=gshieldc_ec(m,k)+ + & grad_shield(m,k)*ehbcorr/fac_shield(k) + gshieldc_ec(m,l)=gshieldc_ec(m,l)+ + & grad_shield(m,l)*ehbcorr/fac_shield(l) + gshieldc_ec(m,k-1)=gshieldc_ec(m,k-1)+ + & grad_shield(m,k)*ehbcorr/fac_shield(k) + gshieldc_ec(m,l-1)=gshieldc_ec(m,l-1)+ + & grad_shield(m,l)*ehbcorr/fac_shield(l) + + enddo + endif + endif endif ehbcorr=ekont*ees return @@ -7883,4 +8600,1172 @@ C----------------------------------------------------------------------- return end C----------------------------------------------------------------------- +C first for shielding is setting of function of side-chains + subroutine set_shield_fac2 + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.IOUNITS' + include 'COMMON.SHIELD' + include 'COMMON.INTERACT' +C this is the squar root 77 devided by 81 the epislion in lipid (in protein) + double precision div77_81/0.974996043d0/, + &div4_81/0.2222222222d0/,sh_frac_dist_grad(3) + +C the vector between center of side_chain and peptide group + double precision pep_side(3),long,side_calf(3), + &pept_group(3),costhet_grad(3),cosphi_grad_long(3), + &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3) +C the line belowe needs to be changed for FGPROC>1 + do i=1,nres-1 + if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle + ishield_list(i)=0 +Cif there two consequtive dummy atoms there is no peptide group between them +C the line below has to be changed for FGPROC>1 + VolumeTotal=0.0 + do k=1,nres + if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle + dist_pep_side=0.0 + dist_side_calf=0.0 + do j=1,3 +C first lets set vector conecting the ithe side-chain with kth side-chain + pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0 +C pep_side(j)=2.0d0 +C and vector conecting the side-chain with its proper calfa + side_calf(j)=c(j,k+nres)-c(j,k) +C side_calf(j)=2.0d0 + pept_group(j)=c(j,i)-c(j,i+1) +C lets have their lenght + dist_pep_side=pep_side(j)**2+dist_pep_side + dist_side_calf=dist_side_calf+side_calf(j)**2 + dist_pept_group=dist_pept_group+pept_group(j)**2 + enddo + dist_pep_side=dsqrt(dist_pep_side) + dist_pept_group=dsqrt(dist_pept_group) + dist_side_calf=dsqrt(dist_side_calf) + do j=1,3 + pep_side_norm(j)=pep_side(j)/dist_pep_side + side_calf_norm(j)=dist_side_calf + enddo +C now sscale fraction + sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield +C print *,buff_shield,"buff" +C now sscale + if (sh_frac_dist.le.0.0) cycle +C If we reach here it means that this side chain reaches the shielding sphere +C Lets add him to the list for gradient + ishield_list(i)=ishield_list(i)+1 +C ishield_list is a list of non 0 side-chain that contribute to factor gradient +C this list is essential otherwise problem would be O3 + shield_list(ishield_list(i),i)=k +C Lets have the sscale value + if (sh_frac_dist.gt.1.0) then + scale_fac_dist=1.0d0 + do j=1,3 + sh_frac_dist_grad(j)=0.0d0 + enddo + else + scale_fac_dist=-sh_frac_dist*sh_frac_dist + & *(2.0d0*sh_frac_dist-3.0d0) + fac_help_scale=6.0d0*(sh_frac_dist-sh_frac_dist**2) + & /dist_pep_side/buff_shield*0.5d0 +C remember for the final gradient multiply sh_frac_dist_grad(j) +C for side_chain by factor -2 ! + do j=1,3 + sh_frac_dist_grad(j)=fac_help_scale*pep_side(j) +C sh_frac_dist_grad(j)=0.0d0 +C scale_fac_dist=1.0d0 +C print *,"jestem",scale_fac_dist,fac_help_scale, +C & sh_frac_dist_grad(j) + enddo + endif +C this is what is now we have the distance scaling now volume... + short=short_r_sidechain(itype(k)) + long=long_r_sidechain(itype(k)) + costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2) + sinthet=short/dist_pep_side*costhet +C now costhet_grad +C costhet=0.6d0 +C sinthet=0.8 + costhet_fac=costhet**3*short**2*(-0.5d0)/dist_pep_side**4 +C sinthet_fac=costhet**2*0.5d0*(short**3/dist_pep_side**4*costhet +C & -short/dist_pep_side**2/costhet) +C costhet_fac=0.0d0 + do j=1,3 + costhet_grad(j)=costhet_fac*pep_side(j) + enddo +C remember for the final gradient multiply costhet_grad(j) +C for side_chain by factor -2 ! +C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1 +C pep_side0pept_group is vector multiplication + pep_side0pept_group=0.0d0 + do j=1,3 + pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j) + enddo + cosalfa=(pep_side0pept_group/ + & (dist_pep_side*dist_side_calf)) + fac_alfa_sin=1.0d0-cosalfa**2 + fac_alfa_sin=dsqrt(fac_alfa_sin) + rkprim=fac_alfa_sin*(long-short)+short +C rkprim=short + +C now costhet_grad + cosphi=1.0d0/dsqrt(1.0d0+rkprim**2/dist_pep_side**2) +C cosphi=0.6 + cosphi_fac=cosphi**3*rkprim**2*(-0.5d0)/dist_pep_side**4 + sinphi=rkprim/dist_pep_side/dsqrt(1.0d0+rkprim**2/ + & dist_pep_side**2) +C sinphi=0.8 + do j=1,3 + cosphi_grad_long(j)=cosphi_fac*pep_side(j) + &+cosphi**3*0.5d0/dist_pep_side**2*(-rkprim) + &*(long-short)/fac_alfa_sin*cosalfa/ + &((dist_pep_side*dist_side_calf))* + &((side_calf(j))-cosalfa* + &((pep_side(j)/dist_pep_side)*dist_side_calf)) +C cosphi_grad_long(j)=0.0d0 + cosphi_grad_loc(j)=cosphi**3*0.5d0/dist_pep_side**2*(-rkprim) + &*(long-short)/fac_alfa_sin*cosalfa + &/((dist_pep_side*dist_side_calf))* + &(pep_side(j)- + &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side) +C cosphi_grad_loc(j)=0.0d0 + enddo +C print *,sinphi,sinthet + VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet)) + & /VSolvSphere_div +C & *wshield +C now the gradient... + do j=1,3 + grad_shield(j,i)=grad_shield(j,i) +C gradient po skalowaniu + & +(sh_frac_dist_grad(j)*VofOverlap +C gradient po costhet + & +scale_fac_dist*VSolvSphere/VSolvSphere_div/4.0d0* + &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*( + & sinphi/sinthet*costhet*costhet_grad(j) + & +sinthet/sinphi*cosphi*cosphi_grad_long(j))) + & )*wshield +C grad_shield_side is Cbeta sidechain gradient + grad_shield_side(j,ishield_list(i),i)= + & (sh_frac_dist_grad(j)*-2.0d0 + & *VofOverlap + & -scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0* + &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*( + & sinphi/sinthet*costhet*costhet_grad(j) + & +sinthet/sinphi*cosphi*cosphi_grad_long(j))) + & )*wshield + + grad_shield_loc(j,ishield_list(i),i)= + & scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0* + &(1.0d0/(dsqrt(1.0d0-sinphi*sinthet))*( + & sinthet/sinphi*cosphi*cosphi_grad_loc(j) + & )) + & *wshield + enddo + VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist + enddo + fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield) +C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) + enddo + return + end +C first for shielding is setting of function of side-chains + subroutine set_shield_fac + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.IOUNITS' + include 'COMMON.SHIELD' + include 'COMMON.INTERACT' +C this is the squar root 77 devided by 81 the epislion in lipid (in protein) + double precision div77_81/0.974996043d0/, + &div4_81/0.2222222222d0/,sh_frac_dist_grad(3) + +C the vector between center of side_chain and peptide group + double precision pep_side(3),long,side_calf(3), + &pept_group(3),costhet_grad(3),cosphi_grad_long(3), + &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3) +C the line belowe needs to be changed for FGPROC>1 + do i=1,nres-1 + if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle + ishield_list(i)=0 +Cif there two consequtive dummy atoms there is no peptide group between them +C the line below has to be changed for FGPROC>1 + VolumeTotal=0.0 + do k=1,nres + if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle + dist_pep_side=0.0 + dist_side_calf=0.0 + do j=1,3 +C first lets set vector conecting the ithe side-chain with kth side-chain + pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0 +C pep_side(j)=2.0d0 +C and vector conecting the side-chain with its proper calfa + side_calf(j)=c(j,k+nres)-c(j,k) +C side_calf(j)=2.0d0 + pept_group(j)=c(j,i)-c(j,i+1) +C lets have their lenght + dist_pep_side=pep_side(j)**2+dist_pep_side + dist_side_calf=dist_side_calf+side_calf(j)**2 + dist_pept_group=dist_pept_group+pept_group(j)**2 + enddo + dist_pep_side=dsqrt(dist_pep_side) + dist_pept_group=dsqrt(dist_pept_group) + dist_side_calf=dsqrt(dist_side_calf) + do j=1,3 + pep_side_norm(j)=pep_side(j)/dist_pep_side + side_calf_norm(j)=dist_side_calf + enddo +C now sscale fraction + sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield +C print *,buff_shield,"buff" +C now sscale + if (sh_frac_dist.le.0.0) cycle +C If we reach here it means that this side chain reaches the shielding sphere +C Lets add him to the list for gradient + ishield_list(i)=ishield_list(i)+1 +C ishield_list is a list of non 0 side-chain that contribute to factor gradient +C this list is essential otherwise problem would be O3 + shield_list(ishield_list(i),i)=k +C Lets have the sscale value + if (sh_frac_dist.gt.1.0) then + scale_fac_dist=1.0d0 + do j=1,3 + sh_frac_dist_grad(j)=0.0d0 + enddo + else + scale_fac_dist=-sh_frac_dist*sh_frac_dist + & *(2.0*sh_frac_dist-3.0d0) + fac_help_scale=6.0*(sh_frac_dist-sh_frac_dist**2) + & /dist_pep_side/buff_shield*0.5 +C remember for the final gradient multiply sh_frac_dist_grad(j) +C for side_chain by factor -2 ! + do j=1,3 + sh_frac_dist_grad(j)=fac_help_scale*pep_side(j) +C print *,"jestem",scale_fac_dist,fac_help_scale, +C & sh_frac_dist_grad(j) + enddo + endif +C if ((i.eq.3).and.(k.eq.2)) then +C print *,i,sh_frac_dist,dist_pep,fac_help_scale,scale_fac_dist +C & ,"TU" +C endif + +C this is what is now we have the distance scaling now volume... + short=short_r_sidechain(itype(k)) + long=long_r_sidechain(itype(k)) + costhet=1.0d0/dsqrt(1.0+short**2/dist_pep_side**2) +C now costhet_grad +C costhet=0.0d0 + costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**4 +C costhet_fac=0.0d0 + do j=1,3 + costhet_grad(j)=costhet_fac*pep_side(j) + enddo +C remember for the final gradient multiply costhet_grad(j) +C for side_chain by factor -2 ! +C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1 +C pep_side0pept_group is vector multiplication + pep_side0pept_group=0.0 + do j=1,3 + pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j) + enddo + cosalfa=(pep_side0pept_group/ + & (dist_pep_side*dist_side_calf)) + fac_alfa_sin=1.0-cosalfa**2 + fac_alfa_sin=dsqrt(fac_alfa_sin) + rkprim=fac_alfa_sin*(long-short)+short +C now costhet_grad + cosphi=1.0d0/dsqrt(1.0+rkprim**2/dist_pep_side**2) + cosphi_fac=cosphi**3*rkprim**2*(-0.5)/dist_pep_side**4 + + do j=1,3 + cosphi_grad_long(j)=cosphi_fac*pep_side(j) + &+cosphi**3*0.5/dist_pep_side**2*(-rkprim) + &*(long-short)/fac_alfa_sin*cosalfa/ + &((dist_pep_side*dist_side_calf))* + &((side_calf(j))-cosalfa* + &((pep_side(j)/dist_pep_side)*dist_side_calf)) + + cosphi_grad_loc(j)=cosphi**3*0.5/dist_pep_side**2*(-rkprim) + &*(long-short)/fac_alfa_sin*cosalfa + &/((dist_pep_side*dist_side_calf))* + &(pep_side(j)- + &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side) + enddo + + VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi) + & /VSolvSphere_div + & *wshield +C now the gradient... +C grad_shield is gradient of Calfa for peptide groups +C write(iout,*) "shield_compon",i,k,VSolvSphere,scale_fac_dist, +C & costhet,cosphi +C write(iout,*) "cosphi_compon",i,k,pep_side0pept_group, +C & dist_pep_side,dist_side_calf,c(1,k+nres),c(1,k),itype(k) + do j=1,3 + grad_shield(j,i)=grad_shield(j,i) +C gradient po skalowaniu + & +(sh_frac_dist_grad(j) +C gradient po costhet + &-scale_fac_dist*costhet_grad(j)/(1.0-costhet) + &-scale_fac_dist*(cosphi_grad_long(j)) + &/(1.0-cosphi) )*div77_81 + &*VofOverlap +C grad_shield_side is Cbeta sidechain gradient + grad_shield_side(j,ishield_list(i),i)= + & (sh_frac_dist_grad(j)*-2.0d0 + & +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet) + & +scale_fac_dist*(cosphi_grad_long(j)) + & *2.0d0/(1.0-cosphi)) + & *div77_81*VofOverlap + + grad_shield_loc(j,ishield_list(i),i)= + & scale_fac_dist*cosphi_grad_loc(j) + & *2.0d0/(1.0-cosphi) + & *div77_81*VofOverlap + enddo + VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist + enddo + fac_shield(i)=VolumeTotal*div77_81+div4_81 +C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) + enddo + 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------------------------------------------------------------------------------------- +C----------------------------------------------------------------------- +C----------------------------------------------------------- +C This subroutine is to mimic the histone like structure but as well can be +C utilizet to nanostructures (infinit) small modification has to be used to +C make it finite (z gradient at the ends has to be changes as well as the x,y +C gradient has to be modified at the ends +C The energy function is Kihara potential +C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6) +C 4eps is depth of well sigma is r_minimum r is distance from center of tube +C and r0 is the excluded size of nanotube (can be set to 0 if we want just a +C simple Kihara potential + subroutine calctube(Etube) + 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' + double precision tub_r,vectube(3),enetube(maxres*2) + Etube=0.0d0 + do i=itube_start,itube_end + enetube(i)=0.0d0 + enetube(i+nres)=0.0d0 + enddo +C first we calculate the distance from tube center +C first sugare-phosphate group for NARES this would be peptide group +C for UNRES + do i=itube_start,itube_end +C lets ommit dummy atoms for now + if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle +C now calculate distance from center of tube and direction vectors + xmin=boxxsize + ymin=boxysize + do j=-1,1 + vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) + vectube(2)=vectube(2)+boxysize*j + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + +C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1) +C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2) + +C as the tube is infinity we do not calculate the Z-vector use of Z +C as chosen axis + vectube(3)=0.0d0 +C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r +C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +C and its 6 power + rdiff6=rdiff**6.0d0 +C for vectorization reasons we will sumup at the end to avoid depenence of previous + enetube(i)=pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6 +C write(iout,*) "TU13",i,rdiff6,enetube(i) +C print *,rdiff,rdiff6,pep_aa_tube +C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +C now we calculate gradient + fac=(-12.0d0*pep_aa_tube/rdiff6- + & 6.0d0*pep_bb_tube)/rdiff6/rdiff +C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i), +C &rdiff,fac + +C now direction of gg_tube vector + do j=1,3 + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0 + gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0 + enddo + enddo +C basically thats all code now we split for side-chains (REMEMBER to sum up at the END) +C print *,gg_tube(1,0),"TU" + + + do i=itube_start,itube_end +C Lets not jump over memory as we use many times iti + iti=itype(i) +C lets ommit dummy atoms for now + if ((iti.eq.ntyp1) +C in UNRES uncomment the line below as GLY has no side-chain... +C .or.(iti.eq.10) + & ) cycle + xmin=boxxsize + ymin=boxysize + do j=-1,1 + vectube(1)=mod((c(1,i+nres)),boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=mod((c(2,i+nres)),boxysize) + vectube(2)=vectube(2)+boxysize*j + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp +C write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2), +C & tubecenter(2) + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + +C as the tube is infinity we do not calculate the Z-vector use of Z +C as chosen axis + vectube(3)=0.0d0 +C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r + +C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +C and its 6 power + rdiff6=rdiff**6.0d0 +C for vectorization reasons we will sumup at the end to avoid depenence of previous + sc_aa_tube=sc_aa_tube_par(iti) + sc_bb_tube=sc_bb_tube_par(iti) + enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6 +C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +C now we calculate gradient + fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff- + & 6.0d0*sc_bb_tube/rdiff6/rdiff +C now direction of gg_tube vector + do j=1,3 + gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac + enddo + enddo + do i=itube_start,itube_end + Etube=Etube+enetube(i)+enetube(i+nres) + enddo +C print *,"ETUBE", etube + return + end +C TO DO 1) add to total energy +C 2) add to gradient summation +C 3) add reading parameters (AND of course oppening of PARAM file) +C 4) add reading the center of tube +C 5) add COMMONs +C 6) add to zerograd + +C----------------------------------------------------------------------- +C----------------------------------------------------------- +C This subroutine is to mimic the histone like structure but as well can be +C utilizet to nanostructures (infinit) small modification has to be used to +C make it finite (z gradient at the ends has to be changes as well as the x,y +C gradient has to be modified at the ends +C The energy function is Kihara potential +C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6) +C 4eps is depth of well sigma is r_minimum r is distance from center of tube +C and r0 is the excluded size of nanotube (can be set to 0 if we want just a +C simple Kihara potential + subroutine calctube2(Etube) + 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' + double precision tub_r,vectube(3),enetube(maxres*2) + Etube=0.0d0 + do i=itube_start,itube_end + enetube(i)=0.0d0 + enetube(i+nres)=0.0d0 + enddo +C first we calculate the distance from tube center +C first sugare-phosphate group for NARES this would be peptide group +C for UNRES + do i=itube_start,itube_end +C lets ommit dummy atoms for now + + if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle +C now calculate distance from center of tube and direction vectors +C vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) +C if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize +C vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) +C if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize + xmin=boxxsize + ymin=boxysize + do j=-1,1 + vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) + vectube(2)=vectube(2)+boxysize*j + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + +C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1) +C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2) + +C as the tube is infinity we do not calculate the Z-vector use of Z +C as chosen axis + vectube(3)=0.0d0 +C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r +C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +C and its 6 power + rdiff6=rdiff**6.0d0 +C THIS FRAGMENT MAKES TUBE FINITE + positi=mod((c(3,i)+c(3,i+1))/2.0d0,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) + print *,positi,bordtubebot,buftubebot,bordtubetop + if ((positi.gt.bordtubebot) + & .and.(positi.lt.bordtubetop)) then +C the energy transfer exist + if (positi.lt.buftubebot) then + fracinbuf=1.0d0- + & ((positi-bordtubebot)/tubebufthick) +C lipbufthick is thickenes of lipid buffore + sstube=sscalelip(fracinbuf) + ssgradtube=-sscagradlip(fracinbuf)/tubebufthick + print *,ssgradtube, sstube,tubetranene(itype(i)) + enetube(i)=enetube(i)+sstube*tubetranenepep +C gg_tube_SC(3,i)=gg_tube_SC(3,i) +C &+ssgradtube*tubetranene(itype(i)) +C gg_tube(3,i-1)= gg_tube(3,i-1) +C &+ssgradtube*tubetranene(itype(i)) +C print *,"doing sccale for lower part" + elseif (positi.gt.buftubetop) then + fracinbuf=1.0d0- + &((bordtubetop-positi)/tubebufthick) + sstube=sscalelip(fracinbuf) + ssgradtube=sscagradlip(fracinbuf)/tubebufthick + enetube(i)=enetube(i)+sstube*tubetranenepep +C gg_tube_SC(3,i)=gg_tube_SC(3,i) +C &+ssgradtube*tubetranene(itype(i)) +C gg_tube(3,i-1)= gg_tube(3,i-1) +C &+ssgradtube*tubetranene(itype(i)) +C print *, "doing sscalefor top part",sslip,fracinbuf + else + sstube=1.0d0 + ssgradtube=0.0d0 + enetube(i)=enetube(i)+sstube*tubetranenepep +C print *,"I am in true lipid" + endif + else +C sstube=0.0d0 +C ssgradtube=0.0d0 + cycle + endif ! if in lipid or buffor + +C for vectorization reasons we will sumup at the end to avoid depenence of previous + enetube(i)=enetube(i)+sstube* + &(pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6) +C write(iout,*) "TU13",i,rdiff6,enetube(i) +C print *,rdiff,rdiff6,pep_aa_tube +C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +C now we calculate gradient + fac=(-12.0d0*pep_aa_tube/rdiff6- + & 6.0d0*pep_bb_tube)/rdiff6/rdiff*sstube +C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i), +C &rdiff,fac + +C now direction of gg_tube vector + do j=1,3 + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0 + gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0 + enddo + gg_tube(3,i)=gg_tube(3,i) + &+ssgradtube*enetube(i)/sstube/2.0d0 + gg_tube(3,i-1)= gg_tube(3,i-1) + &+ssgradtube*enetube(i)/sstube/2.0d0 + + enddo +C basically thats all code now we split for side-chains (REMEMBER to sum up at the END) +C print *,gg_tube(1,0),"TU" + do i=itube_start,itube_end +C Lets not jump over memory as we use many times iti + iti=itype(i) +C lets ommit dummy atoms for now + if ((iti.eq.ntyp1) +C in UNRES uncomment the line below as GLY has no side-chain... + & .or.(iti.eq.10) + & ) cycle + vectube(1)=c(1,i+nres) + vectube(1)=mod(vectube(1),boxxsize) + if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize + vectube(2)=c(2,i+nres) + vectube(2)=mod(vectube(2),boxysize) + if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize + + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) +C THIS FRAGMENT MAKES TUBE FINITE + 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) + print *,positi,bordtubebot,buftubebot,bordtubetop + if ((positi.gt.bordtubebot) + & .and.(positi.lt.bordtubetop)) then +C the energy transfer exist + if (positi.lt.buftubebot) then + fracinbuf=1.0d0- + & ((positi-bordtubebot)/tubebufthick) +C lipbufthick is thickenes of lipid buffore + sstube=sscalelip(fracinbuf) + ssgradtube=-sscagradlip(fracinbuf)/tubebufthick + print *,ssgradtube, sstube,tubetranene(itype(i)) + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) +C gg_tube_SC(3,i)=gg_tube_SC(3,i) +C &+ssgradtube*tubetranene(itype(i)) +C gg_tube(3,i-1)= gg_tube(3,i-1) +C &+ssgradtube*tubetranene(itype(i)) +C print *,"doing sccale for lower part" + elseif (positi.gt.buftubetop) then + fracinbuf=1.0d0- + &((bordtubetop-positi)/tubebufthick) + sstube=sscalelip(fracinbuf) + ssgradtube=sscagradlip(fracinbuf)/tubebufthick + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) +C gg_tube_SC(3,i)=gg_tube_SC(3,i) +C &+ssgradtube*tubetranene(itype(i)) +C gg_tube(3,i-1)= gg_tube(3,i-1) +C &+ssgradtube*tubetranene(itype(i)) +C print *, "doing sscalefor top part",sslip,fracinbuf + else + sstube=1.0d0 + ssgradtube=0.0d0 + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) +C print *,"I am in true lipid" + endif + else +C sstube=0.0d0 +C ssgradtube=0.0d0 + cycle + endif ! if in lipid or buffor +CEND OF FINITE FRAGMENT +C as the tube is infinity we do not calculate the Z-vector use of Z +C as chosen axis + vectube(3)=0.0d0 +C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r +C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +C and its 6 power + rdiff6=rdiff**6.0d0 +C for vectorization reasons we will sumup at the end to avoid depenence of previous + sc_aa_tube=sc_aa_tube_par(iti) + sc_bb_tube=sc_bb_tube_par(iti) + enetube(i+nres)=(sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6) + & *sstube+enetube(i+nres) +C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +C now we calculate gradient + fac=(-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff- + & 6.0d0*sc_bb_tube/rdiff6/rdiff)*sstube +C now direction of gg_tube vector + do j=1,3 + gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac + enddo + gg_tube_SC(3,i)=gg_tube_SC(3,i) + &+ssgradtube*enetube(i+nres)/sstube + gg_tube(3,i-1)= gg_tube(3,i-1) + &+ssgradtube*enetube(i+nres)/sstube + + enddo + do i=itube_start,itube_end + Etube=Etube+enetube(i)+enetube(i+nres) + enddo +C print *,"ETUBE", etube + return + end +C TO DO 1) add to total energy +C 2) add to gradient summation +C 3) add reading parameters (AND of course oppening of PARAM file) +C 4) add reading the center of tube +C 5) add COMMONs +C 6) add to zerograd + + +C#------------------------------------------------------------------------------- +C This subroutine is to mimic the histone like structure but as well can be +C utilizet to nanostructures (infinit) small modification has to be used to +C make it finite (z gradient at the ends has to be changes as well as the x,y +C gradient has to be modified at the ends +C The energy function is Kihara potential +C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6) +C 4eps is depth of well sigma is r_minimum r is distance from center of tube +C and r0 is the excluded size of nanotube (can be set to 0 if we want just a +C simple Kihara potential + subroutine calcnano(Etube) + 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' + double precision tub_r,vectube(3),enetube(maxres*2), + & enecavtube(maxres*2) + Etube=0.0d0 + do i=itube_start,itube_end + enetube(i)=0.0d0 + enetube(i+nres)=0.0d0 + enddo +C first we calculate the distance from tube center +C first sugare-phosphate group for NARES this would be peptide group +C for UNRES + do i=itube_start,itube_end +C lets ommit dummy atoms for now + if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle +C now calculate distance from center of tube and direction vectors + xmin=boxxsize + ymin=boxysize + zmin=boxzsize + + do j=-1,1 + vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) + vectube(2)=vectube(2)+boxysize*j + vectube(3)=mod((c(3,i)+c(3,i+1))/2.0d0,boxzsize) + vectube(3)=vectube(3)+boxzsize*j + + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + zminact=abs(vectube(3)-tubecenter(3)) + + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + if (zmin.gt.zminact) then + zmin=zminact + ztemp=vectube(3) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp + vectube(3)=ztemp + + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + vectube(3)=vectube(3)-tubecenter(3) + +C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1) +C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2) +C as the tube is infinity we do not calculate the Z-vector use of Z +C as chosen axis +C vectube(3)=0.0d0 +C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r + vectube(3)=vectube(3)/tub_r +C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +C and its 6 power + rdiff6=rdiff**6.0d0 +C for vectorization reasons we will sumup at the end to avoid depenence of previous + enetube(i)=pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6 +C write(iout,*) "TU13",i,rdiff6,enetube(i) +C print *,rdiff,rdiff6,pep_aa_tube +C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +C now we calculate gradient + fac=(-12.0d0*pep_aa_tube/rdiff6- + & 6.0d0*pep_bb_tube)/rdiff6/rdiff +C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i), +C &rdiff,fac + if (acavtubpep.eq.0.0d0) then +C go to 667 + enecavtube(i)=0.0 + faccav=0.0 + else + denominator=(1.0+dcavtubpep*rdiff6*rdiff6) + enecavtube(i)= + & (bcavtubpep*rdiff+acavtubpep*sqrt(rdiff)+ccavtubpep) + & /denominator + enecavtube(i)=0.0 + faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/sqrt(rdiff)) + & *denominator-(bcavtubpep*rdiff+acavtubpep*sqrt(rdiff) + & +ccavtubpep)*rdiff6**2.0d0/rdiff*dcavtubpep*12.0d0) + & /denominator**2.0d0 +C faccav=0.0 +C fac=fac+faccav +C 667 continue + endif +C print *,"TUT",i,iti,rdiff,rdiff6,acavtubpep,denominator, +C & enecavtube(i),faccav +C print *,"licz=", +C & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti)) +CX print *,"finene=",enetube(i+nres)+enecavtube(i) + +C now direction of gg_tube vector + do j=1,3 + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0 + gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0 + enddo + enddo + + do i=itube_start,itube_end + enecavtube(i)=0.0 +C Lets not jump over memory as we use many times iti + iti=itype(i) +C lets ommit dummy atoms for now + if ((iti.eq.ntyp1) +C in UNRES uncomment the line below as GLY has no side-chain... +C .or.(iti.eq.10) + & ) cycle + xmin=boxxsize + ymin=boxysize + zmin=boxzsize + do j=-1,1 + vectube(1)=mod((c(1,i+nres)),boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=mod((c(2,i+nres)),boxysize) + vectube(2)=vectube(2)+boxysize*j + vectube(3)=mod((c(3,i+nres)),boxzsize) + vectube(3)=vectube(3)+boxzsize*j + + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + zminact=abs(vectube(3)-tubecenter(3)) + + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + if (zmin.gt.zminact) then + zmin=zminact + ztemp=vectube(3) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp + vectube(3)=ztemp + +C write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2), +C & tubecenter(2) + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + vectube(3)=vectube(3)-tubecenter(3) +C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r + vectube(3)=vectube(3)/tub_r + +C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +C and its 6 power + rdiff6=rdiff**6.0d0 +C for vectorization reasons we will sumup at the end to avoid depenence of previous + sc_aa_tube=sc_aa_tube_par(iti) + sc_bb_tube=sc_bb_tube_par(iti) + enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6 +C enetube(i+nres)=0.0d0 +C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +C now we calculate gradient + fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff- + & 6.0d0*sc_bb_tube/rdiff6/rdiff +C fac=0.0 +C now direction of gg_tube vector +C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12) + if (acavtub(iti).eq.0.0d0) then +C go to 667 + enecavtube(i+nres)=0.0 + faccav=0.0 + else + denominator=(1.0+dcavtub(iti)*rdiff6*rdiff6) + enecavtube(i+nres)= + & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti)) + & /denominator +C enecavtube(i)=0.0 + faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/sqrt(rdiff)) + & *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff) + & +ccavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti)*12.0d0) + & /denominator**2.0d0 +C faccav=0.0 + fac=fac+faccav +C 667 continue + endif +C print *,"TUT",i,iti,rdiff,rdiff6,acavtub(iti),denominator, +C & enecavtube(i),faccav +C print *,"licz=", +C & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti)) +C print *,"finene=",enetube(i+nres)+enecavtube(i) + do j=1,3 + gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac + enddo + enddo +C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12) +C do i=itube_start,itube_end +C enecav(i)=0.0 +C iti=itype(i) +C if (acavtub(iti).eq.0.0) cycle + + + + do i=itube_start,itube_end + Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i) + & +enecavtube(i+nres) + enddo +C print *,"ETUBE", etube + return + end +C TO DO 1) add to total energy +C 2) add to gradient summation +C 3) add reading parameters (AND of course oppening of PARAM file) +C 4) add reading the center of tube +C 5) add COMMONs +C 6) add to zerograd