X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fsrc-HCD-5D%2Fenergy_p_new-sep_barrier.F;h=96f77770fdf54e59677d209122ae088cca018393;hb=de4daf245936d67b5c08eb584fca985adab9928e;hp=c6c68321b27f0d14f2da7536bee354d9e89c6d08;hpb=cf4da5a9f21b849afd840360ba40787a35583c71;p=unres.git diff --git a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F index c6c6832..96f7777 100644 --- a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F @@ -85,6 +85,7 @@ c include 'COMMON.CONTACTS' double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & sigij,r0ij,rcut,sss1,sssgrad1,sqrij double precision sscale,sscagrad + double precision boxshift c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 c do i=iatsc_s,iatsc_e @@ -97,6 +98,7 @@ c do i=iatsc_s,iatsc_e xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + call to_box(xi,yi,zi) C C Calculate SC interaction energy. C @@ -106,9 +108,13 @@ cd & 'iend=',iend(i,iint) c do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle - 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) + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) rij=xj*xj+yj*yj+zj*zj sqrij=dsqrt(rrij) eps0ij=eps(itypi,itypj) @@ -187,6 +193,7 @@ c include 'COMMON.CONTACTS' double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & sigij,r0ij,rcut,sqrij,sss1,sssgrad1 double precision sscale,sscagrad + double precision boxshift c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 c do i=iatsc_s,iatsc_e @@ -199,6 +206,7 @@ c do i=iatsc_s,iatsc_e xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + call to_box(xi,yi,zi) C Change 12/1/95 num_conti=0 C @@ -210,9 +218,13 @@ cd & 'iend=',iend(i,iint) c do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle - 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) + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) C Change 12/1/95 to calculate four-body interactions rij=xj*xj+yj*yj+zj*zj sqrij=dsqrt(rij) @@ -285,6 +297,7 @@ C & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck double precision sscale,sscagrad + double precision boxshift c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 c do i=iatsc_s,iatsc_e @@ -297,6 +310,7 @@ c do i=iatsc_s,iatsc_e xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + call to_box(xi,yi,zi) C C Calculate SC interaction energy. C @@ -304,9 +318,13 @@ c do iint=1,nint_gr(i) c do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle - 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) + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) fac_augm=rrij**expon e_augm=augm(itypi,itypj)*fac_augm @@ -384,6 +402,7 @@ C & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck double precision sscale,sscagrad + double precision boxshift c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 c do i=iatsc_s,iatsc_e @@ -396,6 +415,7 @@ c do i=iatsc_s,iatsc_e xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + call to_box(xi,yi,zi) C C Calculate SC interaction energy. C @@ -403,9 +423,13 @@ c do iint=1,nint_gr(i) c do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle - 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) + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) fac_augm=rrij**expon e_augm=augm(itypi,itypj)*fac_augm @@ -478,6 +502,7 @@ C double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi double precision sss1,sssgrad1 double precision sscale,sscagrad + double precision boxshift c double precision rrsave(maxdim) logical lprn evdw=0.0D0 @@ -499,6 +524,7 @@ c do i=iatsc_s,iatsc_e xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + call to_box(xi,yi,zi) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -523,9 +549,13 @@ c dscj_inv=dsc_inv(itypj) 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) + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -604,6 +634,7 @@ C integer itypi,itypj,itypi1,iint,ind,ikont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi double precision sscale,sscagrad + double precision boxshift c double precision rrsave(maxdim) logical lprn evdw=0.0D0 @@ -625,6 +656,7 @@ c do i=iatsc_s,iatsc_e xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + call to_box(xi,yi,zi) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -649,9 +681,13 @@ c dscj_inv=dsc_inv(itypj) 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) + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -731,6 +767,7 @@ C & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip double precision dist,sscale,sscagrad,sscagradlip,sscalelip double precision subchap,sss1,sssgrad1 + double precision boxshift evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -748,12 +785,8 @@ c do i=iatsc_s,iatsc_e 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 + call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -787,79 +820,27 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) 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-((zi-bordlipbot)/lipbufthick) -C lipbufthick is thickenes of lipid buffore - sslipj=sscalelip(fracinbuf) - ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick - else if (zi.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-zi)/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) - sss1=sscale(1.0d0/rij,r_cut_int) - if (sss1.eq.0.0d0) cycle - sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) - if (sss.lt.1.0d0) then + call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) + 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) + sss1=sscale(1.0d0/rij,r_cut_int) + if (sss1.eq.0.0d0) cycle + sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) + if (sss.lt.1.0d0) then C Calculate angle-dependent terms of energy and contributions to their C derivatives. - sssgrad= + sssgrad= & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa) sssgrad1=sscagrad(1.0d0/rij,r_cut_int) call sc_angular @@ -946,15 +927,13 @@ C include 'COMMON.CONTROL' include "COMMON.SPLITELE" logical lprn - integer xshift,yshift,zshift double precision evdw integer itypi,itypj,itypi1,iint,ind,ikont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, - & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, - & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip + & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip double precision dist,sscale,sscagrad,sscagradlip,sscalelip - double precision subchap + double precision boxshift evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -972,12 +951,8 @@ c do i=iatsc_s,iatsc_e 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 + call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1011,67 +986,15 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) 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-((zi-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-zi)/lipbufthick) - sslipj=sscalelip(fracinbuf) - ssgradlipj=sscagradlip(fracinbuf)/lipbufthick - else - sslipj=1.0d0 - ssgradlipj=0.0 - endif - else - sslipj=0.0d0 - ssgradlipj=0.0 - endif + call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) 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 + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -1191,6 +1114,8 @@ c do i=iatsc_s,iatsc_e xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1217,9 +1142,18 @@ c dscj_inv=dsc_inv(itypj) 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) + call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -1317,6 +1251,7 @@ C & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip double precision dist,sscale,sscagrad,sscagradlip,sscalelip + double precision boxshift evdw=0.0D0 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 @@ -1336,6 +1271,8 @@ c do i=iatsc_s,iatsc_e dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) + call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) c dsci_inv=dsc_inv(itypi) dsci_inv=vbld_inv(i+nres) C @@ -1359,9 +1296,18 @@ c dscj_inv=dsc_inv(itypj) 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) + call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -1612,12 +1558,7 @@ C & .or. itype(i+4).eq.ntyp1 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 + call to_box(xmedi,ymedi,zmedi) num_conti=0 call eelecij_scale(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -1641,12 +1582,7 @@ C & .or. itype(i-1).eq.ntyp1 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 + call to_box(xmedi,ymedi,zmedi) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -1677,12 +1613,7 @@ C & .or. itype(i-1).eq.ntyp1 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 + call to_box(xmedi,ymedi,zmedi) c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend #ifdef FOURBODY num_conti=num_cont_hb(i) @@ -1764,6 +1695,7 @@ C------------------------------------------------------------------------------- double precision sss1,sssgrad1 double precision sscale,sscagrad double precision scalar + double precision boxshift c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT @@ -1796,44 +1728,10 @@ C print *,"WCHODZE2" 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 - + call to_box(xj,yj,zj) + xj=boxshift(xj-xmedi,boxxsize) + yj=boxshift(yj-ymedi,boxysize) + zj=boxshift(zj-zmedi,boxzsize) rij=xj*xj+yj*yj+zj*zj rrmij=1.0D0/rij rij=dsqrt(rij) @@ -2734,6 +2632,7 @@ c write (iout,*) "evdwpp_short" double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp, & dist_temp, dist_init,sss_grad double precision sscale,sscagrad + double precision boxshift integer ikont evdw1=0.0D0 C print *,"WCHODZE" @@ -2754,12 +2653,7 @@ c do i=iatel_s_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 + call to_box(xmedi,ymedi,zmedi) num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i), c & ' ielend',ielend_vdw(i) @@ -2781,43 +2675,10 @@ c do j=ielstart_vdw(i),ielend_vdw(i) 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 + call to_box(xj,yj,zj) + xj=boxshift(xj-xmedi,boxxsize) + yj=boxshift(yj-ymedi,boxysize) + zj=boxshift(zj-zmedi,boxzsize) rij=xj*xj+yj*yj+zj*zj rrmij=1.0D0/rij rij=dsqrt(rij) @@ -2881,14 +2742,12 @@ C logical lprint_short common /shortcheck/ lprint_short double precision ggg(3) - integer xshift,yshift,zshift integer i,iint,j,k,iteli,itypj,subchap double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1, & fac,e1,e2,rij double precision evdw2,evdw2_14,evdwij - double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp, - & dist_temp, dist_init double precision sscale,sscagrad + double precision boxshift integer ikont if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb evdw2=0.0D0 @@ -2907,12 +2766,7 @@ c do i=iatscp_s,iatscp_e 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 + call to_box(xi,yi,zi) c do iint=1,nscp_gr(i) @@ -2928,44 +2782,10 @@ C Uncomment following three lines for Ca-p interactions 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 + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) @@ -3063,6 +2883,7 @@ C & dist_temp, dist_init double precision ggg(3) double precision sscale,sscagrad + double precision boxshift evdw2=0.0D0 evdw2_14=0.0d0 cd print '(a)','Enter ESCP' @@ -3076,12 +2897,7 @@ c & ' iatscp_e=',iatscp_e 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 + call to_box(xi,yi,zi) c if (lprint_short) c & write (iout,*) "i",i," itype",itype(i),itype(i+1), @@ -3102,49 +2918,13 @@ C Uncomment following three lines for Ca-p interactions 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 + xj=c(1,j) + yj=c(2,j) + zj=c(3,j) + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) c sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))