Add neighbors in shell

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 3.1.12
+current_version = 3.2.1
 commit = True
 tag = True
 

.github/workflows/deploy.yml

@@ -14,8 +14,8 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        os: [ubuntu-latest, macos-latest]
-        python-version: [3.9, '3.10', '3.11']
+        os: [ubuntu-latest, macos-13]
+        python-version: ['3.10', '3.11']
 
     steps:
       - uses: actions/checkout@v3

.github/workflows/testing.yml

@@ -10,8 +10,8 @@ jobs:
   build:
     strategy:
       matrix:
-        os: [ubuntu-latest, macos-latest, windows-latest]
-        python-version: [3.9, '3.10', '3.11']
+        os: [ubuntu-latest, macos-13, windows-latest]
+        python-version: ['3.10', '3.11']
 
     runs-on: ${{ matrix.os }}
     steps:
@@ -26,7 +26,7 @@ jobs:
         key:
           ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-${{
           hashFiles('.ci_support/environment.yml') }}
-    - uses: conda-incubator/setup-miniconda@v2
+    - uses: conda-incubator/setup-miniconda@v3
       with:
         activate-environment: pyscal-test
         channel-priority: strict

setup.py

@@ -7,7 +7,7 @@
 
 setup(
     name='pyscal3',
-    version='3.1.12',
+    version='3.2.1',
     author='Sarath Menon',
     author_email='sarath.menon@pyscal.org',
     description='Python library written in C++ for calculation of local atomic structural environment',

src/pyscal3/neighbor.cpp

@@ -384,6 +384,100 @@ void get_all_neighbors_normal(py::dict& atoms,
     atoms[py::str("cutoff")] = cutoff;
 }
 
+void get_all_neighbors_shell_normal(py::dict& atoms,
+    const double dmin,
+    const double dmax,
+    const int triclinic,
+    const vector<vector<double>> rot, 
+    const vector<vector<double>> rotinv,
+    const vector<double> box)
+    {
+
+    double d;
+    double diffx,diffy,diffz;
+    double tempr,temptheta,tempphi;
+    vector<double> diffi, diffj;
+
+    //access positions and put it in an array
+    vector<vector<double>> positions = atoms[py::str("positions")].cast<vector<vector<double>>>();
+    vector<bool> mask_1 = atoms[py::str("mask_1")].cast<vector<bool>>();
+    vector<bool> mask_2 = atoms[py::str("mask_2")].cast<vector<bool>>();
+    //auto positions = atoms[py::str("positions")].cast<py::array_t<double>>();
+
+    int nop = positions.size();
+    vector<vector<int>> neighbors(nop);
+    vector<vector<double>> neighbordist(nop);
+    vector<vector<double>> neighborweight(nop);
+    vector<vector<vector<double>>> diff(nop);
+    vector<vector<double>> r(nop);
+    vector<vector<double>> phi(nop);
+    vector<vector<double>> theta(nop);
+    vector<double> cutoff(nop); 
+
+    //now loop and calculate things
+    for (int ti=0; ti<nop; ti++){
+        //only consider ti if mask_1 is false
+        if (mask_1[ti]) continue;
+        for (int tj=ti+1; tj<nop; tj++){
+            //mask 2 is applied here
+            if (mask_2[tj]) continue;
+            d = get_abs_distance(positions[ti], positions[tj],
+                triclinic, rot, rotinv, box, 
+                diffx, diffy, diffz);
+            if ((d >= dmin) && (d <= dmax)){          
+                neighbors[ti].emplace_back(tj);
+                neighbors[tj].emplace_back(ti);
+
+                neighbordist[ti].emplace_back(d);
+                neighbordist[tj].emplace_back(d);
+
+                neighborweight[ti].emplace_back(1.00);
+                neighborweight[tj].emplace_back(1.00);
+
+                diffi.clear();
+                diffi.emplace_back(diffx);
+                diffi.emplace_back(diffy);
+                diffi.emplace_back(diffz);
+
+                diffj.clear();
+                diffj.emplace_back(-diffx);
+                diffj.emplace_back(-diffy);
+                diffj.emplace_back(-diffz);
+
+                diff[ti].emplace_back(diffi);
+                diff[tj].emplace_back(diffj);
+
+                convert_to_spherical_coordinates(diffx, diffy, diffz, tempr, tempphi, temptheta);
+
+                r[ti].emplace_back(tempr);
+                phi[ti].emplace_back(tempphi);
+                theta[ti].emplace_back(temptheta);
+                //n_neighbors += 1;
+                cutoff[ti] = dmax;
+
+                convert_to_spherical_coordinates(-diffx, -diffy, -diffz, tempr, tempphi, temptheta);
+
+                r[tj].emplace_back(tempr);
+                phi[tj].emplace_back(tempphi);
+                theta[tj].emplace_back(temptheta);
+                //atoms[tj].n_neighbors += 1;
+                cutoff[tj] = dmax;
+
+            }
+        }
+    }
+
+    //calculation over lets assign
+    atoms[py::str("neighbors")] = neighbors;
+    atoms[py::str("neighbordist")] = neighbordist;
+    atoms[py::str("neighborweight")] = neighborweight;
+    atoms[py::str("diff")] = diff;
+    atoms[py::str("r")] = r;
+    atoms[py::str("theta")] = theta;
+    atoms[py::str("phi")] = phi;
+    atoms[py::str("cutoff")] = cutoff;
+}
+
 
 int cell_index(int cx, int cy, int cz, int nx, int ny, int nz){
     return cx*ny*nz + cy*nz + cz;
@@ -603,6 +697,112 @@ void get_all_neighbors_cells(py::dict& atoms,
 }
 
 
+void get_all_neighbors_shell_cells(py::dict& atoms,
+    const double dmin,
+    const double dmax,
+    const int triclinic,
+    const vector<vector<double>> rot, 
+    const vector<vector<double>> rotinv,
+    const vector<double> box){
+    double d;
+    double diffx,diffy,diffz;
+    double tempr,temptheta,tempphi;
+    vector<double> diffi, diffj;
+    int ti, tj;
+
+    //access positions and put it in an array
+    vector<vector<double>> positions = atoms[py::str("positions")].cast<vector<vector<double>>>();
+    vector<bool> mask_1 = atoms[py::str("mask_1")].cast<vector<bool>>();
+    vector<bool> mask_2 = atoms[py::str("mask_2")].cast<vector<bool>>();
+    int nop = positions.size();
+    vector<vector<int>> neighbors(nop);
+    vector<vector<double>> neighbordist(nop);
+    vector<vector<double>> neighborweight(nop);
+    vector<vector<vector<double>>> diff(nop);
+    vector<vector<double>> r(nop);
+    vector<vector<double>> phi(nop);
+    vector<vector<double>> theta(nop);
+    vector<double> cutoff(nop); 
+    vector<cell> cells = set_up_cells(positions, box, dmax);
+    int total_cells = cells.size();
+    int subcell;
+    //now loop to find distance
+    for(int i=0; i<total_cells; i++){
+        //for each member in cell i
+        for(size_t mi=0; mi<cells[i].members.size(); mi++){
+            //now go through the neighbors
+            ti = cells[i].members[mi];
+            if (mask_1[ti]) continue;
+            for(size_t j=0 ; j<cells[i].neighbor_cells.size(); j++){
+               //loop through members of j
+               subcell = cells[i].neighbor_cells[j];
+               for(size_t mj=0; mj<cells[subcell].members.size(); mj++){
+                    //now we have mj -> members/compare with
+                    tj = cells[subcell].members[mj];
+                    if (mask_2[tj]) continue;
+                    //compare ti and tj and add
+                    if (ti < tj){
+                        d = get_abs_distance(positions[ti], positions[tj],
+                            triclinic, rot, rotinv, box, 
+                            diffx, diffy, diffz);
+                        if ((d >= dmin) && (d <= dmax)){
+                            //cout<<"adding "<<ti<<" to "<<tj<<endl;
+                            //cout<<"adding "<<tj<<" to "<<ti<<endl;
+                            neighbors[ti].emplace_back(tj);
+                            neighbors[tj].emplace_back(ti);
+
+                            neighbordist[ti].emplace_back(d);
+                            neighbordist[tj].emplace_back(d);
+
+                            neighborweight[ti].emplace_back(1.00);
+                            neighborweight[tj].emplace_back(1.00);
+
+                            diffi.clear();
+                            diffi.emplace_back(diffx);
+                            diffi.emplace_back(diffy);
+                            diffi.emplace_back(diffz);
+
+                            diffj.clear();
+                            diffj.emplace_back(-diffx);
+                            diffj.emplace_back(-diffy);
+                            diffj.emplace_back(-diffz);
+
+                            diff[ti].emplace_back(diffi);
+                            diff[tj].emplace_back(diffj);
+
+                            convert_to_spherical_coordinates(diffx, diffy, diffz, tempr, tempphi, temptheta);
+
+                            r[ti].emplace_back(tempr);
+                            phi[ti].emplace_back(tempphi);
+                            theta[ti].emplace_back(temptheta);
+                            //n_neighbors += 1;
+                            cutoff[ti] = dmax;
+
+                            convert_to_spherical_coordinates(-diffx, -diffy, -diffz, tempr, tempphi, temptheta);
+
+                            r[tj].emplace_back(tempr);
+                            phi[tj].emplace_back(tempphi);
+                            theta[tj].emplace_back(temptheta);
+                            //atoms[tj].n_neighbors += 1;
+                            cutoff[tj] = dmax;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    //calculation over lets assign
+    atoms[py::str("neighbors")] = neighbors;
+    atoms[py::str("neighbordist")] = neighbordist;
+    atoms[py::str("neighborweight")] = neighborweight;
+    atoms[py::str("diff")] = diff;
+    atoms[py::str("r")] = r;
+    atoms[py::str("theta")] = theta;
+    atoms[py::str("phi")] = phi;
+    atoms[py::str("cutoff")] = cutoff;
+}
+
+
 void get_temp_neighbors_brute(const vector<vector<double>>& positions,
     const vector<bool>& mask_1,
     const vector<bool>& mask_2,

src/pyscal3/operations/identify.py

@@ -1,7 +1,7 @@
 import numpy as np
 import pyscal3.csystem as pc
 
-def find_neighbors(system, method='cutoff', cutoff=0, threshold=2, 
+def find_neighbors(system, method='cutoff', cutoff=0, shell_thickness=0, threshold=2, 
         voroexp=1, padding=1.2, nlimit=6, 
         cells=None, nmax=12, assign_neighbor=True):
     """
@@ -20,6 +20,10 @@ def find_neighbors(system, method='cutoff', cutoff=0, threshold=2,
         the cutoff distance to be used for the `cutoff` based neighbor calculation method described above.
         If the value is specified as 0 or `adaptive`, adaptive method is used.
         If the value is specified as `sann`, sann algorithm is used.
+
+    shell_thickness : float, optional
+        only used with cutoff method. If provided, neighbors are found between cutoff, and cutoff + shell_thickness.
+        Default 0.
 
     threshold : float, optional
         only used if ``cutoff=adaptive``. A threshold which is used as safe limit for calculation of cutoff.
@@ -128,20 +132,31 @@ def find_neighbors(system, method='cutoff', cutoff=0, threshold=2,
                 raise RuntimeError("sann cutoff could not be converged. This is most likely, \
                     due to a low threshold value. Try increasing it.")
 
-        elif cutoff=='adaptive' or cutoff==0:
+        elif cutoff=='adaptive' or (cutoff==0 and shell_thickness==0):
             finished = pc.get_all_neighbors_adaptive(system.atoms, 0.0,
                 system.triclinic, system.rot, system.rotinv,
                 system.boxdims, threshold, nlimit, padding, cells)
             if not bool(finished):
                 raise RuntimeError("Could not find adaptive cutoff")
         else:
-            if cells:
-                pc.get_all_neighbors_cells(system.atoms, cutoff,
-                    system.triclinic, system.rot, system.rotinv, system.boxdims)
+            if (cutoff==0 and shell_thickness>0):
+                cutoff = shell_thickness
+                shell_thickness = 0
+            if shell_thickness == 0:
+                if cells:
+                    pc.get_all_neighbors_cells(system.atoms, cutoff,
+                        system.triclinic, system.rot, system.rotinv, system.boxdims)
+                else:
+                    pc.get_all_neighbors_normal(system.atoms, cutoff,
+                        system.triclinic, system.rot, system.rotinv, system.boxdims)
             else:
-                pc.get_all_neighbors_normal(system.atoms, cutoff,
-                    system.triclinic, system.rot, system.rotinv, system.boxdims)
-
+                if cells:
+                    pc.get_all_neighbors_shell_cells(system.atoms, cutoff, cutoff+shell_thickness,
+                        system.triclinic, system.rot, system.rotinv, system.boxdims)
+                else:
+                    pc.get_all_neighbors_shell_normal(system.atoms, cutoff, cutoff+shell_thickness,
+                        system.triclinic, system.rot, system.rotinv, system.boxdims)
+
     elif method == 'number':
         finished = pc.get_all_neighbors_bynumber(system.atoms, 0.0, 
             system.triclinic, system.rot, system.rotinv,

src/pyscal3/system.h

@@ -84,6 +84,14 @@ void get_all_neighbors_normal(py::dict& atoms,
     const vector<vector<double>> rotinv,
     const vector<double> box);
 
+void get_all_neighbors_shell_normal(py::dict& atoms,
+    const double dmin,
+    const double dmax,
+    const int triclinic,
+    const vector<vector<double>> rot,
+    const vector<vector<double>> rotinv,
+    const vector<double> box);
+
 int cell_index(int, int, int, int, int, int);
 vector<int> cell_periodic(int, int, int, int, int, int);
 
@@ -98,6 +106,14 @@ void get_all_neighbors_cells(py::dict& atoms,
     const vector<vector<double>> rotinv,
     const vector<double> box);
 
+void get_all_neighbors_shell_cells(py::dict& atoms,
+    const double dmin,
+    const double dmax,
+    const int triclinic,
+    const vector<vector<double>> rot,
+    const vector<vector<double>> rotinv,
+    const vector<double> box);
+
 void get_temp_neighbors_brute(const vector<vector<double>>& positions,
     const vector<bool>& mask_1,
     const vector<bool>& mask_2,

src/pyscal3/system_binding.cpp

@@ -22,7 +22,9 @@ PYBIND11_MODULE(csystem, m) {
     m.def("remap_and_displace_atom", &remap_and_displace_atom);
     m.def("reset_all_neighbors", &reset_all_neighbors);
     m.def("get_all_neighbors_normal", &get_all_neighbors_normal);
+    m.def("get_all_neighbors_shell_normal", &get_all_neighbors_shell_normal);
     m.def("get_all_neighbors_cells", &get_all_neighbors_cells);
+    m.def("get_all_neighbors_shell_cells", &get_all_neighbors_shell_cells);
     m.def("get_all_neighbors_bynumber", &get_all_neighbors_bynumber);
     m.def("get_all_neighbors_sann", &get_all_neighbors_sann);
     m.def("get_all_neighbors_adaptive", &get_all_neighbors_adaptive);

tests/test_neighbors.py

@@ -70,4 +70,16 @@ def test_system_init():
  2.7080614376339383,
  2.7080614376339383,
  2.708061437633939])
-    assert np.sum(a1-a2) < 1E-5
+    assert np.sum(a1-a2) < 1E-5
+
+
+def test_neighbor_shell():
+   sys = pc.System.create.element.Cu(repetitions=(5,5,5))
+   sys.find.neighbors(method='cutoff', cutoff=3, shell_thickness=1, cells=False)
+   assert len(sys.atoms.neighbors.distance[0]) == 6
+   assert np.abs(sys.atoms.neighbors.distance[0][0]-3.61) <= 1E-3
+
+   sys = pc.System.create.element.Cu(repetitions=(5,5,5))
+   sys.find.neighbors(method='cutoff', cutoff=3, shell_thickness=1, cells=True)
+   assert len(sys.atoms.neighbors.distance[0]) == 6
+   assert np.abs(sys.atoms.neighbors.distance[0][0]-3.61) <= 1E-3