X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new-sep_barrier.F;h=3815c132dcd8f1e18040d13f6096c33382f8a7ac;hb=3054e5f25dc4adf1fbe8476f0fb3677ccaa81129;hp=7d2d27f1f011e11a17b176ba0cd64a2eb53553f0;hpb=57a36a1183cb12334e9bc9ceacabbf038729d459;p=unres.git diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F index 7d2d27f..3815c13 100644 --- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F +++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F @@ -1,4 +1,33 @@ 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----------------------------------------------------------------------- double precision function sscale(r) double precision r,gamm include "COMMON.SPLITELE" @@ -75,8 +104,8 @@ cd & 'iend=',iend(i,iint) rrij=1.0D0/rij 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 evdw=evdw+(1.0d0-sss)*evdwij C @@ -164,8 +193,8 @@ C Change 12/1/95 to calculate four-body interactions rrij=1.0D0/rij 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 evdw=evdw+sss*evdwij C @@ -248,8 +277,8 @@ C if (sss.lt.1.0d0) then 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 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj) @@ -331,8 +360,8 @@ C if (sss.gt.0.0d0) then 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 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj) @@ -445,16 +474,16 @@ 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 evdw=evdw+evdwij*(1.0d0-sss) 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, @@ -558,16 +587,16 @@ 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 evdw=evdw+evdwij*sss 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, @@ -672,6 +701,33 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) 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- + & ((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 @@ -732,8 +788,8 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 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 @@ -742,8 +798,8 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij*(1.0d0-sss) if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, & epsi,sigm,chi1,chi2,chip1,chip2, @@ -766,6 +822,8 @@ C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=ssgradlipi*evdwij + gg_lipj(3)=ssgradlipj*evdwij C Calculate angular part of the gradient. call sc_grad_scale(1.0d0-sss) endif @@ -854,6 +912,32 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) 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- + & ((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 @@ -914,8 +998,8 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 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 @@ -924,8 +1008,8 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij*sss if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, & epsi,sigm,chi1,chi2,chip1,chip2, @@ -948,6 +1032,8 @@ C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=ssgradlipi*evdwij + gg_lipj(3)=ssgradlipj*evdwij C Calculate angular part of the gradient. call sc_grad_scale(sss) endif @@ -1044,8 +1130,8 @@ 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 @@ -1054,8 +1140,8 @@ c--------------------------------------------------------------- evdwij=evdwij*eps2rt*eps3rt evdw=evdw+(evdwij+e_augm)*(1.0d0-sss) if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0), @@ -1166,8 +1252,8 @@ 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 @@ -1176,8 +1262,8 @@ c--------------------------------------------------------------- evdwij=evdwij*eps2rt*eps3rt evdw=evdw+(evdwij+e_augm)*sss if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0), @@ -1234,10 +1320,10 @@ c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12 enddo c write (iout,*) "gg",(gg(k),k=1,3) 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*scalfac - gvdwx(k,j)=gvdwx(k,j)+gg(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) @@ -1249,8 +1335,8 @@ C C Calculate the components of the gradient in DC and X C do l=1,3 - gvdwc(l,i)=gvdwc(l,i)-gg(l) - gvdwc(l,j)=gvdwc(l,j)+gg(l) + gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l) + gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l) enddo return end @@ -1282,6 +1368,7 @@ C include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' + 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), @@ -1368,8 +1455,8 @@ C do i=iturn3_start,iturn3_end if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 - & .or. itype(i+4).eq.ntyp1 +C & .or. itype(i-1).eq.ntyp1 +C & .or. itype(i+4).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -1395,8 +1482,8 @@ C if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & .or. itype(i+3).eq.ntyp1 & .or. itype(i+4).eq.ntyp1 - & .or. itype(i+5).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +C & .or. itype(i+5).eq.ntyp1 +C & .or. itype(i-1).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -1424,8 +1511,8 @@ c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c do i=iatel_s,iatel_e 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 +C & .or. itype(i-1).eq.ntyp1 &) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -1446,8 +1533,8 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) 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 call eelecij_scale(i,j,ees,evdw1,eel_loc) enddo ! j @@ -1483,6 +1570,7 @@ C------------------------------------------------------------------------------- include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' + 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), @@ -1579,8 +1667,14 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions fac3=ael6i*r6ij fac4=ael3i*r3ij evdwij=ev1+ev2 + if (shield_mode.eq.0) then + fac_shield(i)=1.0 + fac_shield(j)=1.0 + endif el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)) el2=fac4*fac + el1=el1*fac_shield(i)**2*fac_shield(j)**2 + el2=el2*fac_shield(i)**2*fac_shield(j)**2 eesij=el1+el2 C 12/26/95 - for the evaluation of multi-body H-bonding interactions ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) @@ -1612,6 +1706,60 @@ 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 + 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 + c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf @@ -1710,7 +1858,9 @@ C ggg(3)=facvdw*zj 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) + ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* + & fac_shield(i)**2*fac_shield(j)**2 + enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -1728,11 +1878,14 @@ cgrad enddo cgrad enddo do k=1,3 gelc(k,i)=gelc(k,i) - & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) - & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) + & +((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) - & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) - & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) + & +((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 gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo @@ -1928,19 +2081,80 @@ cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij + + 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 + + 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 + C Partial derivatives in virtual-bond dihedral angles gamma 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) + & (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) + & (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) + 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) + 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) + gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l) gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) cgrad ghalf=0.5d0*ggg(l) @@ -1954,14 +2168,22 @@ cgrad enddo cgrad enddo 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) - 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) - 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) - 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) + 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 C Change 12/26/95 to calculate four-body contributions to H-bonding energy @@ -2040,8 +2262,19 @@ c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2) ees0mij=0 endif c ees0mij=0.0D0 + 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 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 @@ -2109,17 +2342,29 @@ cgrad ghalfm=0.5D0*gggm(k) 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 ! wcorr endif ! num_conti.le.maxconts