Initial Work on Reading Crystal and Creation Gaussian Input

crystalpol/__main__.py

@@ -0,0 +1,62 @@
+from crystalpol.polarization import Polarization
+from crystalpol.shared.config import Config
+
+from yaml.loader import SafeLoader
+import yaml
+
+import setproctitle
+import argparse
+import os
+
+
+__VERSION = "v0.0.1"
+os.nice(+19)
+setproctitle.setproctitle("crystalpol-{}".format(__VERSION))
+
+
+def main():
+    """
+    Read and store the arguments passed to the program
+    and set the usage and help messages.
+    """
+
+    parser = argparse.ArgumentParser(prog="CrystalPol")
+    parser.add_argument(
+        "-v", "--version", action="version", version=f"crystalpol-{__VERSION}"
+    )
+    parser.add_argument(
+        "-c", "--config",
+        dest="config",
+        default="config.yml",
+        metavar="INFILE",
+        help="Config file of crystalpol [default = config.yml]"
+    )
+    parser.add_argument(
+        "-i", "--input",
+        dest="infile",
+        default="crystal.xyz",
+        metavar="INFILE",
+        help="Input file of crystalpol [default = crystal.xyz]"
+    )
+    parser.add_argument(
+        "-o", "--output",
+        dest="outfile",
+        default="run.log",
+        metavar="OUTFILE",
+        help="Output file of crystalpol [default = run.log]"
+    )
+    args = parser.parse_args()
+
+    try:
+        with open(args.config) as file:
+            data = yaml.load(file, Loader=SafeLoader)
+            config = Config(**data.get('crystal_pol'))
+    except IOError:
+        raise RuntimeError('Invalid or Missing Config File.')
+
+    pol = Polarization(args.infile, args.outfile, config)
+    pol.run()
+
+
+if __name__ == "__main__":
+    main()

crystalpol/gaussian.py

@@ -0,0 +1,114 @@
+from crystalpol.shared.utils.ptable import atom_symbol
+from crystalpol.shared.system.crystal import Crystal
+from crystalpol.shared.config import Config
+
+from pathlib import Path, PosixPath
+from typing import TextIO
+import subprocess
+import textwrap
+import shutil
+import os
+
+
+class Gaussian:
+    __slots__ = ('qmprog', 'config', 'keywords')
+
+    def __init__(self, config: Config) -> None:
+        self.qmprog = "g16"
+        self.config = config
+
+        self.check_keyword()
+
+    def check_keyword(self):
+        if self.config.pop not in ["chelpg", "mk", "nbo"]:
+            self.config.pop = "chelpg"
+
+    def run(self, cycle: int, crystal: Crystal) -> None:
+
+        file = Path("simfiles", f"crystal-{str(cycle).zfill(2)}.gjf")
+
+        self.make_gaussian_input(cycle, file, crystal)
+
+        if shutil.which("bash") is not None:
+            exit_status = subprocess.call(
+                [
+                    "bash",
+                    "-c",
+                    "exec -a {}-step{} {} {}".format(
+                        self.qmprog, cycle, self.qmprog, file.name
+                    ),
+                ]
+            )
+        else:
+            exit_status = subprocess.call([self.qmprog, file.name])
+
+        if exit_status != 0:
+            raise RuntimeError("Gaussian process did not exit properly")
+
+        return self.read_charges_from_gaussian_output()
+
+    def create_simulation_dir(self):
+        if not os.path.exists(self.config.simulation_dir):
+            os.makedirs(self.config.simulation_dir)
+        else:
+            raise RuntimeError(
+                f"Simulation directory '{self.config.simulation_dir}' already exists. "
+                f"Please remove it before proceeding."
+            )
+
+    def make_gaussian_input(self, cycle: int, file: PosixPath, crystal: Crystal) -> str:
+
+        with open(file, 'w+') as fh:
+
+            fh.write(f"%Mem={self.config.mem}MB\n")
+
+            fh.write(f"%Nprocs={self.config.n_procs}\n")
+
+            kwords_line = f"#P {self.config.level} " \
+                          f"Pop = {self.config.pop} " \
+                          f"Density = Current " \
+                          f"NoSymm "
+
+            if cycle > 1:
+                kwords_line += f"charge"
+
+            fh.write(textwrap.fill(kwords_line, 90))
+            fh.write("\n")
+
+            fh.write(f"\n{self.config.comment} - Cycle number {cycle}\n")
+            fh.write("\n")
+            fh.write(f"{self.config.mult[0]}, {self.config.mult[1]}\n")
+
+            for atom in crystal[0][0]:
+                symbol = atom_symbol[atom.na]
+                fh.write(
+                    f"{symbol:<2s}    "
+                    f"{atom.rx:>10.5f}   "
+                    f"{atom.ry:>10.5f}   "
+                    f"{atom.rz:>10.5f}\n"
+                )
+
+            if cycle > 1:
+                self.make_gaussian_charges(fh, crystal)
+
+            fh.seek(0)
+            return fh.read()
+
+    def make_gaussian_charges(self, fh: TextIO, crystal: Crystal) -> None:
+
+        fh.write("\n")
+
+        for index_cell, cell in enumerate(crystal):
+            for index_mol, molecule in enumerate(cell):
+                if (index_cell == 0 and index_mol != 0) or (index_cell != 0):
+                    for atom in molecule:
+                        symbol = atom_symbol[atom.na]
+                        fh.write(
+                            f"{symbol:<2s}    "
+                            f"{atom.rx:>10.5f}   "
+                            f"{atom.ry:>10.5f}   "
+                            f"{atom.rz:>10.5f}\n"
+                        )
+
+    def read_charges_from_gaussian_output(self) -> None:
+        pass

crystalpol/polarization.py

@@ -0,0 +1,74 @@
+import os
+
+from crystalpol.shared.system.molecule import Molecule
+from crystalpol.shared.system.crystal import Crystal
+from crystalpol.shared.system.atom import Atom
+from crystalpol.shared.config import Config
+from crystalpol.gaussian import Gaussian
+
+from typing import List
+
+
+class Polarization:
+    __slots__ = ('geom_file', 'outfile', 'config', 'gaussian', 'crystal')
+
+    def __init__(self, geom_file: str, outfile: str, config: Config) -> None:
+        self.crystal = None
+
+        self.geom_file = geom_file
+        self.outfile = outfile
+        self.config = config
+
+        self.gaussian = Gaussian(config)
+
+    def run(self):
+
+        self.read_crystal()
+
+        self.gaussian.run(1, self.crystal)
+
+    def read_crystal(self) -> None:
+        with open(self.geom_file, 'r') as geom_file:
+            lines = geom_file.readlines()
+
+        molecules = self._get_molecules_from_lines(lines)
+        structure = self._get_crystal_structure(molecules[0])
+
+        self.crystal = Crystal([structure])
+        for molecule in molecules:
+            self.crystal.add_cell([molecule])
+
+    def _get_molecules_from_lines(self, lines: List[str]) -> List[Molecule]:
+        if (len(lines) % self.config.n_atoms) == 0:
+            molecules: List[Molecule] = []
+
+            for index, molecule in enumerate(split(lines, self.config.n_atoms)):
+                mol = Molecule(f"Molecule-{index}")
+                for atom_line in molecule:
+                    symbol, rx, ry, rz = tuple(atom_line.split())
+                    mol.add_atom(
+                        Atom(
+                            rx, ry, rz,
+                            symbol=symbol.ljust(2),
+                        )
+                    )
+                molecules.append(mol)
+
+            return molecules
+
+        else:
+            raise RuntimeError(
+                "Invalid Geom File, the number of atoms doesn't match the number of lines."
+            )
+
+    @staticmethod
+    def _get_crystal_structure(molecule: Molecule) -> List[str]:
+        structure: List[str] = []
+        for atom in molecule:
+            structure.append(atom.symbol)
+        return structure
+
+
+def split(array: List, partitions: int):
+    for i in range(0, len(array), partitions):
+        yield array[i: i + partitions]

crystalpol/shared/config/__init__.py

@@ -0,0 +1 @@
+from .config import Config

crystalpol/shared/config/config.py

@@ -0,0 +1,42 @@
+from dataclasses import dataclass, field, fields
+
+
+@dataclass
+class Config:
+
+    mem: int
+    level: str
+    n_atoms: int
+
+    n_procs: int = 1
+    pop: str = "chelpg"
+    comment: str = "crystalpol"
+    simulation_dir = "simfiles"
+    mult: list = \
+        field(default_factory=lambda: [0, 1])
+
+    def __post_init__(self):
+        for field in fields(self):
+            value = getattr(self, field.name)
+            if not isinstance(value, field.type):
+                raise ValueError(
+                    f'Expected {field.name} to be {field.type}, '
+                    f'got {repr(value)}'
+                )
+
+        if self.mem is None or self.mem <= 0:
+            raise ValueError(
+                f'Invalid value for mem: {self.mem},'
+                f'Memory must be a integer greater than 0.'
+            )
+
+        if self.level is None:
+            raise ValueError(
+                f'Invalid value for level. Level must not be none.'
+            )
+
+        if self.n_atoms is None or self.n_atoms <= 0:
+            raise ValueError(
+                f'Invalid value for n_atoms: {self.mem},'
+                f'Number of Atoms must be a integer greater than 0.'
+            )

crystalpol/shared/system/atom.py

@@ -0,0 +1,46 @@
+from crystalpol.shared.utils.ptable import atom_mass, atom_symbol
+
+
+class Atom:
+    """
+    Atom class declaration. This class is used throughout the DicePlayer program to represent atoms.
+    Attributes:
+        na (int): Atomic number of the represented atom.
+        symbol (str): Atomic symbol of the represented atom.
+        rx (float): x cartesian coordinates of the represented atom.
+        ry (float): y cartesian coordinates of the represented atom.
+        rz (float): z cartesian coordinates of the represented atom.
+    """
+
+    def __init__(
+            self,
+            rx: float,
+            ry: float,
+            rz: float,
+            na: int = None,
+            symbol: str = None,
+
+    ) -> None:
+        """
+        The constructor function __init__ is used to create new instances of the Atom class.
+        Args:
+            na (int): Atomic number of the represented atom.
+            symbol (str): Atomic symbol of the represented atom.
+            rx (float): x cartesian coordinates of the represented atom.
+            ry (float): y cartesian coordinates of the represented atom.
+            rz (float): z cartesian coordinates of the represented atom.
+        """
+
+        if na is not None:
+            self.na = na
+            self.symbol = atom_symbol[self.na]
+
+        if symbol is not None and symbol in atom_symbol:
+            self.symbol = symbol
+            self.na = atom_symbol.index(self.symbol)
+
+        self.rx = rx
+        self.ry = ry
+        self.rz = rz
+        self.chg = None
+        self.mass = atom_mass[self.na]

crystalpol/shared/system/crystal.py

@@ -0,0 +1,43 @@
+from crystalpol.shared.system.molecule import Molecule
+
+from typing import List
+
+
+class Crystal:
+    """
+    This class represents a crystal, it will be organized in a list structure.
+    Each element is a unitary cell in the crystal. And this unitary cell will have
+    Molecules in it.
+    """
+
+    def __init__(self, structure: List[List[str]]):
+        self.structure = structure
+
+        self.cells = []
+
+    def __iter__(self):
+        for cell in self.cells:
+            yield cell
+
+    def __len__(self):
+        return len(self.cells)
+
+    def __getitem__(self, index):
+        return self.cells[index]
+
+    def add_cell(self, cell: List[Molecule]) -> None:
+        valid = self._is_valid_cell(cell)
+        if not valid:
+            raise ValueError(
+                "This cell does not obey the declared format for this Crystal."
+            )
+        else:
+            self.cells.append(cell)
+
+    def _is_valid_cell(self, cell: List[Molecule]) -> bool:
+        if len(cell) == len(self.structure):
+            for i, molecule in enumerate(cell):
+                if len(molecule.atoms) == len(self.structure[i]) \
+                        and all(atom.symbol == self.structure[i][j] for j, atom in enumerate(molecule.atoms)):
+                    return True
+        return False

crystalpol/shared/system/molecule.py

@@ -0,0 +1,88 @@
+from crystalpol.shared.system.atom import Atom
+from crystalpol.shared.utils.ptable import atom_symbol
+
+from typing import Final, List, TextIO
+import math
+
+
+""" Constants of unit conversion """
+BOHR2ANG: Final[float] = 0.52917721092
+ANG2BOHR: Final[float] = 1 / BOHR2ANG
+
+
+class Molecule:
+    """
+    Molecule class declaration. This class is used throughout the DicePlayer program to represent molecules.
+
+    Atributes:
+        molname (str): The name of the represented molecule
+        atom (List[Atom]): List of atoms of the represented molecule
+        position (NDArray[Any, Any]): The position relative to the internal atoms of the represented molecule
+        energy (NDArray[Any, Any]): The energy of the represented molecule
+        gradient (NDArray[Any, Any]): The first derivative of the energy relative to the position
+        hessian (NDArray[Any, Any]): The second derivative of the energy relative to the position
+        total_mass (int): The total mass of the molecule
+        com (NDArray[Any, Any]): The center of mass of the molecule
+    """
+
+    __slots__ = (
+        'mol_name',
+        'atoms',
+        'position',
+        'energy',
+        'gradient',
+        'hessian',
+        'total_mass',
+        'com'
+    )
+
+    def __init__(self, mol_name: str) -> None:
+        """
+        The constructor function __init__ is used to create new instances of the Molecule class.
+
+        Args:
+            mol_name (str): Molecule name
+        """
+        self.mol_name: str = mol_name
+
+        self.atoms: List[Atom] = []
+
+    def __iter__(self):
+        for atom in self.atoms:
+            yield atom
+
+    def add_atom(self, a: Atom) -> None:
+        """
+        Adds Atom instance to the molecule.
+
+        Args:
+            a (Atom): Atom instance to be added to atom list.
+        """
+
+        self.atoms.append(a)
+
+    def update_charges(self, charges: List[float]) -> None:
+
+        if len(charges) != len(self.atoms):
+            raise ValueError(
+                f"The number of charges ({len(charges)}) does not match the number of atoms ({len(self.atoms)})"
+            )
+
+        for i, atom in enumerate(self.atoms):
+            atom.chg = charges[i]
+
+    def print_mol_info(self, fh: TextIO) -> None:
+        """
+        Prints the Molecule information into a Output File
+
+        Args:
+            fh (TextIO): Output File
+        """
+
+        fh.write("-"*80)
+        fh.write(f"Molecule Name: {self.mol_name}\n")
+
+        fh.write("\n")
+
+        for atom in self.atoms:
+            fh.write(f"{atom_symbol[atom.na]}   r: [{atom.rx}, {atom.ry}, {atom.rz}]    charge: {atom.chg}")

crystalpol/shared/utils/ptable.py

@@ -0,0 +1,34 @@
+# Tuple of atom symbols
+atom_symbol = ("00",
+
+               "H ", "He",
+               "Li", "Be", "B ", "C ", "N ", "O ", "F ", "Ne",
+               "Na", "Mg", "Al", "Si", "P ", "S ", "Cl", "Ar",
+               "K ", "Ca", "Sc", "Ti", "V ", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br",
+               "Kr",
+               "Rb", "Sr", "Y ", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I ",
+               "Xe",
+               "Cs", "Ba",
+               "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu",
+               "Hf", "Ta", "W ", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Ti", "Pb", "Bi", "Po", "At", "Rn",
+               "Fr", "Ra",
+               "Ac", "Th", "Pa", "U ", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr")
+
+# Tuple of atom masses
+atom_mass = (0.0,
+
+             1.0079, 4.0026,
+             6.9410, 9.0122, 10.811, 12.011, 14.007, 15.999, 18.998, 20.180,
+             22.990, 24.305, 26.982, 28.086, 30.974, 32.065, 35.453, 39.948,
+             39.098, 40.078, 44.956, 47.867, 50.942, 51.996, 54.938, 55.845, 58.933, 58.693, 63.546, 65.409, 69.723,
+             72.640, 74.922, 78.960, 79.904, 83.798,
+             85.468, 87.620, 88.906, 91.224, 92.906, 95.940, 98.000, 101.07, 102.91, 106.42, 107.87, 112.41, 114.82,
+             118.71, 121.76, 127.60, 126.90, 131.29,
+             132.91, 137.33,
+             138.91, 140.12, 140.91, 144.24, 145.00, 150.36, 151.96, 157.25, 158.93, 162.50, 164.93, 167.26, 168.93,
+             173.04, 174.97,
+             178.49, 180.95, 183.84, 186.21, 190.23, 192.22, 195.08, 196.97, 200.59, 204.38, 207.20, 208.98, 209.00,
+             210.00, 222.00,
+             223.00, 226.00,
+             227.00, 232.04, 231.04, 238.03, 237.00, 244.00, 243.00, 247.00, 247.00, 251.00, 252.00, 257.00, 258.00,
+             259.00, 262.00)