X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Felecont.f;h=73325f2fa65e9a380e1e9406c1a826983c405528;hb=33ca65593386f7dce46c59efeb0f77b996757cbf;hp=4faf3e584181c3a6a5d1de0b78aaab4e2ce96269;hpb=0bb81c1c3180a2079d70af7dd534295e1e0b1e4c;p=unres.git diff --git a/source/unres/src_MD-M/elecont.f b/source/unres/src_MD-M/elecont.f index 4faf3e5..73325f2 100644 --- a/source/unres/src_MD-M/elecont.f +++ b/source/unres/src_MD-M/elecont.f @@ -52,6 +52,12 @@ c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ xmedi=xi+0.5*dxi ymedi=yi+0.5*dyi zmedi=zi+0.5*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 do 4 j=i+2,nct-1 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4 ind=ind+1 @@ -66,9 +72,49 @@ c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ dxj=c(1,j+1)-c(1,j) dyj=c(2,j+1)-c(2,j) dzj=c(3,j+1)-c(3,j) - xj=c(1,j)+0.5*dxj-xmedi - yj=c(2,j)+0.5*dyj-ymedi - zj=c(3,j)+0.5*dzj-zmedi + xj=c(1,j)+0.5*dxj + yj=c(2,j)+0.5*dyj + zj=c(3,j)+0.5*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-xi)**2+(yj-yi)**2+(zj-zi)**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-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 + 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 + sss=sscale(sqrt(rij)) + sssgrad=sscagrad(sqrt(rij)) rrmij=1.0/(xj*xj+yj*yj+zj*zj) rmij=sqrt(rrmij) r3ij=rrmij*rmij @@ -94,7 +140,7 @@ c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ econt(ncont)=eesij endif ees=ees+eesij - evdw=evdw+evdwij + evdw=evdw+evdwij*sss 4 continue 1 continue if (lprint) then