diff --git a/CMakeLists.txt b/CMakeLists.txt
index 612387e4..56cc1e1f 100644
--- a/CMakeLists.txt
+++ b/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)
diff --git a/cmake_modules/eigen/CMakeLists.txt b/cmake_modules/eigen/CMakeLists.txt
index 391edf08..df6bcfcf 100644
--- a/cmake_modules/eigen/CMakeLists.txt
+++ b/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)
diff --git a/src/pylada/CMakeLists.txt b/src/pylada/CMakeLists.txt
index 7bccf299..a72e0423 100644
--- a/src/pylada/CMakeLists.txt
+++ b/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)
diff --git a/src/pylada/crystal/CMakeLists.txt b/src/pylada/crystal/CMakeLists.txt
index e704e5f5..20282d88 100644
--- a/src/pylada/crystal/CMakeLists.txt
+++ b/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)
diff --git a/src/pylada/crystal/_coordination_shells.pyx b/src/pylada/crystal/_coordination_shells.pyx
index 19f84dd9..cff1dff3 100644
--- a/src/pylada/crystal/_coordination_shells.pyx
+++ b/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])
diff --git a/src/pylada/crystal/_primitive.pyx b/src/pylada/crystal/_primitive.pyx
index 92f094c1..832a0de2 100644
--- a/src/pylada/crystal/_primitive.pyx
+++ b/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
diff --git a/src/pylada/crystal/_space_group.pyx b/src/pylada/crystal/_space_group.pyx
index 91a340f0..6101f3e4 100644
--- a/src/pylada/crystal/_space_group.pyx
+++ b/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 """
diff --git a/src/pylada/crystal/cutilities.pyx b/src/pylada/crystal/cutilities.pyx
index 920cb988..261b6c49 100644
--- a/src/pylada/crystal/cutilities.pyx
+++ b/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
diff --git a/src/pylada/crystal/defects/CMakeLists.txt b/src/pylada/crystal/defects/CMakeLists.txt
index 0c54f549..c5ddea6a 100644
--- a/src/pylada/crystal/defects/CMakeLists.txt
+++ b/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)
diff --git a/src/pylada/crystal/read.py b/src/pylada/crystal/read.py
index d3066128..0796ed58 100644
--- a/src/pylada/crystal/read.py
+++ b/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:
diff --git a/src/pylada/decorations/CMakeLists.txt b/src/pylada/decorations/CMakeLists.txt
index 06d3124a..2dee0bf6 100644
--- a/src/pylada/decorations/CMakeLists.txt
+++ b/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)
diff --git a/src/pylada/ewald/CMakeLists.txt b/src/pylada/ewald/CMakeLists.txt
index f07a93e5..831588f9 100644
--- a/src/pylada/ewald/CMakeLists.txt
+++ b/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)
diff --git a/src/pylada/ipython/launch/__init__.py b/src/pylada/ipython/launch/__init__.py
index 24d8049f..d77ac9a9 100644
--- a/src/pylada/ipython/launch/__init__.py
+++ b/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
diff --git a/src/pylada/vasp/extract/base.py b/src/pylada/vasp/extract/base.py
index ead78ca5..4e497c75 100644
--- a/src/pylada/vasp/extract/base.py
+++ b/src/pylada/vasp/extract/base.py
@@ -269,12 +269,12 @@ def structure(self):
         """ Greps structure and total energy from OUTCAR. """
         if self.nsw == 0 or self.ibrion == -1:
             return self.initial_structure
-
+        
         try:
             result = self._contcar_structure
         except:
             result = self._grep_structure
-
+            
         # tries to find adequate name for structure.
         try:
             name = self.system
@@ -329,6 +329,7 @@ def structure(self):
         else:
             for force, atom in zip(forces, result):
                 atom.force = force
+        
         return result
 
     @property
@@ -447,7 +448,7 @@ def _contcar_structure(self):
         """ Greps structure from CONTCAR. """
         from ...crystal import read
         from quantities import eV
-
+        
         result = read.poscar(self.__contcar__(), self.species)
         result.energy = float(self.total_energy.rescale(eV)) if self.is_dft else 0e0
         return result
@@ -470,8 +471,21 @@ def ions_per_specie(self):
     @property
     @make_cached
     def species(self):
+        import re
         """ Greps species from OUTCAR. """
-        return [ u.group(1) for u in self._search_OUTCAR(r"""VRHFIN\s*=\s*(\S+)\s*:""") ]
+        result = []
+        for u in self._search_OUTCAR(r"""VRHFIN\s*=\s*(.*)$"""):
+            atom = u[1].split(":")[0].strip()
+            match = re.search(r"Z\s*=\s*(\S+)", u[1])
+            if match is None:
+                result.append(atom)
+            else:
+                result.append(atom + ("%3.2f"%(eval(match[1]))).strip(" 0"))
+        
+        # return [ u.group(1) for u in self._search_OUTCAR(r"""VRHFIN\s*=\s*(\S+)\s*:""") ]
+        # return [ u.group(1) for u in self._search_OUTCAR(r"""TITEL\s*=\s*\S*\s*(\S+)""") ]
+
+        return result
 
     @property
     @make_cached
diff --git a/src/pylada/vasp/functional.py b/src/pylada/vasp/functional.py
index 9fdeef70..13c1bddb 100644
--- a/src/pylada/vasp/functional.py
+++ b/src/pylada/vasp/functional.py
@@ -304,7 +304,7 @@ def __init__(self, copy=None, species=None, kpoints=None, **kwargs):
 
             .. _CONTCAR: http://cms.mpi.univie.ac.at/vasp/guide/node60.html
         """
-        self.isym = ChoiceKeyword(values=list(range(-1, 3)))
+        self.isym = ChoiceKeyword(values=list(range(-1, 4)))
         """ Symmetry scheme.
 
             .. seealso:: ISYM_
@@ -1017,8 +1017,9 @@ def bringup(self, structure, outdir, **kwargs):
             # creates POTCAR file
             logger.debug("vasp/functional bringup: files.POTCAR: %s " % files.POTCAR)
             with open(files.POTCAR, 'w') as potcar:
-                self.write_potcar(potcar, structure)
-
+                for s in specieset(structure):
+                    outLines = self.species[s].read_potcar()
+                    potcar.writelines(outLines)
             # Add is running file marker.
             local_path(outdir).join('.pylada_is_running').ensure(file=True)
             self._copy_additional_files(outdir)
@@ -1114,14 +1115,7 @@ def write_kpoints(self, file, structure, kpoints=None):
                 outLine = "{0[0]} {0[1]} {0[2]} {1}\n".format(
                     kpoint, 1 if len(kpoint) == 3 else kpoint[3])
                 file.write(outLine)
-    
-    def write_potcar(self, file, structure):
-        """ Writes the potcar file """
-        from ..crystal import specieset
-        for s in specieset(structure):
-            outLines = self.species[s].read_potcar()
-            file.writelines(outLines)
-    
+
     def __repr__(self, defaults=True, name=None):
         """ Returns representation of this instance """
         from ..tools.uirepr import uirepr
diff --git a/src/pylada/vasp/keywords.py b/src/pylada/vasp/keywords.py
index 1d596769..8998b45f 100644
--- a/src/pylada/vasp/keywords.py
+++ b/src/pylada/vasp/keywords.py
@@ -356,7 +356,10 @@ def value(self, value):
         if value is None:
             self._value = None
             return None
-        from pylada import is_vasp_4
+        try:
+            from pylada import is_vasp_4
+        except:
+            is_vasp_4 = False
         if not hasattr(value, 'lower'):
             raise TypeError("ALGO cannot be set with {0}.".format(value))
         lower = value.lower().rstrip().lstrip()
@@ -871,8 +874,7 @@ def output_map(self, **kwargs):
         from ..error import ValueError
         from ..crystal import write, read, specieset
         from . import files
-        from pylada import is_vasp_4
-        
+
         istruc = self._value
         if istruc is None:
             istruc = 0
@@ -886,7 +888,7 @@ def output_map(self, **kwargs):
             has_restart = vasp.restart.success
         if has_restart:
             structure = vasp.restart.structure
-
+            
         # determines which CONTCAR is the latest, if any exist.
         if istruc in [-1, 1]:
             last_contcar = latest_file(join(outdir, files.CONTCAR))
@@ -896,6 +898,7 @@ def output_map(self, **kwargs):
             # specifically, the order of the atoms of a given specie are the same).
             if last_contcar is not None:
                 other = read.poscar(path=last_contcar, types=specieset(structure))
+                
                 if len(other) != len(structure):
                     raise ValueError('CONTCAR and input structure differ in size.')
                 for type in specieset(other):
@@ -908,6 +911,9 @@ def output_map(self, **kwargs):
                         if a.type != type:
                             continue
                         a.pos = b.pos
+                        if len(getattr(b, 'freeze', '')) != 0:
+                            a.freeze = b.freeze
+    
                 structure.cell = other.cell
                 structure.scale = other.scale
                 print("setting cell and scale!")
@@ -920,7 +926,8 @@ def output_map(self, **kwargs):
             raise ValueError('Structure scale is zero')
         if structure.volume < 1e-8:
             raise ValueError('Structure volume is zero')
-        write.poscar(structure, join(outdir, 'POSCAR'), vasp5=not is_vasp_4)
+        
+        write.poscar(structure, join(outdir, 'POSCAR'))
         return None