X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Fenergy_p_new-sep_barrier.F;h=96f77770fdf54e59677d209122ae088cca018393;hb=de4daf245936d67b5c08eb584fca985adab9928e;hp=93fe9ab1ca2800a5adb87741f575004219ee33a9;hpb=5836ecdab5a8b95f079bbf6e07374dee3fce8a26;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 93fe9ab..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 @@ -81,31 +81,40 @@ C c include 'COMMON.CONTACTS' double precision gg(3) double precision evdw,evdwij - integer i,j,k,itypi,itypj,itypi1,num_conti,iint + integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont 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 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) 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 - do iint=1,nint_gr(i) +c do iint=1,nint_gr(i) cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) - do j=istart(i,iint),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) @@ -138,8 +147,8 @@ C gvdwc(k,j)=gvdwc(k,j)+gg(k) enddo endif - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i do i=1,nct do j=1,3 @@ -180,33 +189,42 @@ C c include 'COMMON.CONTACTS' double precision gg(3) double precision evdw,evdwij - integer i,j,k,itypi,itypj,itypi1,num_conti,iint + integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont 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 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) 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 C Calculate SC interaction energy. C - do iint=1,nint_gr(i) +c do iint=1,nint_gr(i) cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) - do j=istart(i,iint),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) @@ -235,8 +253,8 @@ C gvdwc(k,j)=gvdwc(k,j)+gg(k) enddo endif - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i do i=1,nct do j=1,3 @@ -274,30 +292,39 @@ C include "COMMON.SPLITELE" double precision gg(3) double precision evdw,evdwij - integer i,j,k,itypi,itypj,itypi1,iint + integer i,j,k,itypi,itypj,itypi1,iint,ikont double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & 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 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) 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 - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +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 @@ -340,8 +367,8 @@ C gvdwc(k,j)=gvdwc(k,j)+gg(k) enddo endif - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i do i=1,nct do j=1,3 @@ -370,30 +397,39 @@ C include "COMMON.SPLITELE" double precision gg(3) double precision evdw,evdwij - integer i,j,k,itypi,itypj,itypi1,iint + integer i,j,k,itypi,itypj,itypi1,iint,ikont double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & 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 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) 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 - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +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 @@ -430,8 +466,8 @@ C gvdwc(k,j)=gvdwc(k,j)+gg(k) enddo endif - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i do i=1,nct do j=1,3 @@ -462,10 +498,11 @@ C integer icall common /srutu/ icall double precision evdw - integer itypi,itypj,itypi1,iint,ind + integer itypi,itypj,itypi1,iint,ind,ikont 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 @@ -477,13 +514,17 @@ c else lprn=.false. c endif ind=0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) 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) @@ -492,8 +533,8 @@ c dsci_inv=dsc_inv(itypi) C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -508,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) @@ -559,8 +604,8 @@ C Calculate the angular part of the gradient and sum add the contributions C to the appropriate components of the Cartesian gradient. call sc_grad_scale((1.0d0-sss)*sss1) endif - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i c stop return @@ -586,9 +631,10 @@ C integer icall common /srutu/ icall double precision evdw - integer itypi,itypj,itypi1,iint,ind + 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 @@ -600,13 +646,17 @@ c else lprn=.false. c endif ind=0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) 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) @@ -615,8 +665,8 @@ c dsci_inv=dsc_inv(itypi) C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -631,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) @@ -678,8 +732,8 @@ C Calculate the angular part of the gradient and sum add the contributions C to the appropriate components of the Cartesian gradient. call sc_grad_scale(sss) endif - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i c stop return @@ -706,13 +760,14 @@ C logical lprn integer xshift,yshift,zshift double precision evdw - integer itypi,itypj,itypi1,iint,ind + 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 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 @@ -720,19 +775,18 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon lprn=.false. c if (icall.eq.0) lprn=.false. ind=0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(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 + 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) @@ -743,8 +797,8 @@ c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -766,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 @@ -898,8 +900,8 @@ C Calculate the radial part of the gradient C Calculate angular part of the gradient. call sc_grad_scale((1.0d0-sss)*sss1) endif - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i c write (iout,*) "Number of loop steps in EGB:",ind cccc energy_dec=.false. @@ -925,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 + 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 @@ -941,19 +941,18 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon lprn=.false. c if (icall.eq.0) lprn=.false. ind=0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(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 + 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) @@ -964,8 +963,8 @@ c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -987,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) @@ -1115,8 +1062,8 @@ C Calculate the radial part of the gradient C Calculate angular part of the gradient. call sc_grad_scale(sss) endif - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i c write (iout,*) "Number of loop steps in EGB:",ind cccc energy_dec=.false. @@ -1143,7 +1090,7 @@ C integer icall common /srutu/ icall logical lprn - integer itypi,itypj,itypi1,iint,ind + integer itypi,itypj,itypi1,iint,ind,ikont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij, & xi,yi,zi,fac_augm,e_augm double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, @@ -1157,13 +1104,18 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon lprn=.false. c if (icall.eq.0) lprn=.true. ind=0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) 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) @@ -1172,8 +1124,8 @@ c dsci_inv=dsc_inv(itypi) C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -1190,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) @@ -1257,8 +1218,8 @@ C Calculate the radial part of the gradient C Calculate angular part of the gradient. call sc_grad_scale((1.0d0-sss)*sss1) endif - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i end C----------------------------------------------------------------------------- @@ -1282,7 +1243,7 @@ C integer icall common /srutu/ icall logical lprn - integer itypi,itypj,itypi1,iint,ind + integer itypi,itypj,itypi1,iint,ind,ikont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij, & xi,yi,zi,fac_augm,e_augm double precision evdw @@ -1290,13 +1251,17 @@ 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 lprn=.false. c if (icall.eq.0) lprn=.true. ind=0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -1306,13 +1271,15 @@ c if (icall.eq.0) lprn=.true. 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 C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -1329,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) @@ -1390,8 +1366,8 @@ C Calculate the radial part of the gradient C Calculate angular part of the gradient. call sc_grad_scale(sss) endif - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i end C---------------------------------------------------------------------------- @@ -1481,6 +1457,7 @@ C include 'COMMON.TIME1' include 'COMMON.SHIELD' include "COMMON.SPLITELE" + integer ikont 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), @@ -1581,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) @@ -1610,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 @@ -1629,7 +1596,10 @@ C & .or. itype(i-1).eq.ntyp1 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 +c do i=iatel_s,iatel_e + do ikont=g_listpp_start,g_listpp_end + i=newcontlistppi(ikont) + j=newcontlistppj(ikont) 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 @@ -1643,23 +1613,18 @@ 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) #endif - do j=ielstart(i),ielend(i) +c do j=ielstart(i),ielend(i) 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 +c enddo ! j #ifdef FOURBODY num_cont_hb(i)=num_conti #endif @@ -1730,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 @@ -1762,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) @@ -2700,12 +2632,17 @@ 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" c write (iout,*) "iatel_s_vdw",iatel_s_vdw, c & " iatel_e_vdw",iatel_e_vdw c call flush(iout) - do i=iatel_s_vdw,iatel_e_vdw +c do i=iatel_s_vdw,iatel_e_vdw + do ikont=g_listpp_vdw_start,g_listpp_vdw_end + i=newcontlistpp_vdwi(ikont) + j=newcontlistpp_vdwj(ikont) if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -2716,17 +2653,12 @@ c call flush(iout) 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) c call flush(iout) - do j=ielstart_vdw(i),ielend_vdw(i) +c do j=ielstart_vdw(i),ielend_vdw(i) if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle ind=ind+1 iteli=itel(i) @@ -2743,43 +2675,10 @@ c call flush(iout) 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) @@ -2817,7 +2716,7 @@ C ggg(3)=facvdw*zj gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo endif - enddo ! j +c enddo ! j enddo ! i return end @@ -2843,14 +2742,13 @@ 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 evdw2_14=0.0d0 @@ -2859,22 +2757,20 @@ 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 +c do i=iatscp_s,iatscp_e + do ikont=g_listscp_start,g_listscp_end + i=newcontlistscpi(ikont) + j=newcontlistscpj(ikont) 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 + call to_box(xi,yi,zi) - do iint=1,nscp_gr(i) +c do iint=1,nscp_gr(i) - do j=iscpstart(i,iint),iscpend(i,iint) +c do j=iscpstart(i,iint),iscpend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle C Uncomment following three lines for SC-p interactions @@ -2886,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) @@ -2969,9 +2831,9 @@ c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo endif - enddo +c enddo - enddo ! iint +c enddo ! iint enddo ! i do i=1,nct do j=1,3 @@ -3021,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' @@ -3034,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), @@ -3060,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)))