fix a fatal bug in mk_lspr_sngl, which could affect some of python pr…

…ograms.
ryokbys · Nov 9, 2023 · 782b229 · 782b229
1 parent bf40c9f
commit 782b229
Show file tree

Hide file tree

Showing 2 changed files with 30 additions and 12 deletions.
diff --git a/nappy/adf.py b/nappy/adf.py
@@ -74,10 +74,10 @@ def adf_atom(ia,dang,rcut,nsys,poss,lspr,symbols,sj,sk):
         if rij2 >= rcut2:
             continue
         rij= np.sqrt(rij2)
-        for ki in range(len(lspri)):
+        for ki in range(ji+1,len(lspri)):
             ka= lspri[ki]
-            if ka == ia or ka <= ja:
-                continue
+            # if ka == ia or ka <= ja:
+            #     continue
             ski = symbols[ka]
             if set((sji,ski)) != set((sj,sk)):
                 continue
@@ -96,7 +96,6 @@ def adf_atom(ia,dang,rcut,nsys,poss,lspr,symbols,sj,sk):
             else:
                 rad= np.arccos(cs)
             deg= rad/np.pi *180.0
-
             nda[int(deg/dang)] += 1
     return nda
 

diff --git a/pmd/mod_pairlist.F90 b/pmd/mod_pairlist.F90
@@ -544,7 +544,7 @@ subroutine mk_lspr_sngl(namax,natm,nnmax,tag,ra,rc,h,hi &
 
     integer:: i,j,k,l,m,n
     integer:: mx,my,mz,kux,kuy,kuz,m1x,m1y,m1z,m1,ic,jc,ierr
-    real(8):: xi(3),xij(3),rij(3),rij2
+    real(8):: xi(3),xij(3),rij(3),rij2,shft(3)
 
 !!$    integer,allocatable,save:: lscl(:),lshd(:)
 !!$    real(8),save:: rc2,rcx,rcy,rcz,rcxi,rcyi,rczi,rc1nn2
@@ -611,17 +611,35 @@ subroutine mk_lspr_sngl(namax,natm,nnmax,tag,ra,rc,h,hi &
       mz= min(max(mz,1),lcz)
       m= (mx-1)*lcyz +(my-1)*lcz +mz
       do kux=-1,1
+        shft(1) = 0d0
         m1x = mx +kux
-        if( m1x.lt.1   ) m1x= m1x +lcx
-        if( m1x.gt.lcx ) m1x= m1x -lcx
+        if( m1x.lt.1   ) then
+          m1x= m1x +lcx
+          shft(1) = -1d0
+        else if( m1x.gt.lcx ) then
+          m1x= m1x -lcx
+          shft(1) = +1d0
+        endif
         do kuy=-1,1
+          shft(2) = 0d0
           m1y= my +kuy
-          if( m1y.lt.1   ) m1y= m1y +lcy
-          if( m1y.gt.lcy ) m1y= m1y -lcy
+          if( m1y.lt.1   ) then
+            m1y= m1y +lcy
+            shft(2) = -1d0
+          else if( m1y.gt.lcy ) then
+            m1y= m1y -lcy
+            shft(2) = +1d0
+          endif
           do kuz=-1,1
+            shft(3) = 0d0
             m1z= mz +kuz
-            if( m1z.lt.1   ) m1z= m1z +lcz
-            if( m1z.gt.lcz ) m1z= m1z -lcz
+            if( m1z.lt.1   ) then
+              m1z= m1z +lcz
+              shft(3) = -1d0
+            else if( m1z.gt.lcz ) then
+              m1z= m1z -lcz
+              shft(3) = +1d0
+            endif
             m1= (m1x-1)*lcyz +(m1y-1)*lcz +m1z
             if( lshd(m1).eq.0 ) cycle
             j = lshd(m1)
@@ -630,7 +648,8 @@ subroutine mk_lspr_sngl(namax,natm,nnmax,tag,ra,rc,h,hi &
                 j = lscl(j)
                 cycle
               endif
-              xij(1:3)= ra(1:3,j)-xi(1:3) -anint(ra(1:3,j)-xi(1:3))
+              xij(1:3)= ra(1:3,j) -xi(1:3) +shft(1:3)
+!!$              xij(1:3)= xij(1:3) -anint(xij(1:3))  !! This is wrong in this case.
               rij(1:3)=h(1:3,1)*xij(1) +h(1:3,2)*xij(2) +h(1:3,3)*xij(3)
               rij2= rij(1)**2 +rij(2)**2 +rij(3)**2
               if( rij2.lt.rc2 ) then