Merge pull request #1 from pierre-24/manipulate_geometry

Manipulate geometries

DOCUMENTATION.md

@@ -3,6 +3,7 @@
 1. [Purpose](#purpose)
 2. [Installation](#installation)
 3. [Usage](#usage)
+4. [Geometry manipulation tools](#geometry-manipulation-tools)
 4. [Contribute](#contribute)
 5. [Who?](#who)
 
@@ -80,16 +81,8 @@ You can get info about a slab by runnning `ei-check-slab`:
 ei-check-slab POSCAR
 ```
 Among others, the interslab distance (i.e., the vacuum between two repetition of the slab) should be adjusted (see [10.1021/acs.jctc.5b00170](https://dx.doi.org/10.1021/acs.jctc.5b00170)).
-
-To adjust the interslab distance, you can use `ei-set-vacuum`, which creates a new geometry while enforcing vacuum size (i.e., the size of the last lattice vector).
-For example, to adjust the vacuum size to 25 Å:
-```bash
-mv POSCAR POSCAR_old
-ei-set-vacuum POSCAR_old -v 25.0 -o POSCAR
-```
-The new geometry is saved in `POSCAR`.
-There should be enough vacuum to get an accurate value for the reference (vacuum) potential, which is used to compute the work function.
-Note that the slab is z-centered by the procedure.
+See [below](#geometry-manipulation-tools) for a tool to do so.
+Note that there should be enough vacuum to get an accurate value for the reference (vacuum) potential, which is used to compute the work function.
 
 You can get the number of electron that your system contains with `ei-get-nzc`:
 ```bash
@@ -173,7 +166,7 @@ And then another dataset, with:
 
 Please refer to [10.1039/c9cp06684e](https://doi.org/10.1021/10.1039/c9cp06684e) (and reference therein) for different information that you can extract from those data, such as the surface capacitances, the fukui functions, etc.
 
-### 4. Example
+### 5. Example
 
 See [this archive](https://drive.google.com/file/d/1TkLHsbzXJz_slb6X06r1NbkHko73-fin/view?usp=sharing), which contains an example of calculation on a Li (100) slab using the PBM approach, inspired by [10.1021/acs.jctc.1c01237](https://doi.org/10.1021/acs.jctc.1c01237).
 It gives the following curve:
@@ -183,6 +176,33 @@ It gives the following curve:
 A capacitance of 0.0535 e/V is extracted from this curve using its second derivative.
 Due to the strong anharmonicity (the curve is clearly not symmetric around PZC), the actual value should be a little smaller.
 
+## Geometry manipulation tools
+
+It is generally recommended to use tools such as [ASE](https://wiki.fysik.dtu.dk/ase/) or [pymatgen](https://pymatgen.org/) to manipulate the geometries.
+However, `ec-interface` comes with a few handy tools to perform routine operations.
+Use them with care!
+
++ To adjust the interslab distance, you can use `ei-set-vacuum`, which creates a new geometry while enforcing vacuum size (i.e., the size of the last lattice vector).
+  For example, to adjust the vacuum size to 25 Å:
+  ```bash
+  mv POSCAR POSCAR_old
+  ei-set-vacuum POSCAR_old -v 25.0 -o POSCAR
+  ```
+  The new geometry is saved in `POSCAR`. Note that the slab is z-centered by the procedure.
++ To turn an XYZ geometry into POSCAR, you can use `ei-to-vasp-geometry`:
+  ```bash
+  ei-to-vasp-geometry molecule.xyz --lattice=10,10,10 -o POSCAR
+  ```
+  It comes with a few options, such as `--lattice` to set the lattice vectors and `--sort` to group atoms types (so that it is easier to create the POTCAR).
++ To merge two POSCARs, you can use `ei-merge-poscar`:
+  ```bash
+  ei-to-vasp-geometry POSCAR_cell POSCAR_substrate --shift=5,5,7 -o POSCAR
+  ```
+  It allows to merge two geometries. 
+  The lattice of the first geometry (here `POSCAR_cell`) is used in the final one.
+  The `--shift` option allows to reposition the second molecule in the first.
+  Note that `Selective dynamics` information are kept.
+
 ## Contribute
 
 Contributions, either with [issues](https://github.com/pierre-24/ec-interface/issues) or [pull requests](https://github.com/pierre-24/ec-interface/pulls) are welcomed.

ec_interface/molecular_geometry.py

@@ -0,0 +1,80 @@
+import numpy
+from typing import TextIO, List
+from numpy.typing import NDArray
+
+from ec_interface import vasp_geometry
+
+
+class MolecularGeometry:
+    """Molecular geometry, i.e., a geometry without any PBC
+    """
+
+    def __init__(self, symbols: List[str], positions: NDArray[float]):
+        assert positions.shape == (len(symbols), 3)
+
+        self.symbols = symbols
+        self.positions = positions
+
+    def __len__(self) -> int:
+        return len(self.symbols)
+
+    @classmethod
+    def from_xyz(cls, f: TextIO) -> 'MolecularGeometry':
+        """Read geometry from a XYZ file
+        """
+
+        symbols = []
+        positions = []
+
+        n = int(f.readline())
+        f.readline()
+
+        for i in range(n):
+            data = f.readline().split()
+            symbols.append(data[0])
+            positions.append([float(x) for x in data[1:]])
+
+        return cls(symbols, numpy.array(positions))
+
+    def as_xyz(self, title: str = '') -> str:
+        """Get XYZ representation of this geometry"""
+
+        r = '{}\n{}'.format(len(self), title)
+        for i in range(len(self)):
+            r += '\n{:2} {: .7f} {: .7f} {: .7f}'.format(self.symbols[i], *self.positions[i])
+
+        return r
+
+    def to_vasp(
+            self,
+            lattice_vectors: NDArray,
+            title: str = '',
+            sort: bool = False,
+            shift_positions: NDArray = None
+    ) -> vasp_geometry.Geometry:
+        """Convert to a VASP geometry"""
+
+        symbols = self.symbols.copy()
+        positions = self.positions.copy()
+
+        # sort if any
+        if sort:
+            sorted_symbols = sorted(zip(symbols, range(len(self))))
+            symbols = [x[0] for x in sorted_symbols]
+            positions[list(x[1] for x in sorted_symbols)] = positions
+
+        # count ions
+        ions_types, ions_numbers = vasp_geometry.count_ions(symbols)
+
+        # shift if any
+        if shift_positions is not None:
+            positions += shift_positions
+
+        return vasp_geometry.Geometry(
+            title=title,
+            lattice_vectors=lattice_vectors,
+            ion_types=ions_types,
+            ion_numbers=ions_numbers,
+            positions=positions,
+            is_direct=False
+        )

ec_interface/scripts/__init__.py

@@ -1,6 +1,9 @@
 import argparse
 import pathlib
 
+import numpy
+from numpy._typing import NDArray
+
 from ec_interface.ec_parameters import ECParameters
 
 INPUT_NAME = 'ec_interface.yml'
@@ -29,3 +32,29 @@ def get_ec_parameters(fp: str) -> ECParameters:
 
     with p.open() as f:
         return ECParameters.from_yaml(f)
+
+
+def get_lattice(inp: str) -> NDArray:
+    elmts = inp.split(',')
+    if len(elmts) != 3:
+        raise argparse.ArgumentTypeError('Lattice must have three elements')
+
+    lattice = numpy.zeros((3, 3))
+    for i in range(3):
+        try:
+            lattice[i, i] = float(elmts[i])
+        except ValueError:
+            raise argparse.ArgumentTypeError('Element {} of lattice is not a float'.format(i))
+
+    return lattice
+
+
+def get_vec(inp: str) -> NDArray:
+    elmts = inp.split(',')
+    if len(elmts) != 3:
+        raise argparse.ArgumentTypeError('COM shift must have three elements')
+
+    try:
+        return numpy.array([float(x) for x in elmts])
+    except ValueError:
+        raise argparse.ArgumentTypeError('COM shift must be 3 floats')

ec_interface/scripts/merge_poscar.py

@@ -0,0 +1,36 @@
+"""
+Merge two POSCAR geometries into a single.
+The lattice vector of the first geometry is used in the final one.
+"""
+
+import argparse
+import sys
+
+from ec_interface.scripts import get_vec
+from ec_interface.vasp_geometry import Geometry
+
+
+def get_arguments_parser():
+    parser = argparse.ArgumentParser(description=__doc__)
+
+    parser.add_argument('infile', help='source', type=argparse.FileType('r'))
+    parser.add_argument('additional', help='additional', type=argparse.FileType('r'))
+
+    parser.add_argument('-s', '--shift', help='Shift positions', type=get_vec, default='0,0,0')
+
+    parser.add_argument('-o', '--poscar', type=argparse.FileType('w'), default=sys.stdout)
+    parser.add_argument('-C', '--cartesian', action='store_true', help='Output in cartesian coordinates')
+
+    return parser
+
+
+def main():
+    args = get_arguments_parser().parse_args()
+
+    geometry = Geometry.from_poscar(args.infile)
+    additional = Geometry.from_poscar(args.additional)
+    geometry.merge_with(additional, shift=args.shift).to_poscar(args.poscar, direct=not args.cartesian)
+
+
+if __name__ == '__main__':
+    main()

ec_interface/scripts/to_vasp_geometry.py

@@ -0,0 +1,42 @@
+"""
+Convert a XYZ file to a VASP geometry
+"""
+
+import argparse
+import sys
+import numpy
+
+from ec_interface.molecular_geometry import MolecularGeometry
+from ec_interface.scripts import get_lattice, get_vec
+
+
+def get_arguments_parser():
+    parser = argparse.ArgumentParser(description=__doc__)
+
+    parser.add_argument('infile', help='source interface', type=argparse.FileType('r'))
+
+    parser.add_argument('-l', '--lattice', help='lattice vector size', type=get_lattice, default='10,10,10')
+    parser.add_argument('-s', '--shift', help='Shift positions', type=get_vec, default='0,0,0')
+    parser.add_argument('--sort', help='sort atoms', action='store_true')
+
+    parser.add_argument('-o', '--poscar', type=argparse.FileType('w'), default=sys.stdout)
+    parser.add_argument('-C', '--cartesian', action='store_true', help='Output in cartesian coordinates')
+    parser.add_argument('-S', '--selective', action='store_true', help='Use selective dynamics in output')
+
+    return parser
+
+
+def main():
+    args = get_arguments_parser().parse_args()
+
+    geometry = MolecularGeometry.from_xyz(args.infile)
+    new_geometry = geometry.to_vasp(lattice_vectors=args.lattice, sort=args.sort, shift_positions=args.shift)
+
+    if args.selective:
+        new_geometry.selective_dynamics = numpy.ones((len(new_geometry), 3), dtype=bool)
+
+    new_geometry.to_poscar(args.poscar, direct=not args.cartesian)
+
+
+if __name__ == '__main__':
+    main()

ec_interface/vasp_geometry.py

@@ -1,6 +1,7 @@
 import numpy
+from numpy.typing import NDArray
 
-from typing import List, TextIO, Dict
+from typing import List, TextIO, Dict, Tuple, Optional
 
 
 def get_zvals(f: TextIO) -> Dict[str, float]:
@@ -23,6 +24,32 @@ def get_zvals(f: TextIO) -> Dict[str, float]:
     return zvals
 
 
+def count_ions(symbols: List[str]) -> Tuple[List[str], List[int]]:
+    """Count ions in VASP style (i.e., get a pair `ion_types, ion_numbers` from `symbols`)
+    """
+
+    ions_types = []
+    ions_numbers = []
+
+    if len(symbols) > 0:
+        current_ion = symbols[0]
+        current_number = 1
+        for symbol in symbols[1:]:
+            if symbol == current_ion:
+                current_number += 1
+            else:
+                ions_types.append(current_ion)
+                ions_numbers.append(current_number)
+                current_ion = symbol
+                current_number = 1
+
+        # add last:
+        ions_types.append(current_ion)
+        ions_numbers.append(current_number)
+
+    return ions_types, ions_numbers
+
+
 class Geometry:
     def __init__(
         self,
@@ -213,3 +240,37 @@ def nelect(self, f: TextIO) -> float:
             nelect += n * zvals[symbol]
 
         return nelect
+
+    def merge_with(self, other: 'Geometry', title: str = '', shift: Optional[NDArray] = None) -> 'Geometry':
+
+        # shift `additional` if any
+        if shift is not None:
+            c = other._cartesian_coordinates + shift
+        else:
+            c = other._cartesian_coordinates
+
+        # merge position, ion types and numbers.
+        positions = numpy.vstack([self._cartesian_coordinates, c])
+        ion_types = self.ion_types.copy() + other.ion_types.copy()
+        ion_numbers = self.ion_numbers.copy() + other.ion_numbers.copy()
+
+        selective_dynamics = None
+
+        # keep selective dynamic if any
+        if self.selective_dynamics is not None or other.selective_dynamics is not None:
+            selective_dynamics = numpy.vstack([
+                self.selective_dynamics
+                if self.selective_dynamics is not None else numpy.ones((len(self), 3), dtype=bool),
+                other.selective_dynamics
+                if other.selective_dynamics is not None else numpy.ones((len(other), 3), dtype=bool)
+            ])
+
+        return Geometry(
+            title,
+            lattice_vectors=self.lattice_vectors.copy(),
+            ion_types=ion_types,
+            ion_numbers=ion_numbers,
+            positions=positions,
+            is_direct=False,
+            selective_dynamics=selective_dynamics
+        )

pyproject.toml

@@ -55,6 +55,8 @@ dev = [
 'ei-make-directories' = 'ec_interface.scripts.make_directories:main'
 'ei-extract-data' = 'ec_interface.scripts.extract_data:main'
 'ei-compute-fee' = 'ec_interface.scripts.compute_fee:main'
+'ei-merge-poscar' = 'ec_interface.scripts.merge_poscar:main'
+'ei-to-vasp-geometry' = 'ec_interface.scripts.to_vasp_geometry:main'
 'ei-xy-average' = 'ec_interface.scripts.make_xy_average:main'
 
 [tool.setuptools]

tests/THF.xyz

@@ -0,0 +1,15 @@
+13
+C4H8O
+C        0.729192000      0.990407000     -0.229613000
+C        1.162040000     -0.425591000      0.132806000
+C       -0.729191000      0.990408000      0.229613000
+H       -1.519126000     -0.476343000     -1.169724000
+H       -1.944522000     -0.820810000      0.520504000
+C       -1.162040000     -0.425589000     -0.132806000
+H       -0.788206000      1.144203000      1.311812000
+H       -1.340555000      1.751411000     -0.259436000
+O       -0.000001000     -1.246840000      0.000000000
+H        1.944521000     -0.820813000     -0.520504000
+H        1.519125000     -0.476346000      1.169724000
+H        1.340558000      1.751409000      0.259435000
+H        0.788207000      1.144201000     -1.311812000