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=1f00b2be0cab7209eb77c47484020dfa2f69ca5a;hb=34d3ad3987785642be58fb2f26557d3314215577;hp=815ca5ad7a55c1de3a2223f453a8d70523f47714;hpb=f690e8b70bab14132839afebf080d4a28363b226;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 815ca5a..1f00b2b 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" @@ -13,6 +42,21 @@ C----------------------------------------------------------------------- return end C----------------------------------------------------------------------- +C----------------------------------------------------------------------- + double precision function sscagrad(r) + double precision r,gamm + include "COMMON.SPLITELE" + if(r.lt.r_cut-rlamb) then + sscagrad=0.0d0 + else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then + gamm=(r-(r_cut-rlamb))/rlamb + sscagrad=gamm*(6*gamm-6.0d0)/rlamb + else + sscagrad=0.0d0 + endif + return + end +C----------------------------------------------------------------------- subroutine elj_long(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains @@ -37,7 +81,7 @@ c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e itypi=itype(i) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) @@ -50,7 +94,7 @@ cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) itypj=itype(j) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -60,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 @@ -123,7 +167,7 @@ c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e itypi=itype(i) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) @@ -138,7 +182,7 @@ cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) itypj=itype(j) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -149,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 @@ -209,7 +253,7 @@ c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e itypi=itype(i) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) @@ -220,7 +264,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) itypj=itype(j) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -233,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) @@ -292,7 +336,7 @@ c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e itypi=itype(i) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) @@ -303,7 +347,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) itypj=itype(j) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -316,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) @@ -384,7 +428,7 @@ c endif ind=0 do i=iatsc_s,iatsc_e itypi=itype(i) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) @@ -401,7 +445,7 @@ C do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=itype(j) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) chi1=chi(itypi,itypj) @@ -430,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, @@ -497,7 +541,7 @@ c endif ind=0 do i=iatsc_s,iatsc_e itypi=itype(i) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) @@ -514,7 +558,7 @@ C do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=itype(j) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) chi1=chi(itypi,itypj) @@ -543,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, @@ -598,6 +642,7 @@ C include 'COMMON.CALC' include 'COMMON.CONTROL' logical lprn + integer xshift,yshift,zshift evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -607,11 +652,17 @@ c if (icall.eq.0) lprn=.false. ind=0 do i=iatsc_s,iatsc_e itypi=itype(i) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + 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 dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -626,7 +677,7 @@ C do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=itype(j) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, @@ -642,16 +693,80 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) alf1=alp(itypi) alf2=alp(itypj) alf12=0.5D0*(alf1+alf2) - 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) + 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- + & ((positi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + sslipj=sscalelip(fracinbuf) + ssgradlipj=sscagradlip(fracinbuf)/lipbufthick + else + sslipj=1.0d0 + ssgradlipj=0.0 + endif + else + sslipj=0.0d0 + ssgradlipj=0.0 + endif + aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 + & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + 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) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj))) - + sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj)) if (sss.lt.1.0d0) then C Calculate angle-dependent terms of energy and contributions to their @@ -674,8 +789,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 @@ -684,8 +799,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, @@ -702,11 +817,14 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac + fac=fac+evdwij/(1.0-sss)*(-sssgrad)/sigmaii(itypi,itypj)*rij c fac=0.0d0 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 @@ -736,6 +854,7 @@ C include 'COMMON.CALC' include 'COMMON.CONTROL' logical lprn + integer xshift,yshift,zshift evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -745,11 +864,17 @@ c if (icall.eq.0) lprn=.false. ind=0 do i=iatsc_s,iatsc_e itypi=itype(i) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + 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 dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -764,7 +889,7 @@ C do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=itype(j) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, @@ -780,16 +905,79 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) alf1=alp(itypi) alf2=alp(itypj) alf12=0.5D0*(alf1+alf2) - 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) + 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- + & ((positi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + sslipj=sscalelip(fracinbuf) + ssgradlipj=sscagradlip(fracinbuf)/lipbufthick + else + sslipj=1.0d0 + ssgradlipj=0.0 + endif + else + sslipj=0.0d0 + ssgradlipj=0.0 + endif + aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 + & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + 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) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj))) - + sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj)) if (sss.gt.0.0d0) then C Calculate angle-dependent terms of energy and contributions to their @@ -812,8 +1000,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 @@ -822,8 +1010,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, @@ -840,11 +1028,14 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac + fac=fac+evdwij/sss*sssgrad/sigmaii(itypi,itypj)*rij c fac=0.0d0 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 @@ -882,7 +1073,7 @@ c if (icall.eq.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e itypi=itype(i) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) @@ -899,7 +1090,7 @@ C do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=itype(j) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) @@ -941,8 +1132,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 @@ -951,8 +1142,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), @@ -1004,7 +1195,7 @@ c if (icall.eq.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e itypi=itype(i) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) @@ -1021,7 +1212,7 @@ C do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=itype(j) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) @@ -1063,8 +1254,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 @@ -1073,8 +1264,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), @@ -1131,10 +1322,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)) @@ -1146,8 +1337,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 @@ -1179,6 +1370,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), @@ -1208,6 +1400,7 @@ cd write(iout,*) 'EE',EE(:,:,i) cd enddo cd call check_vecgrad cd stop +C print *,"WCHODZE3" if (icheckgrad.eq.1) then do i=1,nres-1 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i))) @@ -1262,8 +1455,11 @@ C C Loop over i,i+2 and i,i+3 pairs of the peptide groups C do i=iturn3_start,iturn3_end - if (itype(i).eq.21 .or. itype(i+1).eq.21 - & .or. itype(i+2).eq.21 .or. itype(i+3).eq.21) cycle + if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 + & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).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) dzi=dc(3,i) @@ -1273,15 +1469,24 @@ C xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=0 call eelecij_scale(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) num_cont_hb(i)=num_conti enddo do i=iturn4_start,iturn4_end - if (itype(i).eq.21 .or. itype(i+1).eq.21 - & .or. itype(i+3).eq.21 - & .or. itype(i+4).eq.21) cycle + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 + & .or. itype(i+3).eq.ntyp1 + & .or. itype(i+4).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) dzi=dc(3,i) @@ -1291,9 +1496,15 @@ C xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=num_cont_hb(i) call eelecij_scale(i,i+3,ees,evdw1,eel_loc) - if (wturn4.gt.0.0d0 .and. itype(i+2).ne.21) + if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & call eturn4(i,eello_turn4) num_cont_hb(i)=num_conti enddo ! i @@ -1301,7 +1512,10 @@ c 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.21 .or. itype(i+1).eq.21) cycle + 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 dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -1311,10 +1525,19 @@ c xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize 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.21 .or. itype(j+1).eq.21) cycle + if (itype(j).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 num_cont_hb(i)=num_conti @@ -1349,10 +1572,13 @@ C------------------------------------------------------------------------------- include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' + include 'COMMON.SHIELD' + integer xshift,yshift,zshift 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), - & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4) + & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4), + & gmuij2(4),gmuji2(4) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 @@ -1369,6 +1595,7 @@ C 13-go grudnia roku pamietnego... & 0.0d0,0.0d0,1.0d0/ c time00=MPI_Wtime() cd write (iout,*) "eelecij",i,j +C print *,"WCHODZE2" ind=ind+1 iteli=itel(i) itelj=itel(j) @@ -1383,16 +1610,54 @@ cd write (iout,*) "eelecij",i,j dx_normj=dc_norm(1,j) dy_normj=dc_norm(2,j) dz_normj=dc_norm(3,j) - xj=c(1,j)+0.5D0*dxj-xmedi - yj=c(2,j)+0.5D0*dyj-ymedi - zj=c(3,j)+0.5D0*dzj-zmedi + xj=c(1,j)+0.5D0*dxj + yj=c(2,j)+0.5D0*dyj + zj=c(3,j)+0.5D0*dzj + 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 + dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + isubchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + isubchap=1 + endif + enddo + enddo + enddo + if (isubchap.eq.1) then + xj=xj_temp-xmedi + yj=yj_temp-ymedi + zj=zj_temp-zmedi + else + xj=xj_safe-xmedi + yj=yj_safe-ymedi + zj=zj_safe-zmedi + endif + rij=xj*xj+yj*yj+zj*zj rrmij=1.0D0/rij rij=dsqrt(rij) rmij=1.0D0/rij c For extracting the short-range part of Evdwpp sss=sscale(rij/rpp(iteli,itelj)) - + sssgrad=sscagrad(rij/rpp(iteli,itelj)) r3ij=rrmij*rmij r6ij=r3ij*r3ij cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj @@ -1406,8 +1671,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) @@ -1439,6 +1710,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 @@ -1457,9 +1782,9 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj + ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj) + ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj) + ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj) c do k=1,3 c ghalf=0.5D0*ggg(k) c gvdwpp(k,i)=gvdwpp(k,i)+ghalf @@ -1511,9 +1836,12 @@ cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo c 9/28/08 AL Gradient compotents will be summed only at the end - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj +C ggg(1)=facvdw*xj +C ggg(2)=facvdw*yj +C ggg(3)=facvdw*zj + ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj) + ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj) + ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj) do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) @@ -1534,7 +1862,9 @@ c 9/28/08 AL Gradient compotents will be summed only at the end 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) @@ -1552,11 +1882,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 @@ -1582,6 +1915,14 @@ C do l=1,2 kkk=kkk+1 muij(kkk)=mu(k,i)*mu(l,j) +#ifdef NEWCORR + gmuij1(kkk)=gtb1(k,i+1)*mu(l,j) +c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j) + gmuij2(kkk)=gUb2(k,i)*mu(l,j) + gmuji1(kkk)=mu(k,i)*gtb1(l,j+1) +c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i) + gmuji2(kkk)=mu(k,i)*gUb2(l,j) +#endif enddo enddo cd write (iout,*) 'EELEC: i',i,' j',j @@ -1752,19 +2093,128 @@ 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 -C Partial derivatives in virtual-bond dihedral angles gamma + + 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 + +#ifdef NEWCORR + geel_loc_ij=(a22*gmuij1(1) + & +a23*gmuij1(2) + & +a32*gmuij1(3) + & +a33*gmuij1(4)) + & *fac_shield(i)*fac_shield(j) +c write(iout,*) "derivative over thatai" +c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), +c & a33*gmuij1(4) + gloc(nphi+i,icg)=gloc(nphi+i,icg)+ + & geel_loc_ij*wel_loc +c write(iout,*) "derivative over thatai-1" +c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3), +c & a33*gmuij2(4) + geel_loc_ij= + & a22*gmuij2(1) + & +a23*gmuij2(2) + & +a32*gmuij2(3) + & +a33*gmuij2(4) + gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ + & geel_loc_ij*wel_loc + & *fac_shield(i)*fac_shield(j) + +c Derivative over j residue + geel_loc_ji=a22*gmuji1(1) + & +a23*gmuji1(2) + & +a32*gmuji1(3) + & +a33*gmuji1(4) +c write(iout,*) "derivative over thataj" +c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3), +c & a33*gmuji1(4) + + gloc(nphi+j,icg)=gloc(nphi+j,icg)+ + & geel_loc_ji*wel_loc + & *fac_shield(i)*fac_shield(j) + + geel_loc_ji= + & +a22*gmuji2(1) + & +a23*gmuji2(2) + & +a32*gmuji2(3) + & +a33*gmuji2(4) +c write(iout,*) "derivative over thataj-1" +c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), +c & a33*gmuji2(4) + gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ + & geel_loc_ji*wel_loc + & *fac_shield(i)*fac_shield(j) +#endif +cC 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) @@ -1778,14 +2228,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 @@ -1864,8 +2322,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 @@ -1933,17 +2402,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 @@ -1989,18 +2470,21 @@ C include 'COMMON.VECTORS' include 'COMMON.FFIELD' dimension ggg(3) + integer xshift,yshift,zshift c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ #else double precision scal_el /0.5d0/ #endif +c write (iout,*) "evdwpp_short" evdw1=0.0D0 +C print *,"WCHODZE" c write (iout,*) "iatel_s_vdw",iatel_s_vdw, c & " iatel_e_vdw",iatel_e_vdw - call flush(iout) +c call flush(iout) do i=iatel_s_vdw,iatel_e_vdw - if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle + if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2010,12 +2494,18 @@ c & " iatel_e_vdw",iatel_e_vdw xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i), c & ' ielend',ielend_vdw(i) - call flush(iout) +c call flush(iout) do j=ielstart_vdw(i),ielend_vdw(i) - if (itype(j).eq.21 .or. itype(j+1).eq.21) cycle + if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle ind=ind+1 iteli=itel(i) itelj=itel(j) @@ -2028,13 +2518,53 @@ c & ' ielend',ielend_vdw(i) dx_normj=dc_norm(1,j) dy_normj=dc_norm(2,j) dz_normj=dc_norm(3,j) - xj=c(1,j)+0.5D0*dxj-xmedi - yj=c(2,j)+0.5D0*dyj-ymedi - zj=c(3,j)+0.5D0*dzj-zmedi + xj=c(1,j)+0.5D0*dxj + yj=c(2,j)+0.5D0*dyj + zj=c(3,j)+0.5D0*dzj + 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 + dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + isubchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + isubchap=1 + endif + enddo + enddo + enddo + if (isubchap.eq.1) then + xj=xj_temp-xmedi + yj=yj_temp-ymedi + zj=zj_temp-zmedi + else + xj=xj_safe-xmedi + yj=yj_safe-ymedi + zj=zj_safe-zmedi + endif rij=xj*xj+yj*yj+zj*zj rrmij=1.0D0/rij rij=dsqrt(rij) - sss=sscale(rij/rpp(iteli,itelj)) +c sss=sscale(rij/rpp(iteli,itelj)) +c sssgrad=sscagrad(rij/rpp(iteli,itelj)) + sss=sscale(rij) + sssgrad=sscagrad(rij) if (sss.gt.0.0d0) then rmij=1.0D0/rij r3ij=rrmij*rmij @@ -2048,13 +2578,18 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss endif evdw1=evdw1+evdwij*sss + if (energy_dec) write (iout,'(a10,2i5,0pf7.3)') + & 'evdw1_sum',i,j,evdw1 C C Calculate contributions to the Cartesian gradient. C facvdw=-6*rrmij*(ev1+evdwij)*sss - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj + ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj) + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj) + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj) +C ggg(1)=facvdw*xj +C ggg(2)=facvdw*yj +C ggg(3)=facvdw*zj do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) @@ -2082,37 +2617,90 @@ C include 'COMMON.FFIELD' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' + logical lprint_short + common /shortcheck/ lprint_short dimension ggg(3) + integer xshift,yshift,zshift + if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb evdw2=0.0D0 evdw2_14=0.0d0 -cd print '(a)','Enter ESCP' +CD print '(a)','Enter ESCP KURWA' cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e +c if (lprint_short) +c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s, +c & ' iatscp_e=',iatscp_e do i=iatscp_s,iatscp_e - if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) - + 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 do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) itypj=itype(j) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi c zj=c(3,nres+j)-zi C Uncomment following three lines for Ca-p interactions - xj=c(1,j)-xi - yj=c(2,j)-yi - zj=c(3,j)-zi + xj=c(1,j) + yj=c(2,j) + zj=c(3,j) +c corrected by AL + 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 +c end correction + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif + rrij=1.0D0/(xj*xj+yj*yj+zj*zj) sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) - + sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) + if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij), + & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss if (sss.lt.1.0d0) then - fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) e2=fac*bad(itypj,iteli) @@ -2128,7 +2716,9 @@ C Uncomment following three lines for Ca-p interactions C C Calculate contributions to the gradient in the virtual-bond and SC vectors. C + fac=-(evdwij+e1)*rrij*(1.0d0-sss) + fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli) ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac @@ -2183,35 +2773,101 @@ C include 'COMMON.FFIELD' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' + integer xshift,yshift,zshift + logical lprint_short + common /shortcheck/ lprint_short dimension ggg(3) evdw2=0.0D0 evdw2_14=0.0d0 cd print '(a)','Enter ESCP' -cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e +c if (lprint_short) +c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s, +c & ' iatscp_e=',iatscp_e + if (energy_dec) write (iout,*) "escp_short:",r_cut,rlamb do i=iatscp_s,iatscp_e - if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) + 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 +c if (lprint_short) +c & write (iout,*) "i",i," itype",itype(i),itype(i+1), +c & " nscp_gr",nscp_gr(i) do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) itypj=itype(j) - if (itypj.eq.21) cycle +c if (lprint_short) +c & write (iout,*) "j",j," itypj",itypj + if (itypj.eq.ntyp1) cycle C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi c zj=c(3,nres+j)-zi C Uncomment following three lines for Ca-p interactions - xj=c(1,j)-xi - yj=c(2,j)-yi - zj=c(3,j)-zi + xj=c(1,j) + yj=c(2,j) + zj=c(3,j) +c corrected by AL + 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 +c end correction + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 +c if (lprint_short) then +c write (iout,*) i,j,xi,yi,zi,xj,yj,zj +c write (iout,*) "dist_init",dsqrt(dist_init) +c endif + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo +c if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp) + 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 rrij=1.0D0/(xj*xj+yj*yj+zj*zj) - - sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) - +c sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) +c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) + sss=sscale(1.0d0/(dsqrt(rrij))) + sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) + if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij), + & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss +c if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij), +c & " subchap",subchap," sss",sss if (sss.gt.0.0d0) then fac=rrij**expon2 @@ -2230,6 +2886,7 @@ C C Calculate contributions to the gradient in the virtual-bond and SC vectors. C fac=-(evdwij+e1)*rrij*sss + fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli) ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac