chore: better project structure

diceplayer/environment/__init__.py

@@ -0,0 +1,6 @@
+from .atom import Atom
+from .molecule import Molecule
+from .system import System
+
+
+__all__ = ["Atom", "Molecule", "System"]

diceplayer/environment/atom.py

@@ -0,0 +1,52 @@
+from diceplayer.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/environment/molecule.py

@@ -0,0 +1,385 @@
+from __future__ import annotations
+
+from typing_extensions import TYPE_CHECKING
+
+
+if TYPE_CHECKING:
+    from nptyping import Float, NDArray, Shape
+
+from diceplayer import logger
+from diceplayer.environment.atom import Atom
+from diceplayer.utils.misc import BOHR2ANG
+from diceplayer.utils.ptable import ghost_number
+
+import numpy as np
+from numpy.linalg import linalg
+from typing_extensions import Any, List, Tuple, Union
+
+import math
+from copy import deepcopy
+
+
+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
+        """
+        for atom in self.atom:
+            atom.rx -= self.com[0]
+            atom.ry -= self.com[1]
+            atom.rz -= self.com[2]
+
+        self.center_of_mass()
+
+    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 index1, atom1 in enumerate(self.atom):
+            for index2, atom2 in enumerate(self.atom):
+                if index1 != index2:
+                    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 - 1)
+
+    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 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 ValueError:
+            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 update_charges(self, charges: NDArray) -> int:
+        """
+        Updates the charges of the atoms of the molecule and
+        returns the max difference between the new and old charges
+        """
+        diff = 0
+        for i, atom in enumerate(self.atom):
+            diff = max(diff, abs(atom.chg - charges[i]))
+            atom.chg = charges[i]
+
+        return diff
+
+    # @staticmethod
+    # def update_hessian(
+    #         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:
+            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
+        """
+
+        logger.info(
+            "    Center of mass = ( {:>10.4f} , {:>10.4f} , {:>10.4f} )".format(
+                self.com[0], self.com[1], self.com[2]
+            )
+        )
+        self.inertia_tensor()
+        evals, evecs = self.principal_axes()
+
+        logger.info(
+            "    Moments of inertia =  {:>9E}  {:>9E}  {:>9E}".format(
+                evals[0], evals[1], evals[2]
+            )
+        )
+
+        logger.info(
+            "    Major principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )".format(
+                evecs[0, 0], evecs[1, 0], evecs[2, 0]
+            )
+        )
+        logger.info(
+            "    Inter principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )".format(
+                evecs[0, 1], evecs[1, 1], evecs[2, 1]
+            )
+        )
+        logger.info(
+            "    Minor principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )".format(
+                evecs[0, 2], evecs[1, 2], evecs[2, 2]
+            )
+        )
+
+        sizes = self.sizes_of_molecule()
+        logger.info(
+            "    Characteristic lengths = ( {:>6.2f} , {:>6.2f} , {:>6.2f} )".format(
+                sizes[0], sizes[1], sizes[2]
+            )
+        )
+        logger.info("    Total mass = {:>8.2f} au".format(self.total_mass))
+
+        chg_dip = self.charges_and_dipole()
+        logger.info("    Total charge = {:>8.4f} e".format(chg_dip[0]))
+        logger.info(
+            "    Dipole moment = ( {:>9.4f} , {:>9.4f} , {:>9.4f} )     Total = {:>9.4f} Debye".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/environment/system.py

@@ -0,0 +1,238 @@
+from diceplayer import logger
+from diceplayer.environment.molecule import Molecule
+from diceplayer.utils.misc import BOHR2ANG
+
+import numpy as np
+from numpy import linalg
+from typing_extensions import List, Tuple
+
+import math
+from copy import deepcopy
+
+
+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 an empty system object that will be populated afterwards
+        """
+        self.nmols: List[int] = []
+        self.molecule: List[Molecule] = []
+
+    def add_type(self, m: Molecule) -> None:
+        """
+        Adds a new molecule type to the system
+
+        Args:
+            m (Molecule): The instance of the new type of molecule
+        """
+        if isinstance(m, Molecule) is False:
+            raise TypeError("Error: molecule is not a Molecule instance")
+        self.molecule.append(m)
+
+    def update_molecule(self, position: np.ndarray) -> None:
+        """Updates the position of the molecule in the Output file
+
+        Args:
+            position (np.ndarray): numpy position vector
+        """
+
+        position_in_ang = (position * BOHR2ANG).tolist()
+        self.add_type(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)
+
+        logger.info("Projected new conformation of reference molecule with RMSD fit")
+        logger.info(f"RMSD = {rmsd:>8.5f} Angstrom")
+
+    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 Exception as err:
+            raise RuntimeError(
+                "Error: diagonalization of RR matrix did not converge"
+            ) from err
+
+        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 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 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 print_charges_and_dipole(self, cycle: int) -> None:
+        """
+        Print the charges and dipole of the molecule in the Output file
+
+        Args:
+            cycle (int): Number of the cycle
+            fh (TextIO): Output file
+        """
+
+        logger.info("Cycle # {}\n".format(cycle))
+        logger.info("Number of site: {}\n".format(len(self.molecule[0].atom)))
+
+        chargesAndDipole = self.molecule[0].charges_and_dipole()
+
+        logger.info(
+            "{:>10.6f}  {:>10.6f}  {:>10.6f}  {:>10.6f}  {:>10.6f}\n".format(
+                chargesAndDipole[0],
+                chargesAndDipole[1],
+                chargesAndDipole[2],
+                chargesAndDipole[3],
+                chargesAndDipole[4],
+            )
+        )