A few fixes to Pylada

CMakeLists.txt

@@ -1,4 +1,4 @@
-CMAKE_MINIMUM_REQUIRED(VERSION 3.12)
+CMAKE_MINIMUM_REQUIRED(VERSION 3.20.1)
 project(Pylada
     VERSION 1.0.1
     DESCRIPTION "A pythonic computation material science platform"
@@ -49,21 +49,14 @@ endif()
 
 set(OLD_CMAKE_FIND_FRAMEWORK "${CMAKE_FIND_FRAMEWORK}")
 set(CMAKE_FIND_FRAMEWORK LAST)
-find_package(Python3 COMPONENTS Development Interpreter)
+find_package(Python COMPONENTS Development.Module Interpreter REQUIRED)
 set(CMAKE_FIND_FRAMEWORK "${OLD_CMAKE_FIND_FRAMEWORK}")
-if(NOT Python3_FOUND)
+if(Python2_VERSION_MAJOR EQUAL 2)
     set(CMAKE_FIND_FRAMEWORK LAST)
-    find_package(Python2 COMPONENTS Development Interpreter)
     set(CMAKE_FIND_FRAMEWORK "${OLD_CMAKE_FIND_FRAMEWORK}")
-endif()
-if(NOT Python3_FOUND AND NOT Python2_FOUND)
-    message(FATAL_ERROR "Could find neither python3 nor python2")
-endif()
-
-if(Python3_FOUND)
-    set(PYVER PY3)
-else()
     set(PYVER PY2)
+else()
+    set(PYVER PY3)
 endif()
 
 if(NOT SKBUILD)
@@ -85,7 +78,6 @@ if(NOT SKBUILD)
     list(APPEND CMAKE_MODULE_PATH "${SKBUILD_LOCATION}/resources/cmake")
 endif()
 
-find_package(PythonExtensions REQUIRED)
 find_package(Cython REQUIRED)
 find_package(NumPy REQUIRED)
 find_package(Eigen3 NO_MODULE)

cmake_modules/eigen/CMakeLists.txt

@@ -1,5 +1,5 @@
 project(external CXX)
-cmake_minimum_required(VERSION 3.12)
+cmake_minimum_required(VERSION 3.20.1)
 include("ExternalProject")
 
 if(NOT EIGEN_TAR_ARCHIVE)

src/pylada/CMakeLists.txt

@@ -8,7 +8,6 @@ file(WRITE "${CMAKE_CURRENT_BINARY_DIR}/_version.py"
 )
 install(FILES "${CMAKE_CURRENT_BINARY_DIR}/_version.py" DESTINATION ${PY_ROOT_DIR})
 
-include_directories("${PYTHON_INCLUDE_DIR}")
 add_subdirectory(crystal)
 add_subdirectory(decorations)
 add_subdirectory(ewald)

src/pylada/crystal/CMakeLists.txt

@@ -30,18 +30,21 @@ set(cython_sources
 foreach(filename ${cython_sources})
     get_filename_component(name "${filename}" NAME_WE)
     add_cython_target(cython_${name} ${filename} CXX ${PYVER})
-    add_library(${name} MODULE ${cython_${name}})
-    python_extension_module(${name})
-    target_link_libraries(${name} Eigen3::Eigen)
-    target_include_directories(${name} PRIVATE "${PY_HEADER_DIR}" "${NumPy_INCLUDE_DIRS}")
+
+    python_add_library(${name} MODULE WITH_SOABI ${cython_${name}})
+
+    target_link_libraries(${name} PRIVATE Eigen3::Eigen)
+    target_include_directories(${name} PRIVATE "${PY_HEADER_DIR}" "${NumPy_INCLUDE_DIRS}" "${Python3_INCLUDE_DIRS}")
     install(TARGETS ${name} LIBRARY DESTINATION ${PY_ROOT_DIR}/crystal)
 endforeach()
 
 add_cython_target(cython_cutilities cutilities.pyx CXX ${PYVER})
-add_library(cutilities MODULE ${cython_cutilities} smith_normal_form.cc gruber.cc noopt.cc)
-python_extension_module(cutilities)
-target_include_directories(cutilities PRIVATE "${PY_HEADER_DIR}" "${NumPy_INCLUDE_DIRS}")
-target_link_libraries(cutilities Eigen3::Eigen)
+
+python_add_library(cutilities MODULE WITH_SOABI ${cython_cutilities} smith_normal_form.cc gruber.cc noopt.cc)
+
+
+target_include_directories(cutilities PRIVATE "${PY_HEADER_DIR}" "${NumPy_INCLUDE_DIRS}" "${Python3_INCLUDE_DIRS}")
+target_link_libraries(cutilities PRIVATE Eigen3::Eigen)
 install(TARGETS cutilities LIBRARY DESTINATION ${PY_ROOT_DIR}/crystal)
 
 add_subdirectory(defects)

src/pylada/crystal/_coordination_shells.pyx

@@ -65,7 +65,7 @@ def coordination_shells(structure, int nshells, center, tolerance=1e-12, natoms=
     from . import gruber, into_voronoi
 
     cdef int N = len(structure)
-    cell = gruber(structure.cell, 3e0*tolerance)
+    cell = gruber(structure.cell, tolerance = 3e0*tolerance)
     invcell = inv(cell)
     cdef double volume = structure.volume
     cdef int maxatoms = __natoms(natoms, nshells)
@@ -154,7 +154,7 @@ def __neighbors(structure, int nmax, center, tolerance=1e-12, natoms=0):
     from . import gruber, into_voronoi
 
     cdef int N = len(structure)
-    cell = gruber(structure.cell, 3e0*tolerance)
+    cell = gruber(structure.cell, tolerance = 3e0*tolerance)
     invcell = inv(cell)
     cdef double volume = structure.volume
     cdef int maxatoms = max([natoms, nmax + 2])

src/pylada/crystal/_primitive.pyx

@@ -20,10 +20,10 @@
 #  <http://www.gnu.org/licenses/>.
 ###############################
 
-cdef __translations(structure, double tolerance):
+def __translations(structure, double tolerance):
     """ Looks for internal translations """
     from numpy.linalg import inv
-    from numpy import all, abs, allclose
+    from numpy import all, abs, allclose, mod
     from . import gruber, into_cell, into_voronoi
 
     cell = gruber(structure.cell)
@@ -32,32 +32,36 @@ cdef __translations(structure, double tolerance):
     front_type = structure[0].type
     center = structure[0].pos
     translations = []
+
     for site in structure:
         if front_type != site.type:
             continue
 
-        translation = into_voronoi(site.pos - center, cell, invcell)
+        translation = into_voronoi(site.pos - center + tolerance, cell, invcell) - tolerance
+
         if all(abs(translation) < tolerance):
             continue
 
         for mapping in structure:
-            pos = into_cell(mapping.pos + translation, cell, invcell)
+            pos = into_cell(mapping.pos + translation + tolerance, cell, invcell) - tolerance
             for mappee in structure:
-                if mapping.type == mappee.type and allclose(mappee.pos, pos, tolerance):
+                if mapping.type == mappee.type and allclose(into_cell(mappee.pos + tolerance, cell, invcell)-tolerance, pos, atol=tolerance):
                     break
             else:
                 break
         else:
-            translations.append(into_voronoi(translation, cell, invcell))
+            yield into_voronoi(translation, cell, invcell)
 
-    return translations
+    yield cell[:, 0]
+    yield cell[:, 1]
+    yield cell[:, 2]
 
 
-def primitive(structure, double tolerance=1e-8):
+def primitive(structure, double tolerance=1e-8, bint recursive=False):
     """ Tries to compute the primitive cell of the input structure """
     # The tolerance is the absolute tolerance on the translations
     from numpy.linalg import inv, det
-    from numpy import all, abs, array, dot, allclose, round
+    from numpy import all, abs, array, dot, allclose, round, mod
     from . import gruber, into_cell, into_voronoi, into_cell
     from .. import error, logger
 
@@ -68,32 +72,29 @@ def primitive(structure, double tolerance=1e-8):
     cell = gruber(result.cell)
     invcell = inv(cell)
     for atom in result:
-        atom.pos = into_cell(atom.pos, cell, invcell)
+        atom.pos = into_cell(atom.pos + tolerance, cell, invcell) - tolerance
 
-    translations = __translations(result, tolerance)
-    if len(translations) == 0:
+    for i,t in enumerate(__translations(result, tolerance)):
+        if i > 2:
+            break
+    else:
         logger.debug("Found no inner translations: structure is primitive")
         for a, b in zip(result, structure):
             a.pos[:] = b.pos[:]
         return result
-
-    # adds original translations.
-    translations.append(cell[:, 0])
-    translations.append(cell[:, 1])
-    translations.append(cell[:, 2])
-
+
     # Looks for cell with smallest volume 
     new_cell = result.cell.copy()
     volume = abs(det(new_cell))
-    for i, first in enumerate(translations):
-        for j, second in enumerate(translations):
+    for i, first in enumerate(__translations(result, tolerance)):
+        for j, second in enumerate(__translations(result, tolerance)):
             if i == j:
                 continue
-            for k, third in enumerate(translations):
+            for k, third in enumerate(__translations(result, tolerance)):
                 if i == k or j == k:
                     continue
                 trial = array([first, second, third]).T
-                if abs(det(trial)) < abs(volume)/len(structure)*(1 - 3*tolerance): # The volume of the cell cannot be smaller than the specific volume 
+                if abs(det(trial)) < abs(result.volume)/len(structure)*(1 - 3*tolerance): # The volume of the cell cannot be smaller than the specific volume
                     continue
                 if abs(det(trial)) > volume*(1 - 3.0 * tolerance):
                     continue
@@ -106,30 +107,39 @@ def primitive(structure, double tolerance=1e-8):
                         logger.error(msg)
                         raise error.RuntimeError(msg)
                 integer_cell = dot(inv(trial), cell)
-                if allclose(integer_cell, round(integer_cell + 1e-7), tolerance): # This needs to be a supercell within the specified tolerance
+                if allclose(integer_cell, round(integer_cell + 1e-7), atol = tolerance): # This needs to be a supercell within the specified tolerance
                     new_cell = trial
                     volume = abs(det(trial))
-
+                    if recursive:
+                        break
+            else:
+                continue
+            break
+        else:
+            continue
+        break
+
     # Found the new cell with smallest volume (e.g. primivite)
     if abs(structure.volume - volume) < 1e-12: # It actually did not change
         msg = "Found translation but no primitive cell."
         logger.error(msg)
         raise error.RuntimeError(msg)
-
+    
     # now creates new lattice.
     result.clear()
     logger.debug("Found potential cell {!r}".format(new_cell))
     result.cell = gruber(new_cell)
     invcell = inv(result.cell)
     for site in structure:
-        pos = into_cell(site.pos, result.cell, invcell)
+        pos = result.cell.dot(mod(invcell.dot(site.pos) + tolerance, 1) - tolerance)
+        # pos = into_cell(site.pos + tolerance, result.cell, invcell)-tolerance
         for unique in result:
-            if site.type == unique.type and allclose(unique.pos, pos, abs(structure.volume / result.volume) * tolerance): # The difference between pos and unique pos is that of site.pos and the image of unique.pos
+            if site.type == unique.type and allclose(unique.pos, pos, atol = abs(structure.volume / result.volume) * tolerance): # The difference between pos and unique pos is that of site.pos and the image of unique.pos
                 break
         else:
             result.append(site.copy())
             result[-1].pos = pos
-
+            
     if len(structure) % len(result) != 0:
         msg = "Nb of atoms in output not multiple of input."
         logger.error(msg)
@@ -141,7 +151,11 @@ def primitive(structure, double tolerance=1e-8):
         raise error.RuntimeError(msg)
 
     logger.debug("Primitive structure found with %i/%i atoms" % (len(result), len(structure)))
-    return result;
+
+    if recursive:
+        return primitive(result, tolerance=tolerance, recursive=True)
+    else:
+        return result
 
 def is_primitive(structure, double tolerance = 1e-12):
     """ True if the lattice is primitive
@@ -165,4 +179,8 @@ def is_primitive(structure, double tolerance = 1e-12):
     for atom in result:
         atom.pos = into_cell(atom.pos, cell, invcell)
 
-    return len(__translations(result, tolerance)) == 0
+    for i,t in enumerate(__translations(result, tolerance)):
+        if i > 2:
+            return False
+    else:
+        return True

src/pylada/crystal/_space_group.pyx

@@ -21,6 +21,7 @@
 ###############################
 import numpy as np
 cimport numpy as np
+np.import_array()
 
 cpdef __gvectors(double[:, ::1] cell, double tolerance):
     """ Computes all gvectors in prolate defined by cell """

src/pylada/crystal/cutilities.pyx

@@ -1,5 +1,6 @@
 import numpy as np
 cimport numpy as np
+np.import_array()
 
 cdef extern from "pylada/crystal/types.h" namespace "pylada::types":
     ctypedef int t_int

src/pylada/crystal/defects/CMakeLists.txt

@@ -24,11 +24,12 @@
 ###############################
 
 add_cython_target(defects_cython _defects.pyx CXX ${PYVER})
-add_library(_defects MODULE ${defects_cython} third_order.cc)
-python_extension_module(_defects)
-target_link_libraries(_defects Eigen3::Eigen)
+
+python_add_library(_defects MODULE WITH_SOABI ${defects_cython} third_order.cc)
+
+target_link_libraries(_defects PRIVATE Eigen3::Eigen)
 target_include_directories(
     _defects PRIVATE
-    "${PY_HEADER_DIR}" "${NumPy_INCLUDE_DIRS}" "${PYTHON_INCLUDE_DIR}"
+    "${PY_HEADER_DIR}" "${NumPy_INCLUDE_DIRS}" "${Python3_INCLUDE_DIR}"
 )
 install(TARGETS _defects LIBRARY DESTINATION ${PY_ROOT_DIR}/crystal/defects)

src/pylada/crystal/read.py

@@ -89,12 +89,23 @@ def poscar(path="POSCAR", types=None):
         for i in line:
             if not re.match(r"[A-Z][a-z]?", i):
                 is_vasp_5 = False
-                break
+                break  
         if is_vasp_5:
             text_types = deepcopy(line)
-            if types is not None and not set(text_types).issubset(set(types)):
-                raise error.ValueError("Unknown species in poscar: {0} not in {1}."
-                                       .format(text_types, types))
+
+            if types is not None:
+                # If partial H (e.g. H.75) are specified they will replace H
+                # in the order they were specified
+                for s in types:
+                    if s[0] == "H" and "." in s:
+                        for i,c in enumerate(text_types):
+                            if c == "H":
+                                text_types[i] = s
+                                break
+
+                if not set(text_types).issubset(set(types)):
+                    raise error.ValueError("Unknown species in poscar: {0} not in {1}."
+                                           .format(text_types, types))
             types = text_types
             line = poscar.readline().split()
         if types is None:

src/pylada/decorations/CMakeLists.txt

@@ -21,10 +21,11 @@
 ###############################
 
 add_cython_target(decorations_cython _decorations.pyx CXX ${PYVER}})
-add_library(_decorations MODULE ${decorations_cython})
-python_extension_module(_decorations)
+
+python_add_library(_decorations MODULE WITH_SOABI ${decorations_cython})
+
 target_include_directories(
     _decorations PRIVATE
-    "${PY_HEADER_DIR}" "${NumPy_INCLUDE_DIRS}" "${PYTHON_INCLUDE_DIR}"
+    "${PY_HEADER_DIR}" "${NumPy_INCLUDE_DIRS}" "${Python3_INCLUDE_DIR}"
 )
 install(TARGETS _decorations LIBRARY DESTINATION ${PY_ROOT_DIR}/decorations)

src/pylada/ewald/CMakeLists.txt

@@ -20,8 +20,8 @@
 #  <http://www.gnu.org/licenses/>.
 ###############################
 add_cython_target(ewald_cython ewald.pyx CXX ${PYVER})
-add_library(ewald MODULE
-    ${ewald_cython} ep_com.f90 erfc.cc ewaldf.f90 ewald.cc)
-python_extension_module(ewald)
-target_include_directories(ewald PRIVATE "${PY_HEADER_DIR}")
+
+python_add_library(ewald MODULE WITH_SOABI ${ewald_cython} ep_com.f90 erfc.cc ewaldf.f90 ewald.cc)
+
+target_include_directories(ewald PRIVATE "${PY_HEADER_DIR}" "${Python3_INCLUDE_DIR}")
 install(TARGETS ewald LIBRARY DESTINATION ${PY_ROOT_DIR}/ewald)

src/pylada/ipython/launch/__init__.py

@@ -167,7 +167,6 @@ def mppalloc(job):
             """ Returns number of processes for this job. """
             natom = len(job.structure)  # number of atoms.
             # Round down to a multiple of ppn
-#vladan            nnode = max(1, natom / event.ppn)
             nnode = max(1, natom // event.ppn)
             nproc = nnode * event.ppn
             return nproc