Read Potentials and Initial Testing

diceplayer/shared/config/dice_dto.py

@@ -0,0 +1,54 @@
+from diceplayer.shared.utils.dataclass_protocol import Dataclass
+
+from dataclasses import dataclass
+from dacite import from_dict
+from typing import List
+
+
+@dataclass
+class DiceDTO(Dataclass):
+
+    ljname: str
+    outname: str
+    ncores: int
+    dens: float
+    nmol: List[int]
+    nstep: List[int]
+
+    upbuf = 360
+    combrule = "*"
+    isave: int = 1000
+    press: float = 1.0
+    temp: float = 300.0
+    randominit: str = 'first'
+
+    def __post_init__(self):
+
+        if self.ljname is None:
+            raise ValueError(
+                "Error: 'ljname' keyword not specified in config file"
+            )
+
+        if self.outname is None:
+            raise ValueError(
+                "Error: 'outname' keyword not specified in config file"
+            )
+
+        if self.dens is None:
+            raise ValueError(
+                "Error: 'dens' keyword not specified in config file"
+            )
+
+        if self.nmol == 0:
+            raise ValueError(
+                "Error: 'nmol' keyword not defined appropriately in config file"
+            )
+
+        if self.nstep == 0:
+            raise ValueError(
+                "Error: 'nstep' keyword not defined appropriately in config file"
+            )
+
+    @classmethod
+    def from_dict(cls, param: dict):
+        return from_dict(DiceDTO, param)

diceplayer/shared/config/gaussian_dto.py

@@ -0,0 +1,28 @@
+from diceplayer.shared.utils.dataclass_protocol import Dataclass
+
+from dataclasses import dataclass
+from dacite import from_dict
+
+
+@dataclass
+class GaussianDTO(Dataclass):
+    level: str
+    qmprog: str
+    keywords: str
+
+    chgmult = [0, 1]
+    pop: str = 'chelpg'
+
+    def __post_init__(self):
+        if self.qmprog not in ("g03", "g09", "g16"):
+            raise ValueError(
+                "Error: invalid qmprog value."
+            )
+        if self.level is None:
+            raise ValueError(
+                "Error: 'level' keyword not specified in config file."
+            )
+
+    @classmethod
+    def from_dict(cls, param: dict):
+        return from_dict(GaussianDTO, param)

diceplayer/shared/config/player_dto.py

@@ -0,0 +1,24 @@
+from diceplayer.shared.utils.dataclass_protocol import Dataclass
+
+from dataclasses import dataclass
+from dacite import from_dict
+
+
+@dataclass
+class PlayerDTO(Dataclass):
+    opt: bool
+    maxcyc: int
+    nprocs: int
+
+    qmprog: str = 'g16'
+    altsteps: int = 20000
+    simulation_dir = 'simfiles'
+
+    def __post_init__(self):
+        MIN_STEP = 20000
+        # altsteps value is always the nearest multiple of 1000
+        self.altsteps = round(max(MIN_STEP, self.altsteps) / 1000) * 1000
+
+    @classmethod
+    def from_dict(cls, param: dict):
+        return from_dict(PlayerDTO, param)

diceplayer/shared/environment/atom.py

@@ -0,0 +1,52 @@
+from diceplayer.shared.utils.ptable import atommass
+
+
+class Atom:
+    """
+    Atom class declaration. This class is used throughout the DicePlayer program to represent atoms.
+
+    Atributes:
+        lbl (int): Dice derived variable used to represent atoms with identical energies and simetric positions.
+        na (int): Atomic number 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.
+        chg (float): charge of the represented atom.
+        eps (float): quantum number epsilon of the represented atom.
+        sig (float): quantum number sigma of the represented atom.
+    """
+
+    def __init__(
+            self,
+            lbl: int,
+            na: int,
+            rx: float,
+            ry: float,
+            rz: float,
+            chg: float,
+            eps: float,
+            sig: float,
+    ) -> None:
+        """
+        The constructor function __init__ is used to create new instances of the Atom class.
+
+        Args:
+            lbl (int): Dice derived variable used to represent atoms with identical energies and simetric positions.
+            na (int): Atomic number 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.
+            chg (float): charge of the represented atom.
+            eps (float): quantum number epsilon of the represented atom.
+            sig (float): quantum number sigma of the represented atom.
+        """
+
+        self.lbl = lbl
+        self.na = na
+        self.rx = rx
+        self.ry = ry
+        self.rz = rz
+        self.chg = chg
+        self.eps = eps
+        self.sig = sig
+        self.mass = atommass[self.na]

diceplayer/shared/environment/molecule.py

@@ -0,0 +1,406 @@
+import logging
+import math
+from copy import deepcopy
+from typing import List, Any, Tuple, Final, Union
+
+import numpy as np
+from nptyping import NDArray, Shape, Float
+from numpy.linalg import linalg
+
+from diceplayer.shared.environment.atom import Atom
+from diceplayer.shared.utils.misc import BOHR2ANG
+from diceplayer.shared.utils.ptable import ghost_number
+
+
+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
+    """
+
+    def __init__(self, molname: str) -> None:
+        """
+        The constructor function __init__ is used to create new instances of the Molecule class.
+
+        Args:
+            molname (str): Molecule name
+        """
+        self.molname: str = molname
+
+        self.atom: List[Atom] = []
+        self.position: NDArray[Any, Any]
+        self.energy: NDArray[Any, Any]
+        self.gradient: NDArray[Any, Any]
+        self.hessian: NDArray[Any, Any]
+
+        self.ghost_atoms: List[Atom] = []
+        self.lp_atoms: List[Atom] = []
+
+        self.total_mass: int = 0
+        self.com: Union[None, NDArray[Any, Any]] = None
+
+    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.atom.append(a)
+        self.total_mass += a.mass
+
+        self.center_of_mass()
+
+    def center_of_mass(self) -> NDArray[Any, Any]:
+        """
+        Calculates the center of mass of the molecule
+        """
+
+        self.com = np.zeros(3)
+
+        for atom in self.atom:
+            self.com += atom.mass * np.array([atom.rx, atom.ry, atom.rz])
+
+        self.com = self.com / self.total_mass
+
+        return self.com
+
+    def center_of_mass_to_origin(self) -> None:
+        """
+        Updated positions based on the center of mass of the molecule
+        """
+
+        self.center_of_mass()
+
+        for atom in self.atom:
+            atom.rx -= self.com[0]
+            atom.ry -= self.com[1]
+            atom.rz -= self.com[2]
+
+    def charges_and_dipole(self) -> List[float]:
+        """
+        Calculates the charges and dipole of the molecule atoms
+
+        Returns:
+            List[float]: Respectivly magnitude of the: charge magnitude, first dipole,
+            second dipole, third dipole and total dipole.
+        """
+
+        eA_to_Debye = 1 / 0.20819434
+        charge = 0
+        dipole = np.zeros(3)
+        for atom in self.atom:
+            position = np.array([atom.rx, atom.ry, atom.rz])
+            dipole += atom.chg * position
+            charge += atom.chg
+
+        dipole *= eA_to_Debye
+        total_dipole = math.sqrt(dipole[0] ** 2 + dipole[1] ** 2 + dipole[2] ** 2)
+
+        return [charge, dipole[0], dipole[1], dipole[2], total_dipole]
+
+    def distances_between_atoms(self) -> NDArray[Shape["Any,Any"], Float]:
+        """
+        Calculates distances between the atoms of the molecule
+
+        Returns:
+           NDArray[Shape["Any,Any"],Float]: distances between the atoms.
+        """
+
+        distances = []
+        dim = len(self.atom)
+        for atom1 in self.atom:
+            if atom1.na != ghost_number:
+                for atom2 in self.atom:
+                    if atom2.na != ghost_number:
+                        dx = atom1.rx - atom2.rx
+                        dy = atom1.ry - atom2.ry
+                        dz = atom1.rz - atom2.rz
+                        distances.append(math.sqrt(dx ** 2 + dy ** 2 + dz ** 2))
+
+        return np.array(distances).reshape(dim, dim)
+
+    def inertia_tensor(self) -> NDArray[Shape["3, 3"], Float]:
+        """
+        Calculates the inertia tensor of the molecule.
+
+        Returns:
+            NDArray[Shape["3, 3"], Float]: inertia tensor of the molecule.
+        """
+
+        self.center_of_mass()
+        Ixx = Ixy = Ixz = Iyy = Iyz = Izz = 0.0
+
+        for atom in self.atom:
+            dx = atom.rx - self.com[0]
+            dy = atom.ry - self.com[1]
+            dz = atom.rz - self.com[2]
+
+            Ixx += atom.mass * (dy ** 2 + dz ** 2)
+            Iyy += atom.mass * (dz ** 2 + dx ** 2)
+            Izz += atom.mass * (dx ** 2 + dy ** 2)
+
+            Ixy += atom.mass * dx * dy * -1
+            Ixz += atom.mass * dx * dz * -1
+            Iyz += atom.mass * dy * dz * -1
+
+        return np.array([[Ixx, Ixy, Ixz], [Ixy, Iyy, Iyz], [Ixz, Iyz, Izz]])
+
+    def axes(self) -> NDArray[Shape["3, 3"], Float]:
+        """
+        Calculates the axes of the molecule
+
+        Returns:
+           NDArray[Shape["3, 3"], Float]: Returns the axes of molecule
+        """
+
+        eixos = np.zeros(3)
+        if len(self.atom) == 2:
+
+            position1 = np.array([self.atom[0].rx, self.atom[0].ry, self.atom[0].rz])
+            position2 = np.array([self.atom[1].rx, self.atom[1].ry, self.atom[1].rz])
+            eixos = position2 - position1
+            eixos /= linalg.norm(eixos)
+
+        elif len(self.atom) > 2:
+
+            position1 = np.array([self.atom[0].rx, self.atom[0].ry, self.atom[0].rz])
+            position2 = np.array([self.atom[1].rx, self.atom[1].ry, self.atom[1].rz])
+            position3 = np.array([self.atom[2].rx, self.atom[2].ry, self.atom[2].rz])
+            v1 = position2 - position1
+            v2 = position3 - position1
+            v3 = np.cross(v1, v2)
+            v2 = np.cross(v1, v3)
+            v1 /= linalg.norm(v1)
+            v2 /= linalg.norm(v2)
+            v3 /= linalg.norm(v3)
+            eixos = np.array(
+                [[v1[0], v1[1], v1[2]], [v2[0], v2[1], v2[2]], [v3[0], v3[1], v3[2]]]
+            )
+
+        return eixos
+
+    def principal_axes(self) -> Tuple[np.ndarray, np.ndarray]:
+        """
+        Calculates the principal axes of the molecule
+
+        Returns:
+            Tuple[np.ndarray, np.ndarray]: Tuple where the first value is the Eigen Values and the second is the Eigen Vectors,
+             representing the principal axes of the molecule.
+        """
+
+        try:
+            evals, evecs = linalg.eigh(self.inertia_tensor())
+        except:
+            raise RuntimeError("Error: diagonalization of inertia tensor did not converge")
+
+        return evals, evecs
+
+    def read_position(self) -> np.ndarray:
+        """Reads the position of the molecule from the position values of the atoms
+
+        Returns:
+            np.ndarray: internal position relative to atoms of the molecule
+        """
+
+        position_list = []
+        for atom in self.atom:
+            position_list.extend([atom.rx, atom.ry, atom.rz])
+        position = np.array(position_list)
+        position *= BOHR2ANG
+
+        return position
+
+    def updateCharges(self, charges: List[float]) -> None:
+
+        for i, atom in enumerate(self.atom):
+            atom.chg = charges[i]
+
+    def update_hessian(
+            self,
+            step: np.ndarray,
+            cur_gradient: np.ndarray,
+            old_gradient: np.ndarray,
+            hessian: np.ndarray,
+    ) -> np.ndarray:
+        """
+        Updates the Hessian of the molecule based on the current hessian, the current gradient and the previous gradient
+
+        Args:
+            step (np.ndarray): step value of the iteration
+            cur_gradient (np.ndarray): current gradient
+            old_gradient (np.ndarray): previous gradient
+            hessian (np.ndarray): current hessian
+
+        Returns:
+            np.ndarray: updated hessian of the molecule
+        """
+
+        dif_gradient = cur_gradient - old_gradient
+
+        mat1 = 1 / np.dot(dif_gradient, step) * np.matmul(dif_gradient.T, dif_gradient)
+        mat2 = 1 / np.dot(step, np.matmul(hessian, step.T).T)
+        mat2 *= np.matmul(np.matmul(hessian, step.T), np.matmul(step, hessian))
+
+        return hessian + mat1 - mat2
+
+    def sizes_of_molecule(self) -> List[float]:
+        """
+        Calculates sides of the smallest box that the molecule could fit
+
+        Returns:
+            List[float]: list of the sizes of the molecule
+        """
+
+        x_list = []
+        y_list = []
+        z_list = []
+
+        for atom in self.atom:
+            if atom.na != ghost_number:
+                x_list.append(atom.rx)
+                y_list.append(atom.ry)
+                z_list.append(atom.rz)
+
+        x_max = max(x_list)
+        x_min = min(x_list)
+        y_max = max(y_list)
+        y_min = min(y_list)
+        z_max = max(z_list)
+        z_min = min(z_list)
+
+        sizes = [x_max - x_min, y_max - y_min, z_max - z_min]
+
+        return sizes
+
+    def standard_orientation(self) -> None:
+        """
+        Rotates the molecule to the standard orientation
+        """
+
+        self.center_of_mass_to_origin()
+        evals, evecs = self.principal_axes()
+
+        if round(linalg.det(evecs)) == -1:
+            evecs[0, 2] *= -1
+            evecs[1, 2] *= -1
+            evecs[2, 2] *= -1
+
+        if round(linalg.det(evecs)) != 1:
+            raise RuntimeError(
+                "Error: could not make a rotation matrix while adopting the standard orientation"
+            )
+
+        rot_matrix = evecs.T
+
+        for atom in self.atom:
+            position = np.array([atom.rx, atom.ry, atom.rz])
+            new_position = np.matmul(rot_matrix, position.T).T
+
+            atom.rx = new_position[0]
+            atom.ry = new_position[1]
+            atom.rz = new_position[2]
+
+    def translate(self, vector: np.ndarray) -> "Molecule":
+        """
+        Creates a new Molecule object where its' atoms has been translated by a vector
+
+        Args:
+            vector (np.ndarray): translation vector
+
+        Returns:
+            Molecule: new Molecule object translated by a vector
+        """
+
+        new_molecule = deepcopy(self)
+
+        for atom in new_molecule.atom:
+            atom.rx += vector[0]
+            atom.ry += vector[1]
+            atom.rz += vector[2]
+
+        return new_molecule
+
+    def print_mol_info(self) -> None:
+        """
+        Prints the Molecule information into a Output File
+        """
+
+        logging.info(
+            "    Center of mass = ( {:>10.4f} , {:>10.4f} , {:>10.4f} )\n".format(
+                self.com[0], self.com[1], self.com[2]
+            )
+        )
+        inertia = self.inertia_tensor()
+        evals, evecs = self.principal_axes()
+
+        logging.info(
+            "    Moments of inertia =  {:>9E}  {:>9E}  {:>9E}\n".format(
+                evals[0], evals[1], evals[2]
+            )
+        )
+
+        logging.info(
+            "    Major principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
+                evecs[0, 0], evecs[1, 0], evecs[2, 0]
+            )
+        )
+        logging.info(
+            "    Inter principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
+                evecs[0, 1], evecs[1, 1], evecs[2, 1]
+            )
+        )
+        logging.info(
+            "    Minor principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
+                evecs[0, 2], evecs[1, 2], evecs[2, 2]
+            )
+        )
+
+        sizes = self.sizes_of_molecule()
+        logging.info(
+            "    Characteristic lengths = ( {:>6.2f} , {:>6.2f} , {:>6.2f} )\n".format(
+                sizes[0], sizes[1], sizes[2]
+            )
+        )
+        logging.info("    Total mass = {:>8.2f} au\n".format(self.total_mass))
+
+        chg_dip = self.charges_and_dipole()
+        logging.info("    Total charge = {:>8.4f} e\n".format(chg_dip[0]))
+        logging.info(
+            "    Dipole moment = ( {:>9.4f} , {:>9.4f} , {:>9.4f} )     Total = {:>9.4f} Debye\n\n".format(
+                chg_dip[1], chg_dip[2], chg_dip[3], chg_dip[4]
+            )
+        )
+
+    def minimum_distance(self, molec: 'Molecule') -> float:
+        """
+        Return the minimum distance between two molecules
+
+        Args:
+            molec (Molecule): Molecule object to be compared
+
+        Returns:
+            float: minimum distance between the two molecules
+        """
+
+        distances = []
+        for atom1 in self.atom:
+            if atom1.na != ghost_number:
+                for atom2 in molec.atom:
+                    if atom2.na != ghost_number:
+                        dx = atom1.rx - atom2.rx
+                        dy = atom1.ry - atom2.ry
+                        dz = atom1.rz - atom2.rz
+                        distances.append(math.sqrt(dx ** 2 + dy ** 2 + dz ** 2))
+
+        return min(distances)

diceplayer/shared/environment/system.py

@@ -0,0 +1,236 @@
+import math
+from copy import deepcopy
+from typing import List, Tuple, TextIO
+
+import numpy as np
+from numpy import linalg
+
+from diceplayer.shared.environment.molecule import Molecule
+from diceplayer.shared.utils.misc import BOHR2ANG
+from diceplayer.shared.utils.ptable import atomsymb
+
+
+class System:
+    """
+    System class declaration. This class is used throughout the DicePlayer program to represent the system containing the molecules.
+
+    Atributes:
+        molecule (List[Molecule]): List of molecules of the system
+        nmols (List[int]): List of number of molecules in the system
+    """
+
+    def __init__(self) -> None:
+        """
+        Initializes a empty system object that will be populated afterwards
+        """
+
+        self.molecule: List[Molecule] = []
+        self.nmols: List[int] = []
+
+    def add_type(self, nmols: int, m: Molecule) -> None:
+        """
+        Adds a new molecule type to the system
+
+        Args:
+            nmols (int): Number of molecules of the new type in the system
+            m (Molecule): The instance of the new type of molecule
+        """
+        self.molecule.append(m)
+        self.nmols.append(nmols)
+
+    def center_of_mass_distance(self, a: int, b: int) -> float:
+        """
+        Calculates the distance between the center of mass of two molecules
+
+        Args:
+            a (Molecule): First Molecule Instance
+            b (Molecule): Second Molecule Instance
+
+        Returns:
+            float: module of the distance between the two center of masses
+        """
+
+        com1 = self.molecule[a].center_of_mass()
+        com2 = self.molecule[b].center_of_mass()
+        dx = com1[0] - com2[0]
+        dy = com1[1] - com2[1]
+        dz = com1[2] - com2[2]
+        distance = math.sqrt(dx**2 + dy**2 + dz**2)
+
+        return distance
+
+    def rmsd_fit(self, p_index: int, r_index: int) -> Tuple[float, Molecule]:
+
+        projecting_mol = self.molecule[p_index]
+        reference_mol = self.molecule[r_index]
+
+        if len(projecting_mol.atom) != len(reference_mol.atom):
+            raise RuntimeError(
+                "Error in RMSD fit procedure: molecules have different number of atoms"
+            )
+        dim = len(projecting_mol.atom)
+
+        new_projecting_mol = deepcopy(projecting_mol)
+        new_reference_mol = deepcopy(reference_mol)
+
+        new_projecting_mol.center_of_mass_to_origin()
+        new_reference_mol.center_of_mass_to_origin()
+
+        x = []
+        y = []
+
+        for atom in new_projecting_mol.atom:
+            x.extend([atom.rx, atom.ry, atom.rz])
+
+        for atom in new_reference_mol.atom:
+            y.extend([atom.rx, atom.ry, atom.rz])
+
+        x = np.array(x).reshape(dim, 3)
+        y = np.array(y).reshape(dim, 3)
+
+        r = np.matmul(y.T, x)
+        rr = np.matmul(r.T, r)
+
+        try:
+            evals, evecs = linalg.eigh(rr)
+        except:
+            raise RuntimeError("Error: diagonalization of RR matrix did not converge")
+
+        a1 = evecs[:, 2].T
+        a2 = evecs[:, 1].T
+        a3 = np.cross(a1, a2)
+
+        A = np.array([a1[0], a1[1], a1[2], a2[0], a2[1], a2[2], a3[0], a3[1], a3[2]])
+        A = A.reshape(3, 3)
+
+        b1 = np.matmul(r, a1.T).T  # or np.dot(r, a1)
+        b1 /= linalg.norm(b1)
+        b2 = np.matmul(r, a2.T).T  # or np.dot(r, a2)
+        b2 /= linalg.norm(b2)
+        b3 = np.cross(b1, b2)
+
+        B = np.array([b1[0], b1[1], b1[2], b2[0], b2[1], b2[2], b3[0], b3[1], b3[2]])
+        B = B.reshape(3, 3).T
+
+        rot_matrix = np.matmul(B, A)
+        x = np.matmul(rot_matrix, x.T).T
+
+        rmsd = 0
+        for i in range(dim):
+            rmsd += (
+                (x[i, 0] - y[i, 0]) ** 2
+                + (x[i, 1] - y[i, 1]) ** 2
+                + (x[i, 2] - y[i, 2]) ** 2
+            )
+        rmsd = math.sqrt(rmsd / dim)
+
+        for i in range(dim):
+            new_projecting_mol.atom[i].rx = x[i, 0]
+            new_projecting_mol.atom[i].ry = x[i, 1]
+            new_projecting_mol.atom[i].rz = x[i, 2]
+
+        reference_mol.center_of_mass()
+
+        projected_mol = new_projecting_mol.translate(reference_mol.com)
+
+        return rmsd, projected_mol
+
+    def update_molecule(self, position: np.ndarray, fh: TextIO) -> None:
+        """Updates the position of the molecule in the Output file
+
+        Args:
+            position (np.ndarray): numpy position vector
+            fh (TextIO): Output file
+        """
+
+        position_in_ang = (position * BOHR2ANG).tolist()
+        self.add_type(self.nmols[0], deepcopy(self.molecule[0]))
+
+        for atom in self.molecule[-1].atom:
+
+            atom.rx = position_in_ang.pop(0)
+            atom.ry = position_in_ang.pop(0)
+            atom.rz = position_in_ang.pop(0)
+
+        rmsd, self.molecule[0] = self.rmsd_fit(-1, 0)
+        self.molecule.pop(-1)
+
+        fh.write("\nProjected new conformation of reference molecule with RMSD fit\n")
+        fh.write("RMSD = {:>8.5f} Angstrom\n".format(rmsd))
+
+    def nearest_image(
+        self,
+        index_r: int,
+        index_m: int,
+        lx: float,
+        ly: float,
+        lz: float,
+        criterium=None,
+    ) -> Tuple[float, Molecule]:
+
+        if criterium in None:
+            criterium = "com"
+
+        if criterium != "com" and criterium != "min":
+            raise RuntimeError("Error in value passed to function nearest_image")
+
+        min_dist = 1e20
+
+        for i in range(-1, 2):
+            for j in range(-1, 2):
+                for k in range(-1, 2):
+
+                    tr_vector = [i * lx, j * ly, k * lz]
+                    self.add_molecule(self.molecule[index_m].translate(tr_vector))
+
+                    if criterium == "com":
+                        dist = self.center_of_mass_distance(index_r, -1)
+                    else:
+                        dist = self.minimum_distance(index_r, -1)
+
+                    if dist < min_dist:
+                        min_dist = dist
+                        nearestmol = deepcopy(self.molecule[-1])
+
+        self.molecule.pop(-1)
+
+        return min_dist, nearestmol
+
+    def print_geom(self, cycle: int, fh: TextIO) -> None:
+        """
+        Print the geometry of the molecule in the Output file
+
+        Args:
+            cycle (int): Number of the cycle
+            fh (TextIO): Output file
+        """
+
+        fh.write("Cycle # {}\n".format(cycle))
+        fh.write("Number of site: {}\n".format(len(self.molecule[0].atom)))
+        for atom in self.molecule[0].atom:
+            symbol = atomsymb[atom.na]
+            fh.write(
+                "{:<2s}    {:>10.6f}  {:>10.6f}  {:>10.6f}\n".format(
+                    symbol, atom.rx, atom.ry, atom.rz
+                )
+            )
+
+    def printChargesAndDipole(self, cycle: int, fh: TextIO) -> None:
+        """
+        Print the charges and dipole of the molecule in the Output file
+
+        Args:
+            cycle (int): Number of the cycle
+            fh (TextIO): Output file
+        """
+
+        fh.write("Cycle # {}\n".format(cycle))
+        fh.write("Number of site: {}\n".format(len(self.molecule[0].atom)))
+
+        chargesAndDipole = self.molecule[0].charges_and_dipole()
+
+        fh.write(
+            "{:>10.6f}  {:>10.6f}  {:>10.6f}  {:>10.6f}  {:>10.6f}\n".format(
+                chargesAndDipole[0], chargesAndDipole[1], chargesAndDipole[2], chargesAndDipole[3], chargesAndDipole[4]
+            )
+        )

diceplayer/shared/external/__external.py

@@ -0,0 +1,26 @@
+from diceplayer.shared.utils.dataclass_protocol import Dataclass
+
+from abc import ABC, abstractmethod
+
+
+class External(ABC):
+    __slots__ = [
+        'config'
+    ]
+
+    @abstractmethod
+    def __init__(self, data: dict):
+        pass
+
+    @staticmethod
+    @abstractmethod
+    def set_config(data: dict) -> Dataclass:
+        pass
+
+    @abstractmethod
+    def start(self):
+        pass
+
+    @abstractmethod
+    def reset(self):
+        pass

diceplayer/shared/external/dice.py

@@ -0,0 +1,22 @@
+from diceplayer.shared.utils.dataclass_protocol import Dataclass
+from diceplayer.shared.external.__external import External
+from diceplayer.shared.config.dice_dto import DiceDTO
+
+
+class Dice(External):
+
+    def __init__(self, data: dict):
+        self.config: DiceDTO = self.set_config(data)
+
+    @staticmethod
+    def set_config(data: dict) -> DiceDTO:
+        return DiceDTO.from_dict(data)
+
+    def configure(self):
+        pass
+
+    def start(self):
+        pass
+
+    def reset(self):
+        pass

diceplayer/shared/external/gaussian.py

@@ -0,0 +1,23 @@
+from diceplayer.shared.utils.dataclass_protocol import Dataclass
+from diceplayer.shared.config.gaussian_dto import GaussianDTO
+from diceplayer.shared.external.__external import External
+
+
+class Gaussian(External):
+
+    def __init__(self, data: dict):
+        self.config: GaussianDTO = self.set_config(data)
+
+    @staticmethod
+    def set_config(data: dict) -> GaussianDTO:
+        return GaussianDTO.from_dict(data)
+
+    def configure(self):
+        pass
+
+    def start(self):
+        pass
+
+    def reset(self):
+        pass
+

diceplayer/shared/utils/dataclass_protocol.py

@@ -0,0 +1,6 @@
+from typing import runtime_checkable, Protocol
+
+
+@runtime_checkable
+class Dataclass(Protocol):
+    __dataclass_fields__: dict

diceplayer/shared/utils/misc.py

@@ -0,0 +1,64 @@
+import gzip
+import os
+import shutil
+import sys
+import time
+from typing import Final
+
+#######################################  constants  ######################################
+
+
+BOHR2ANG: Final[float] = 0.52917721092
+ANG2BOHR: Final[float] = 1 / BOHR2ANG
+
+
+#######################################  functions  ######################################
+
+def weekday_date_time():
+    return time.strftime("%A, %d %b %Y at %H:%M:%S")
+
+
+def date_time():
+    return time.strftime("%d %b %Y at %H:%M:%S")
+
+
+def compress_files_1mb(path):
+    working_dir = os.getcwd()
+    os.chdir(path)
+
+    files = filter(os.path.isfile, os.listdir(os.curdir))
+    for file in files:
+        if os.path.getsize(file) > 1024 * 1024:  ## If bigger than 1MB
+            filegz = file + ".gz"
+            try:
+                with open(file, 'rb') as f_in:
+                    with gzip.open(filegz, 'wb') as f_out:
+                        shutil.copyfileobj(f_in, f_out)
+            except:
+                sys.exit("Error: cannot compress file {}".format(file))
+
+    os.chdir(working_dir)
+
+    return
+
+
+def make_step_dir(cycle):
+    sim_dir = "simfiles"
+    step_dir = "step{:02d}".format(cycle)
+    path = sim_dir + os.sep + step_dir
+    if os.path.exists(path):
+        sys.exit("Error: a file or directory {} already exists".format(step_dir))
+    try:
+        os.makedirs(path)
+    except:
+        sys.exit("Error: cannot make directory {}".format(step_dir))
+
+
+def make_qm_dir(cycle):
+    sim_dir = "simfiles"
+    step_dir = "step{:02d}".format(cycle)
+    path = sim_dir + os.sep + step_dir + os.sep + "qm"
+    try:
+        os.makedirs(path)
+    except:
+        sys.exit("Error: cannot make directory {}".format(path))

diceplayer/shared/utils/ptable.py

@@ -0,0 +1,35 @@
+####  Label used in Dice for a ghost atom
+dice_ghost_label = "Xx"
+
+####  Tuple of atom symbols
+atomsymb = ( "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",
+               dice_ghost_label )
+               
+####  Tuple of atom masses  
+atommass = ( 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,
+                     0.000 )
+
+####  Number of the ghost atom
+ghost_number = len(atomsymb) - 1

diceplayer/shared/utils/validations.py

@@ -0,0 +1,15 @@
+def NotNull(requiredArgs=[]):
+    def _NotNull(function):
+        def wrapper(*args, **kwargs):
+            for arg in requiredArgs:
+                try:
+                    assert (
+                        kwargs.get(arg) is not None
+                    ), "Invalid Config File. Keyword {} is required".format(arg)
+                except AssertionError as err:
+                    print(err)
+            return function(*args, **kwargs)
+
+        return wrapper
+
+    return _NotNull