Merge branch 'devel-folder-rework' into devel

control.in

@@ -1,22 +0,0 @@
-# diceplayer
-maxcyc = 3
-opt = no
-nprocs = 4
-qmprog = g16
-lps = no
-ghosts = no
-altsteps = 20000
-
-# dice
-ncores = 3
-nmol = 1 50
-dens = 0.75
-nstep = 4000 6000
-isave = 1000
-ljname = phb.ljc
-outname = phb
-randominit = first
-
-# Gaussian
-level = MP2/aug-cc-pVDZ
-keywords = freq=intmodes

control.yml

@@ -0,0 +1,22 @@
+diceplayer:
+    maxcyc: 3
+    opt: no
+    nprocs: 4
+    qmprog: 'g16'
+    lps: no
+    ghosts: no
+    altsteps: 20000
+
+dice:
+    ncores: 3
+    nmol: [1, 50]
+    dens: 0.75
+    nstep: [4000, 6000]
+    isave: 1000
+    ljname: 'phb.ljc'
+    outname: 'phb'
+    randominit: 'first'
+
+gaussian:
+    level: 'MP2/aug-cc-pVDZ'
+    keywords: 'freq=intmodes'

diceplayer/DPpack/Environment/Atom.py

@@ -0,0 +1,52 @@
+from diceplayer.DPpack.Utils.PTable import *
+from diceplayer.DPpack.Utils.Misc import *
+
+
+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/DPpack/Environment/Molecule.py

@@ -0,0 +1,418 @@
+from diceplayer.DPpack.Utils.PTable import *
+from diceplayer.DPpack.Utils.Misc import *
+
+from diceplayer.DPpack.Environment.Atom import Atom
+
+from typing import IO, Any, Final, Tuple, List, TextIO
+from nptyping import Float, NDArray, Shape
+
+from numpy import linalg
+import numpy as np
+
+from copy import deepcopy
+import sys, math
+import sys
+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
+        totalMass (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.totalMass: int = 0
+        self.com: 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
+
+        if a.na == ghost_number:
+
+            self.ghost_atoms.append(self.atom.index(a))
+
+        self.center_of_mass()
+
+    def center_of_mass(self) -> None:
+        """
+        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
+
+    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:
+            sys.exit("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_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:
+
+            sys.exit(
+                "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, fh: TextIO) -> None:
+        """
+        Prints the Molecule information into a Output File
+
+        Args:
+            fh (TextIO): Output File
+        """
+
+        fh.write(
+            "    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()
+
+        fh.write(
+            "    Moments of inertia =  {:>9E}  {:>9E}  {:>9E}\n".format(
+                evals[0], evals[1], evals[2]
+            )
+        )
+
+        fh.write(
+            "    Major principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
+                evecs[0, 0], evecs[1, 0], evecs[2, 0]
+            )
+        )
+        fh.write(
+            "    Inter principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
+                evecs[0, 1], evecs[1, 1], evecs[2, 1]
+            )
+        )
+        fh.write(
+            "    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()
+        fh.write(
+            "    Characteristic lengths = ( {:>6.2f} , {:>6.2f} , {:>6.2f} )\n".format(
+                sizes[0], sizes[1], sizes[2]
+            )
+        )
+        fh.write("    Total mass = {:>8.2f} au\n".format(self.total_mass))
+
+        chg_dip = self.charges_and_dipole()
+        fh.write("    Total charge = {:>8.4f} e\n".format(chg_dip[0]))
+        fh.write(
+            "    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/DPpack/Environment/System.py

@@ -0,0 +1,224 @@
+from diceplayer.DPpack.Utils.PTable import *
+from diceplayer.DPpack.Utils.Misc import *
+
+from diceplayer.DPpack.Environment.Molecule import ANG2BOHR, BOHR2ANG, Molecule
+from diceplayer.DPpack.Environment.Atom import Atom
+
+from typing import IO, Final, Tuple, List, TextIO
+
+from numpy import linalg
+import numpy as np
+
+from copy import deepcopy
+import sys, math
+import sys
+import math
+
+BOHR2ANG: Final[float] = 0.52917721092
+ANG2BOHR: Final[float] = 1 / BOHR2ANG
+
+
+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: Molecule, b: Molecule) -> 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):
+            sys.exit(
+                "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:
+            sys.exit("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":
+            sys.exit("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
+                )
+            )

diceplayer/DPpack/External/Dice.py

@@ -0,0 +1,786 @@
+from dataclasses import dataclass
+from diceplayer.DPpack.Utils.PTable import *
+from diceplayer.DPpack.Utils.Misc import *
+
+from diceplayer.DPpack.Environment.Molecule import Molecule
+from diceplayer.DPpack.Environment.Atom import Atom
+
+from typing import IO, Final, Tuple, List, TextIO, Union
+
+from numpy.core.numeric import partition
+from numpy import random
+
+from multiprocessing import Process, connection
+import subprocess
+import setproctitle
+import os
+import sys
+import shutil
+
+from diceplayer.DPpack.Utils.Validations import NotNull
+
+
+DICE_END_FLAG: Final[str] = "End of simulation"
+DICE_FLAG_LINE: Final[int] = -2
+UMAANG3_TO_GCM3: Final[float] = 1.6605
+
+MAX_SEED: Final[int] = 4294967295
+
+class Dice:
+
+    title = "Diceplayer run"
+    progname = "dice"
+    
+    path = None
+    
+    nprocs: int = None
+    randominit = "first"
+    combrule = "*"
+    ncores = 1
+    
+    temp = 300.0
+    press = 1.0
+    isave = 1000
+    dens = None
+    ljname = None
+    outname = None
+    nmol: List[int] = []
+    nstep: List[int] = []
+    upbuf = 360
+
+    def __init__(self, infile: TextIO, outfile: TextIO) -> None:
+
+        self.infile = infile
+        self.outfile = outfile
+
+    @NotNull(requiredArgs = [
+        "ncores", 
+        "nmol",
+        "dens", 
+        "nstep",
+        "ljname", 
+        "outname"
+    ])
+    def updateKeywords(self, **data):
+        self.__dict__.update(data)
+
+    def __new_density(self, cycle: int, proc: int) -> float:
+
+        sim_dir = "simfiles"
+        step_dir = "step{:02d}".format(cycle - 1)
+        proc_dir = "p{:02d}".format(proc)
+        path = sim_dir + os.sep + step_dir + os.sep + proc_dir
+        file = path + os.sep + "last.xyz"
+
+        if not os.path.isfile(file):
+            sys.exit(
+                "Error: cannot find the xyz file {} in main directory".format(file)
+            )
+        try:
+            with open(file) as fh:
+                xyzfile = fh.readlines()
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        box = xyzfile[1].split()
+        volume = float(box[-3]) * float(box[-2]) * float(box[-1])
+
+        total_mass = 0
+        for i in range(len(self.molecule)):
+
+            total_mass += self.molecule[i].total_mass * self.nmol[i]
+
+        density = (total_mass / volume) * UMAANG3_TO_GCM3
+
+        return density
+
+    def __print_last_config(self, cycle: int, proc: int) -> None:
+
+        sim_dir = "simfiles"
+        step_dir = "step{:02d}".format(cycle)
+        proc_dir = "p{:02d}".format(proc)
+        path = sim_dir + os.sep + step_dir + os.sep + proc_dir
+        file = path + os.sep + "phb.xyz"
+        if not os.path.isfile(file):
+            sys.exit("Error: cannot find the xyz file {}".format(file))
+        try:
+            with open(file) as fh:
+                xyzfile = fh.readlines()
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        nsites = len(self.molecule[0].atom) * self.nmol[0]
+        for i in range(1, len(self.nmol)):
+            nsites += self.nmol[i] * len(self.molecule[i].atom)
+
+        nsites += 2
+
+        nsites *= -1
+        xyzfile = xyzfile[nsites:]
+
+        file = path + os.sep + "last.xyz"
+        fh = open(file, "w")
+        for line in xyzfile:
+            fh.write(line)
+
+    def __make_dice_inputs(self, cycle: int, proc: int) -> None:
+
+        sim_dir = "simfiles"
+        step_dir = "step{:02d}".format(cycle)
+        proc_dir = "p{:02d}".format(proc)
+        path = sim_dir + os.sep + step_dir + os.sep + proc_dir
+
+        num = time.time()
+        num = (num - int(num)) * 1e6
+        
+        num = int((num - int(num)) * 1e6)
+        random.seed((os.getpid() * num) % (MAX_SEED + 1))
+
+        if self.randominit == "first" and cycle > self.initcyc:
+            last_step_dir = "step{:02d}".format(cycle - 1)
+            last_path = sim_dir + os.sep + last_step_dir + os.sep + proc_dir
+            xyzfile = last_path + os.sep + "last.xyz"
+            self.__make_init_file(path, xyzfile)
+
+        if len(self.nstep) == 2:
+
+            self.__make_nvt_ter(cycle, path)
+            self.__make_nvt_eq(path)
+
+        elif len(self.nstep) == 3:
+
+            if self.randominit == "first" and cycle > self.initcyc:
+                self.dens = self.__new_density(cycle, proc)
+            else:
+                self.__make_nvt_ter(cycle, path)
+
+            self.__make_npt_ter(cycle, path)
+            self.__make_npt_eq(path)
+
+        else:
+            sys.exit("Error: bad number of entries for 'nstep'")
+
+        self.__make_potential(path)
+
+    def __make_nvt_ter(self, cycle: int, path: str) -> None:
+
+        file = path + os.sep + "NVT.ter"
+        try:
+            fh = open(file, "w")
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        fh.write("title = {} - NVT Thermalization\n".format(self.title))
+        fh.write("ncores = {}\n".format(self.ncores))
+        fh.write("ljname = {}\n".format(self.ljname))
+        fh.write("outname = {}\n".format(self.outname))
+
+        string = " ".join(str(x) for x in self.nmol)
+        fh.write("nmol = {}\n".format(string))
+
+        fh.write("dens = {}\n".format(self.dens))
+        fh.write("temp = {}\n".format(self.temp))
+
+        if self.randominit == "first" and cycle > self.initcyc:
+            fh.write("init = yesreadxyz\n")
+            fh.write("nstep = {}\n".format(self.altsteps))
+        else:
+            fh.write("init = yes\n")
+            fh.write("nstep = {}\n".format(self.nstep[0]))
+
+        fh.write("vstep = 0\n")
+        fh.write("mstop = 1\n")
+        fh.write("accum = no\n")
+        fh.write("iprint = 1\n")
+        fh.write("isave = 0\n")
+        fh.write("irdf = 0\n")
+
+        seed = int(1e6 * random.random())
+        fh.write("seed = {}\n".format(seed))
+        fh.write("upbuf = {}".format(self.upbuf))
+
+        fh.close()
+
+    def __make_nvt_eq(self, path: str) -> None:
+
+        file = path + os.sep + "NVT.eq"
+        try:
+            fh = open(file, "w")
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        fh.write("title = {} - NVT Production\n".format(self.title))
+        fh.write("ncores = {}\n".format(self.ncores))
+        fh.write("ljname = {}\n".format(self.ljname))
+        fh.write("outname = {}\n".format(self.outname))
+
+        string = " ".join(str(x) for x in self.nmol)
+        fh.write("nmol = {}\n".format(string))
+
+        fh.write("dens = {}\n".format(self.dens))
+        fh.write("temp = {}\n".format(self.temp))
+        fh.write("init = no\n")
+        fh.write("nstep = {}\n".format(self.nstep[1]))
+        fh.write("vstep = 0\n")
+        fh.write("mstop = 1\n")
+        fh.write("accum = no\n")
+        fh.write("iprint = 1\n")
+        fh.write("isave = {}\n".format(self.isave))
+        fh.write("irdf = {}\n".format(10 * self.nprocs))
+
+        seed = int(1e6 * random.random())
+        fh.write("seed = {}\n".format(seed))
+
+        fh.close()
+
+    def __make_npt_ter(self, cycle: int, path: str) -> None:
+
+        file = path + os.sep + "NPT.ter"
+        try:
+            fh = open(file, "w")
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        fh.write("title = {} - NPT Thermalization\n".format(self.title))
+        fh.write("ncores = {}\n".format(self.ncores))
+        fh.write("ljname = {}\n".format(self.ljname))
+        fh.write("outname = {}\n".format(self.outname))
+
+        string = " ".join(str(x) for x in self.nmol)
+        fh.write("nmol = {}\n".format(string))
+
+        fh.write("press = {}\n".format(self.press))
+        fh.write("temp = {}\n".format(self.temp))
+
+        if self.randominit == "first" and cycle > self.initcyc:
+            fh.write("init = yesreadxyz\n")
+            fh.write("dens = {:<8.4f}\n".format(self.dens))
+            fh.write("vstep = {}\n".format(int(self.altsteps / 5)))
+        else:
+            fh.write("init = no\n")
+            fh.write("vstep = {}\n".format(int(self.nstep[1] / 5)))
+
+        fh.write("nstep = 5\n")
+        fh.write("mstop = 1\n")
+        fh.write("accum = no\n")
+        fh.write("iprint = 1\n")
+        fh.write("isave = 0\n")
+        fh.write("irdf = 0\n")
+
+        seed = int(1e6 * random.random())
+        fh.write("seed = {}\n".format(seed))
+
+        fh.close()
+
+    def __make_npt_eq(self, path: str) -> None:
+
+        file = path + os.sep + "NPT.eq"
+        try:
+            fh = open(file, "w")
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        fh.write("title = {} - NPT Production\n".format(self.title))
+        fh.write("ncores = {}\n".format(self.ncores))
+        fh.write("ljname = {}\n".format(self.ljname))
+        fh.write("outname = {}\n".format(self.outname))
+
+        string = " ".join(str(x) for x in self.nmol)
+        fh.write("nmol = {}\n".format(string))
+
+        fh.write("press = {}\n".format(self.press))
+        fh.write("temp = {}\n".format(self.temp))
+
+        fh.write("nstep = 5\n")
+
+        fh.write("vstep = {}\n".format(int(self.nstep[2] / 5)))
+        fh.write("init = no\n")
+        fh.write("mstop = 1\n")
+        fh.write("accum = no\n")
+        fh.write("iprint = 1\n")
+        fh.write("isave = {}\n".format(self.isave))
+        fh.write("irdf = {}\n".format(10 * self.nprocs))
+
+        seed = int(1e6 * random.random())
+        fh.write("seed = {}\n".format(seed))
+
+        fh.close()
+
+    def __make_init_file(self, path: str, file: TextIO) -> None:
+
+        if not os.path.isfile(file):
+            sys.exit(
+                "Error: cannot find the xyz file {} in main directory".format(file)
+            )
+        try:
+            with open(file) as fh:
+                xyzfile = fh.readlines()
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        nsites_mm = 0
+        for i in range(1, len(self.nmol)):
+            nsites_mm += self.nmol[i] * len(self.molecule[i].atom)
+
+        nsites_mm *= -1
+        
+        xyzfile = xyzfile[nsites_mm:]
+
+        file = path + os.sep + self.outname + ".xy"
+
+        try:
+            fh = open(file, "w", 1)
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        for atom in self.molecule[0].atom:
+            fh.write(
+                "{:>10.6f}  {:>10.6f}  {:>10.6f}\n".format(atom.rx, atom.ry, atom.rz)
+            )
+
+        for line in xyzfile:
+            atom = line.split()
+            rx = float(atom[1])
+            ry = float(atom[2])
+            rz = float(atom[3])
+            fh.write("{:>10.6f}  {:>10.6f}  {:>10.6f}\n".format(rx, ry, rz))
+
+        fh.write("$end")
+
+        fh.close()
+
+    def __make_potential(self, path: str) -> None:
+
+        fstr = "{:<3d} {:>3d}  {:>10.5f} {:>10.5f} {:>10.5f}  {:>10.6f} {:>9.5f} {:>7.4f}\n"
+
+        file = path + os.sep + self.ljname
+        try:
+            fh = open(file, "w")
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        fh.write("{}\n".format(self.combrule))
+        fh.write("{}\n".format(len(self.nmol)))
+
+        nsites_qm = (
+            len(self.molecule[0].atom)
+            + len(self.molecule[0].ghost_atoms)
+            + len(self.molecule[0].lp_atoms)
+        )
+
+        fh.write("{} {}\n".format(nsites_qm, self.molecule[0].molname))
+        for atom in self.molecule[0].atom:
+            fh.write(
+                fstr.format(
+                    atom.lbl,
+                    atom.na,
+                    atom.rx,
+                    atom.ry,
+                    atom.rz,
+                    atom.chg,
+                    atom.eps,
+                    atom.sig,
+                )
+            )
+
+        ghost_label = self.molecule[0].atom[-1].lbl + 1
+        for i in self.molecule[0].ghost_atoms:
+            fh.write(
+                fstr.format(
+                    ghost_label,
+                    ghost_number,
+                    self.molecule[0].atom[i].rx,
+                    self.molecule[0].atom[i].ry,
+                    self.molecule[0].atom[i].rz,
+                    self.molecule[0].atom[i].chg,
+                    0,
+                    0,
+                )
+            )
+
+        ghost_label += 1
+        for lp in self.molecule[0].lp_atoms:
+            fh.write(
+                fstr.format(
+                    ghost_label,
+                    ghost_number,
+                    lp["rx"],
+                    lp["ry"],
+                    lp["rz"],
+                    lp["chg"],
+                    0,
+                    0,
+                )
+            )
+
+        for mol in self.molecule[1:]:
+            fh.write("{} {}\n".format(len(mol.atom), mol.molname))
+            for atom in mol.atom:
+                fh.write(
+                    fstr.format(
+                        atom.lbl,
+                        atom.na,
+                        atom.rx,
+                        atom.ry,
+                        atom.rz,
+                        atom.chg,
+                        atom.eps,
+                        atom.sig,
+                    )
+                )
+
+    def __make_proc_dir(self, cycle: int, proc: int) -> None:
+
+        sim_dir = "simfiles"
+        step_dir = "step{:02d}".format(cycle)
+        proc_dir = "p{:02d}".format(proc)
+        path = sim_dir + os.sep + step_dir + os.sep + proc_dir
+        try:
+            os.makedirs(path)
+        except:
+            sys.exit("Error: cannot make directory {}".format(path))
+
+    def __run_dice(self, cycle: int, proc: int, fh: str) -> None:
+
+        sim_dir = "simfiles"
+        step_dir = "step{:02d}".format(cycle)
+        proc_dir = "p{:02d}".format(proc)
+
+        try:
+            fh.write(
+                "Simulation process {} initiated with pid {}\n".format(
+                    sim_dir + os.sep + step_dir + os.sep + proc_dir, os.getpid()
+                )
+            )
+
+        except Exception as err:
+            print("I/O error({0}): {1}".format(err))
+
+        path = sim_dir + os.sep + step_dir + os.sep + proc_dir
+        working_dir = os.getcwd()
+        os.chdir(path)
+
+        if len(self.nstep) == 2:
+
+            if self.randominit == "first" and cycle > self.initcyc:
+                string_tmp = "previous"
+            else:
+                string_tmp = "random"
+
+            string = "(from " + string_tmp + " configuration)"
+            fh.write(
+                "p{:02d}> NVT thermalization finished {} on {}\n".format(
+                    proc, string, date_time()
+                )
+            )
+
+            infh = open("NVT.ter")
+            outfh = open("NVT.ter.out", "w")
+
+            if shutil.which("bash") != None:
+                exit_status = subprocess.call(
+                    [
+                        "bash",
+                        "-c",
+                        "exec -a dice-step{}-p{} {} < {} > {}".format(
+                            cycle, proc, self.progname, infh.name, outfh.name
+                        ),
+                    ]
+                )
+            else:
+                exit_status = subprocess.call(
+                    self.progname, stin=infh.name, stout=outfh.name
+                )
+
+            infh.close()
+            outfh.close()
+
+            if os.getppid() == 1:
+                sys.exit()
+
+            if exit_status != 0:
+                sys.exit(
+                    "Dice process step{:02d}-p{:02d} did not exit properly".format(
+                        cycle, proc
+                    )
+                )
+            else:
+                outfh = open("NVT.ter.out")
+                flag = outfh.readlines()[DICE_FLAG_LINE].strip()
+                outfh.close()
+                if flag != DICE_END_FLAG:
+                    sys.exit(
+                        "Dice process step{:02d}-p{:02d} did not exit properly".format(
+                            cycle, proc
+                        )
+                    )
+
+            fh.write(
+                "p{:02d}> NVT production initiated on {}\n".format(proc, date_time())
+            )
+
+            infh = open("NVT.eq")
+            outfh = open("NVT.eq.out", "w")
+
+            if shutil.which("bash") != None:
+                exit_status = subprocess.call(
+                    [
+                        "bash",
+                        "-c",
+                        "exec -a dice-step{}-p{} {} < {} > {}".format(
+                            cycle, proc, self.progname, infh.name, outfh.name
+                        ),
+                    ]
+                )
+            else:
+                exit_status = subprocess.call(
+                    self.progname, stin=infh.name, stout=outfh.name
+                )
+
+            infh.close()
+            outfh.close()
+
+            if os.getppid() == 1:
+                sys.exit()
+
+            if exit_status != 0:
+                sys.exit(
+                    "Dice process step{:02d}-p{:02d} did not exit properly".format(
+                        cycle, proc
+                    )
+                )
+            else:
+                outfh = open("NVT.eq.out")
+                flag = outfh.readlines()[DICE_FLAG_LINE].strip()
+                outfh.close()
+                if flag != DICE_END_FLAG:
+                    sys.exit(
+                        "Dice process step{:02d}-p{:02d} did not exit properly".format(
+                            cycle, proc
+                        )
+                    )
+
+            fh.write(
+                "p{:02d}> ----- NVT production finished on {}\n".format(
+                    proc, date_time()
+                )
+            )
+
+        elif len(self.nstep) == 3:
+            if (
+                self.randominit == "always"
+                or (self.randominit == "first" and cycle == 1)
+                or self.continued
+            ):
+                string = "(from random configuration)"
+                fh.write(
+                    "p{:02d}> NVT thermalization initiated {} on {}\n".format(
+                        proc, string, date_time()
+                    )
+                )
+                infh = open("NVT.ter")
+                outfh = open("NVT.ter.out", "w")
+
+                if shutil.which("bash") != None:
+                    exit_status = subprocess.call(
+                        [
+                            "bash",
+                            "-c",
+                            "exec -a dice-step{}-p{} {} < {} > {}".format(
+                                cycle, proc, self.progname, infh.name, outfh.name
+                            ),
+                        ]
+                    )
+                else:
+                    exit_status = subprocess.call(
+                        self.progname, stin=infh.name, stout=outfh.name
+                    )
+
+                infh.close()
+                outfh.close()
+
+                if os.getppid() == 1:
+                    sys.exit()
+
+                if exit_status != 0:
+                    sys.exit(
+                        "Dice process step{:02d}-p{:02d} did not exit properly".format(
+                            cycle, proc
+                        )
+                    )
+                else:
+                    outfh = open("NVT.ter.out")
+                    flag = outfh.readlines()[DICE_FLAG_LINE].strip()
+                    outfh.close()
+                    if flag != DICE_END_FLAG:
+                        sys.exit(
+                            "Dice process step{:02d}-p{:02d} did not exit properly".format(
+                                cycle, proc
+                            )
+                        )
+
+            if not self.randominit == "always" or (
+                (self.randominit == "first" and cycle > self.initcyc)
+            ):
+                string = " (from previous configuration) "
+            else:
+                string = " "
+            fh.write(
+                "p{:02d}> NPT thermalization finished {} on {}\n".format(
+                    proc, string, date_time()
+                )
+            )
+
+            infh = open("NPT.ter")
+            outfh = open("NPT.ter.out", "w")
+
+            if shutil.which("bash") != None:
+                exit_status = subprocess.call(
+                    [
+                        "bash",
+                        "-c",
+                        "exec -a dice-step{}-p{} {} < {} > {}".format(
+                            cycle, proc, self.progname, infh.name, outfh.name
+                        ),
+                    ]
+                )
+            else:
+                exit_status = subprocess.call(
+                    self.progname, stin=infh.name, stout=outfh.name
+                )
+
+            infh.close()
+            outfh.close()
+
+            if os.getppid() == 1:
+                sys.exit()
+
+            if exit_status != 0:
+                sys.exit(
+                    "Dice process step{:02d}-p{:02d} did not exit properly".format(
+                        cycle, proc
+                    )
+                )
+            else:
+                outfh = open("NPT.ter.out")
+                flag = outfh.readlines()[DICE_FLAG_LINE].strip()
+                outfh.close()
+                if flag != DICE_END_FLAG:
+                    sys.exit(
+                        "Dice process step{:02d}-p{:02d} did not exit properly".format(
+                            cycle, proc
+                        )
+                    )
+
+            fh.write(
+                "p{:02d}> NPT production initiated on {}\n".format(proc, date_time())
+            )
+
+            infh = open("NPT.eq")
+            outfh = open("NPT.eq.out", "w")
+
+            if shutil.which("bash") != None:
+                exit_status = subprocess.call(
+                    [
+                        "bash",
+                        "-c",
+                        "exec -a dice-step{}-p{} {} < {} > {}".format(
+                            cycle, proc, self.progname, infh.name, outfh.name
+                        ),
+                    ]
+                )
+            else:
+                exit_status = subprocess.call(
+                    self.progname, stin=infh.name, stout=outfh.name
+                )
+
+            infh.close()
+            outfh.close()
+
+            if os.getppid() == 1:
+                sys.exit()
+
+            if exit_status != 0:
+                sys.exit(
+                    "Dice process step{:02d}-p{:02d} did not exit properly".format(
+                        cycle, proc
+                    )
+                )
+            else:
+                outfh = open("NPT.eq.out")
+                flag = outfh.readlines()[DICE_FLAG_LINE].strip()
+                outfh.close()
+                if flag != DICE_END_FLAG:
+                    sys.exit(
+                        "Dice process step{:02d}-p{:02d} did not exit properly".format(
+                            cycle, proc
+                        )
+                    )
+
+            fh.write(
+                "p{:02d}> ----- NPT production finished on {}\n".format(
+                    proc, date_time()
+                )
+            )
+
+        os.chdir(working_dir)
+
+    def __simulation_process(self, cycle: int, proc: int):
+        setproctitle.setproctitle("diceplayer-step{:0d}-p{:0d}".format(cycle, proc))
+
+        try:
+            self.__make_proc_dir(cycle, proc)
+            self.__make_dice_inputs(cycle, proc)
+            self.__run_dice(cycle, proc, self.outfile)
+        except Exception as err:
+            sys.exit(err)
+
+    def configure(
+        self,
+        initcyc: int,
+        nprocs: int,
+        altsteps: int,
+        nmol: List[int],
+        molecule: List[Molecule],
+    ):
+
+        self.initcyc = initcyc
+
+        self.nprocs = nprocs
+        self.altsteps = altsteps
+        self.nmol = nmol
+        self.molecule = molecule
+
+    def start(self, cycle: int) -> None:
+
+        procs = []
+        sentinels = []
+
+        for proc in range(1, self.nprocs + 1):
+
+            p = Process(target=self.__simulation_process, args=(cycle, proc))
+            p.start()
+
+            procs.append(p)
+            sentinels.append(p.sentinel)
+
+        while procs:
+            finished = connection.wait(sentinels)
+            for proc_sentinel in finished:
+                i = sentinels.index(proc_sentinel)
+                status = procs[i].exitcode
+                procs.pop(i)
+                sentinels.pop(i)
+                if status != 0:
+                    for p in procs:
+                        p.terminate()
+                    sys.exit(status)
+
+        for proc in range(1, self.nprocs + 1):
+            self.__print_last_config(cycle, proc)
+
+    def reset(self):
+
+        del self.nprocs
+        del self.altsteps
+        del self.molecule

diceplayer/DPpack/External/Gaussian.py

@@ -0,0 +1,688 @@
+from turtle import position
+from diceplayer.DPpack.Environment.Atom import *
+from diceplayer.DPpack.Utils.PTable import *
+from diceplayer.DPpack.Utils.Misc import *
+
+from diceplayer.DPpack.External.Dice import *
+
+from diceplayer.DPpack.Environment.Molecule import Molecule
+from diceplayer.DPpack.Environment.Atom import Atom
+
+from typing import IO, List, TextIO
+
+import numpy as np
+
+import subprocess
+import os
+import sys
+import shutil
+import textwrap
+import types
+
+
+class Gaussian:
+
+    qmprog = "g09"
+    path = None
+    mem = None
+    keywords = None
+    chgmult = [0, 1]
+    gmiddle = None  # In each case, if a filename is given, its content will be
+    gbottom = None  # inserted in the gaussian input
+    pop = "chelpg"
+    level = None
+
+    def __init__(self) -> None:
+        pass
+
+    @NotNull(requiredArgs = [
+        "level"
+    ])
+    def updateKeywords(self, **data):
+        self.__dict__.update(data)
+
+    def run_formchk(self, cycle: int, fh: TextIO):
+
+        simdir = "simfiles"
+        stepdir = "step{:02d}".format(cycle)
+        path = simdir + os.sep + stepdir + os.sep + "qm"
+
+        work_dir = os.getcwd()
+        os.chdir(path)
+
+        fh.write("Formatting the checkpoint file... \n")
+
+        exit_status = subprocess.call(["formchk", "asec.chk"], stdout=fh)
+
+        fh.write("Done\n")
+
+        os.chdir(work_dir)
+
+    def read_forces_fchk(self, file: str, fh: TextIO) -> np.ndarray:
+
+        forces = []
+        try:
+            with open(file) as tmpfh:
+                fchkfile = tmpfh.readlines()
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        start = fchkfile.pop(0).strip()
+        while start.find("Cartesian Gradient") != 0:  # expression in begining of line
+            start = fchkfile.pop(0).strip()
+
+        degrees = 3 * len(self.molecule[0].atom)
+        count = 0
+        while len(forces) < degrees:
+            values = fchkfile.pop(0).split()
+            forces.extend([float(x) for x in values])
+            count += len(values)
+            if count >= degrees:
+                forces = forces[:degrees]
+                break
+
+        gradient = np.array(forces)
+
+        fh.write("\nGradient read from file {}:\n".format(file))
+        fh.write(
+            "-----------------------------------------------------------------------\n"
+            "Center     Atomic                     Forces (Hartree/Bohr)\n"
+            "Number     Number              X                Y                Z\n"
+            "-----------------------------------------------------------------------\n"
+        )
+        for i in range(len(self.molecule[0].atom)):
+            fh.write(
+                "  {:>5d}     {:>3d}        {:>14.9f}   {:>14.9f}   {:>14.9f}\n".format(
+                    i + 1,
+                    self.molecule[0].atom[i].na,
+                    forces.pop(0),
+                    forces.pop(0),
+                    forces.pop(0),
+                )
+            )
+
+        fh.write(
+            "-----------------------------------------------------------------------\n"
+        )
+
+        force_max = np.amax(np.absolute(gradient))
+        force_rms = np.sqrt(np.mean(np.square(gradient)))
+
+        fh.write(
+            "  Max Force = {:>14.9f}      RMS Force = {:>14.9f}\n\n".format(
+                force_max, force_rms
+            )
+        )
+
+        return gradient
+
+    def read_hessian_fchk(self, file: str) -> np.ndarray:
+
+        force_const = []
+        try:
+            with open(file) as tmpfh:
+                fchkfile = tmpfh.readlines()
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        start = fchkfile.pop(0).strip()
+        while start.find("Cartesian Force Constants") != 0:
+            start = fchkfile.pop(0).strip()
+
+        degrees = 3 * len(self.molecule[0].atom)
+        last = round(degrees * (degrees + 1) / 2)
+        count = 0
+
+        while len(force_const) < last:
+
+            value = fchkfile.pop(0).split()
+            force_const.extend([float(x) for x in value])
+
+        # while len(force_const) < last:
+
+        # 	values = fchkfile.pop(0).split()
+        # 	force_const.extend([ float(x) for x in values ])
+        # 	count += len(values)
+        # 	if count >= last:
+        # 		force_const = force_const[:last]
+        # 		break
+
+        hessian = np.zeros((degrees, degrees))
+        for i in range(degrees):
+            for j in range(i + 1):
+                hessian[i, j] = force_const.pop(0)
+                hessian[j, i] = hessian[i, j]
+
+        return hessian
+
+    def read_hessian_log(self, file: str) -> np.ndarray:
+
+        try:
+            with open(file) as tmpfh:
+                logfile = tmpfh.readlines()
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        start = logfile.pop(0).strip()
+        while start.find("The second derivative matrix:") != 0:
+            start = logfile.pop(0).strip()
+
+        degrees = 3 * len(self.molecule[0].atom)
+        hessian = np.zeros((degrees, degrees))
+
+        k = 0
+        while k < degrees:
+            logfile.pop(0)
+            for i in range(k, degrees):
+                values = logfile.pop(0).split()[1:]
+                for j in range(k, min(i + 1, k + 5)):
+                    hessian[i, j] = float(values.pop(0))
+                    hessian[j, i] = hessian[i, j]
+            k += 5
+
+        return hessian
+
+    def print_grad_hessian(
+        self, cycle: int, cur_gradient: np.ndarray, hessian: np.ndarray
+    ) -> None:
+
+        try:
+            fh = open("grad_hessian.dat", "w")
+        except:
+            sys.exit("Error: cannot open file grad_hessian.dat")
+
+        fh.write("Optimization cycle: {}\n".format(cycle))
+        fh.write("Cartesian Gradient\n")
+        degrees = 3 * len(self.molecule[0].atom)
+        for i in range(degrees):
+            fh.write(" {:>11.8g}".format(cur_gradient[i]))
+            if (i + 1) % 5 == 0 or i == degrees - 1:
+                fh.write("\n")
+
+        fh.write("Cartesian Force Constants\n")
+        n = int(np.sqrt(2 * degrees))
+        last = degrees * (degrees + 1) / 2
+        count = 0
+        for i in range(n):
+            for j in range(i + 1):
+                count += 1
+                fh.write(" {:>11.8g}".format(hessian[i, j]))
+                if count % 5 == 0 or count == last:
+                    fh.write("\n")
+
+        fh.close()
+
+    # Change the name to make_gaussian_input
+    def make_gaussian_input(self, cycle: int, asec_charges=None) -> None:
+
+        simdir = "simfiles"
+        stepdir = "step{:02d}".format(cycle)
+        path = simdir + os.sep + stepdir + os.sep + "qm"
+
+        file = path + os.sep + "asec.gjf"
+
+        try:
+            fh = open(file, "w")
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        fh.write("%Chk=asec.chk\n")
+        if self.mem != None:
+            fh.write("%Mem={}MB\n".format(self.mem))
+        fh.write("%Nprocs={}\n".format(self.nprocs * self.ncores))
+
+        kword_line = "#P " + str(self.level)
+
+        if self.keywords != None:
+            kword_line += " " + self.keywords
+
+        if self.opt == "yes":
+            kword_line += " Force"
+
+        # kword_line += " Charge"
+        kword_line += " NoSymm"
+        kword_line += " Pop={} Density=Current".format(self.pop)
+
+        if cycle > 1:
+            kword_line += " Guess=Read"
+
+        fh.write(textwrap.fill(kword_line, 90))
+        fh.write("\n")
+
+        fh.write("\nForce calculation - Cycle number {}\n".format(cycle))
+        fh.write("\n")
+        fh.write("{},{}\n".format(self.chgmult[0], self.chgmult[1]))
+
+        for atom in self.molecule[0].atom:
+            symbol = atomsymb[atom.na]
+            fh.write(
+                "{:<2s}    {:>10.5f}   {:>10.5f}   {:>10.5f}\n".format(
+                    symbol, atom.rx, atom.ry, atom.rz
+                )
+            )
+
+        # ## If also performing charge fit in the same calculation
+        # if cycle >= self.player.switchcyc:
+        # 	for ghost in ghost_atoms:
+        # 		fh.write("Bq    {:>10.5f}   {:>10.5f}   {:>10.5f}\n".format(
+        # 											ghost['rx'], ghost['ry'], ghost['rz']))
+
+        # 	for lp in lp_atoms:
+        # 		fh.write("Bq    {:>10.5f}   {:>10.5f}   {:>10.5f}\n".format(
+        # 														lp['rx'], lp['ry'], lp['rz']))
+
+        # fh.write("\n")
+
+        # If gmiddle file was informed, write its contents in asec.gjf
+        # if self.gmiddle != None:
+        # 	if not os.path.isfile(self.gmiddle):
+        # 		sys.exit("Error: cannot find file {} in main directory".format(
+        # 																self.gmiddle))
+        # 	try:
+        # 		with open(self.gmiddle) as gmiddlefile:
+        # 			gmiddle = gmiddlefile.readlines()
+        # 	except:
+        # 		sys.exit("Error: cannot open file {}".format(self.gmiddle))
+
+        # 	for line in gmiddle:
+        # 		fh.write(line)
+
+        # 	fh.write("\n")
+
+        # ## Write the ASEC:
+        # for charge in asec_charges:
+        # 	fh.write("{:>10.5f}   {:>10.5f}   {:>10.5f}     {:>11.8f}\n".format(
+        # 							charge['rx'], charge['ry'], charge['rz'], charge['chg']))
+
+        fh.write("\n")
+
+        # ## If gbottom file was informed, write its contents in asec.gjf
+        # if self.gbottom != None:
+        # 	if not os.path.isfile(self.gbottom):
+        # 		sys.exit("Error: cannot find file {} in main directory".format(
+        # 																self.gbottom))
+        # 	try:
+        # 		with open(self.gbottom) as gbottomfile:
+        # 			gbottom = gbottomfile.readlines()
+        # 	except:
+        # 		sys.exit("Error: cannot open file {}".format(self.gbottom))
+
+        # 	for line in gbottom:
+        # 		fh.write(line)
+
+        # fh.write("\n")
+
+        # fh.close()
+
+    def read_charges(self, file: str, fh: TextIO) -> None:
+
+        try:
+            with open(file) as tmpfh:
+                glogfile = tmpfh.readlines()
+        except:
+            sys.exit("Error: cannot open file {}".format(file))
+
+        start = glogfile.pop(0).strip()
+        while start != "Fitting point charges to electrostatic potential":
+            start = glogfile.pop(0).strip()
+
+        glogfile = glogfile[3:]  # Consume 3 more lines
+
+        fh.write("\nAtomic charges:\n")
+        fh.write("------------------------------------\n")
+        for atom in self.molecule[0].atom:
+            line = glogfile.pop(0).split()
+            atom_str = line[1]
+            charge = float(line[2])
+            atom.chg = charge
+            fh.write(" {:<2s}      {:>10.6f}\n".format(atom_str, charge))
+
+        # if self.pop == "chelpg":
+        # 	for ghost in ghost_atoms:
+        # 		line = glogfile.pop(0).split()
+        # 		atom_str = line[1]
+        # 		charge = float(line[2])
+        # 		ghost['chg'] = charge
+        # 		fh.write(" {:<2s}      {:>10.6f}\n".format(atom_str, charge))
+
+        # 	for lp in lp_atoms:
+        # 		line = glogfile.pop(0).split()
+        # 		atom_str = line[1]
+        # 		charge = float(line[2])
+        # 		lp['chg'] = charge
+        # 		fh.write(" {:<2s}      {:>10.6f}\n".format(atom_str, charge))
+
+        fh.write("------------------------------------\n")
+
+    def run_gaussian(self, cycle: int, type: str, fh: TextIO) -> None:
+
+        simdir = "simfiles"
+        stepdir = "step{:02d}".format(cycle)
+        path = simdir + os.sep + stepdir + os.sep + "qm"
+        work_dir = os.getcwd()
+        os.chdir(path)
+
+        # if type == "force":
+        # 	infile = "asec.gjf"
+        # elif type == "charge":
+        # 	infile = "asec2.gjf"
+
+        infile = "asec.gjf"
+
+        fh.write(
+            "\nCalculation of {}s initiated with Gaussian on {}\n".format(
+                type, date_time()
+            )
+        )
+
+        if shutil.which("bash") != None:
+            exit_status = subprocess.call(
+                [
+                    "bash",
+                    "-c",
+                    "exec -a {}-step{} {} {}".format(
+                        self.qmprog, cycle, self.qmprog, infile
+                    ),
+                ]
+            )
+        else:
+            exit_status = subprocess.call([self.qmprog, infile])
+
+        if exit_status != 0:
+            sys.exit("Gaussian process did not exit properly")
+
+        fh.write("Calculation of {}s finished on {}\n".format(type, date_time()))
+
+        os.chdir(work_dir)
+
+    # def calculate_step(
+    #     self, cycle: int, gradient: np.ndarray, hessian: np.ndarray
+    # ) -> np.ndarray:
+
+    #     invhessian = np.linalg.inv(hessian)
+    #     pre_step = -1 * np.matmul(invhessian, gradient.T).T
+    #     maxstep = np.amax(np.absolute(pre_step))
+    #     factor = min(1, self.player.maxstep / maxstep)
+    #     step = factor * pre_step
+
+    #     self.outfile.write("\nCalculated step-{}:\n".format(cycle))
+    #     pre_step_list = pre_step.tolist()
+
+    #     self.outfile.write(
+    #         "-----------------------------------------------------------------------\n"
+    #         "Center     Atomic                          Step (Bohr)\n"
+    #         "Number     Number              X                Y                Z\n"
+    #         "-----------------------------------------------------------------------\n"
+    #     )
+    #     for i in range(len(self.system.molecule[0].atom)):
+    #         self.outfile.write(
+    #             "  {:>5d}     {:>3d}        {:>14.9f}   {:>14.9f}   {:>14.9f}\n".format(
+    #                 i + 1,
+    #                 self.system.molecule[0].atom[i].na,
+    #                 pre_step_list.pop(0),
+    #                 pre_step_list.pop(0),
+    #                 pre_step_list.pop(0),
+    #             )
+    #         )
+
+    #     self.outfile.write(
+    #         "-----------------------------------------------------------------------\n"
+    #     )
+
+    #     self.outfile.write("Maximum step is {:>11.6}\n".format(maxstep))
+    #     self.outfile.write("Scaling factor = {:>6.4f}\n".format(factor))
+    #     self.outfile.write("\nFinal step (Bohr):\n")
+    #     step_list = step.tolist()
+
+    #     self.outfile.write(
+    #         "-----------------------------------------------------------------------\n"
+    #         "Center     Atomic                          Step (Bohr)\n"
+    #         "Number     Number              X                Y                Z\n"
+    #         "-----------------------------------------------------------------------\n"
+    #     )
+    #     for i in range(len(self.system.molecule[0].atom)):
+    #         self.outfile.write(
+    #             "  {:>5d}     {:>3d}        {:>14.9f}   {:>14.9f}   {:>14.9f}\n".format(
+    #                 i + 1,
+    #                 self.system.molecule[0].atom[i].na,
+    #                 step_list.pop(0),
+    #                 step_list.pop(0),
+    #                 step_list.pop(0),
+    #             )
+    #         )
+
+    #     self.outfile.write(
+    #         "-----------------------------------------------------------------------\n"
+    #     )
+
+    #     step_max = np.amax(np.absolute(step))
+    #     step_rms = np.sqrt(np.mean(np.square(step)))
+
+    #     self.outfile.write(
+    #         "  Max Step = {:>14.9f}      RMS Step = {:>14.9f}\n\n".format(
+    #             step_max, step_rms
+    #         )
+    #     )
+
+    #     return step
+
+
+    def configure(
+        self,
+        initcyc: int,
+        nprocs: int,
+        ncores: int,
+        altsteps: int,
+        switchcyc: int,
+        opt: str,
+        nmol: List[int],
+        molecule: List[Molecule],
+    ):
+
+        self.initcyc = initcyc
+
+        self.nprocs = nprocs
+        self.ncores = ncores
+        self.altsteps = altsteps
+        self.switchcyc = switchcyc
+        self.opt = opt
+
+        self.molecule = molecule
+
+    def start(self, cycle: int, outfile: TextIO, readhessian: str) -> np.ndarray:
+
+        make_qm_dir(cycle)
+
+        if self.qmprog in ("g03", "g09", "g16"):
+
+            if cycle > 1:
+
+                src = (
+                    "simfiles"
+                    + os.sep
+                    + "step{:02d}".format(cycle - 1)
+                    + os.sep
+                    + "qm"
+                    + os.sep
+                    + "asec.chk"
+                )
+                dst = (
+                    "simfiles"
+                    + os.sep
+                    + "step{:02d}".format(cycle)
+                    + os.sep
+                    + "qm"
+                    + os.sep
+                    + "asec.chk"
+                )
+                shutil.copyfile(src, dst)
+
+            self.make_gaussian_input(cycle)
+            self.run_gaussian(cycle, "force", outfile)
+            self.run_formchk(cycle, outfile)
+
+            ## Read the gradient
+            file = (
+                "simfiles"
+                + os.sep
+                + "step{:02d}".format(cycle)
+                + os.sep
+                + "qm"
+                + os.sep
+                + "asec.fchk"
+            )
+
+            try:
+                gradient
+                old_gradient = gradient
+            except:
+                pass
+
+            gradient = self.read_forces_fchk(file, outfile)
+
+            # If 1st step, read the hessian
+            if cycle == 1:
+
+                if readhessian == "yes":
+
+                    file = "grad_hessian.dat"
+                    outfile.write(
+                        "\nReading the hessian matrix from file {}\n".format(file)
+                    )
+                    hessian = self.read_hessian_log(file)
+
+                else:
+
+                    file = (
+                        "simfiles"
+                        + os.sep
+                        + "step01"
+                        + os.sep
+                        + "qm"
+                        + os.sep
+                        + "asec.fchk"
+                    )
+                    outfile.write(
+                        "\nReading the hessian matrix from file {}\n".format(file)
+                    )
+                    hessian = self.read_hessian_fchk(file)
+
+            # From 2nd step on, update the hessian
+            else:
+                outfile.write(
+                    "\nUpdating the hessian matrix using the BFGS method... "
+                )
+                hessian = self.molecule[0].update_hessian(
+                    step, gradient, old_gradient, hessian
+                )
+                outfile.write("Done\n")
+
+            # Save gradient and hessian
+            self.print_grad_hessian(cycle, gradient, hessian)
+
+            # Calculate the step and update the position
+            step = self.calculate_step(cycle, gradient, hessian)
+
+            position += step
+
+            ## If needed, calculate the charges
+            if cycle < self.switchcyc:
+
+                # internal.gaussian.make_charge_input(cycle, asec_charges)
+                self.run_gaussian(cycle, "charge", outfile)
+
+                file = (
+                    "simfiles"
+                    + os.sep
+                    + "step{:02d}".format(cycle)
+                    + os.sep
+                    + "qm"
+                    + os.sep
+                    + "asec2.log"
+                )
+                self.read_charges(file, outfile)
+            else:
+                file = (
+                    "simfiles"
+                    + os.sep
+                    + "step{:02d}".format(cycle)
+                    + os.sep
+                    + "qm"
+                    + os.sep
+                    + "asec.log"
+                )
+                self.read_charges(file, outfile)
+
+            ## Print new info for molecule[0]
+            self.outfile.write("\nNew values for molecule type 1:\n\n")
+            self.molecule[0].print_mol_info(outfile)
+
+            ##
+            ##  Molcas block
+            ##
+            # if player['qmprog'] == "molcas":
+
+        # elif player['opt'] == "ts":
+
+        ##
+        ##  Gaussian block
+        ##
+        # if player['qmprog'] in ("g03", "g09", "g16"):
+
+        ##
+        ##  Molcas block
+        ##
+        # if player['qmprog'] == "molcas":
+
+        # else:  ## Only relax the charge distribution
+
+        #     if internal.player.qmprog in ("g03", "g09", "g16"):
+
+        #         if cycle > 1:
+        #             src = (
+        #                 "simfiles"
+        #                 + os.sep
+        #                 + "step{:02d}".format(cycle - 1)
+        #                 + os.sep
+        #                 + "qm"
+        #                 + os.sep
+        #                 + "asec.chk"
+        #             )
+        #             dst = (
+        #                 "simfiles"
+        #                 + os.sep
+        #                 + "step{:02d}".format(cycle)
+        #                 + os.sep
+        #                 + "qm"
+        #                 + os.sep
+        #                 + "asec.chk"
+        #             )
+        #             shutil.copyfile(src, dst)
+
+        #         # internal.gaussian.make_charge_input(cycle, asec_charges)
+        #         internal.gaussian.run_gaussian(cycle, "charge", internal.outfile)
+
+        #         file = (
+        #             "simfiles"
+        #             + os.sep
+        #             + "step{:02d}".format(cycle)
+        #             + os.sep
+        #             + "qm"
+        #             + os.sep
+        #             + "asec2.log"
+        #         )
+        #         internal.read_charges(file)
+
+        #         ## Print new info for molecule[0]
+        #         internal.outfile.write("\nNew values for molecule type 1:\n\n")
+        #         internal.system.molecule[0].print_mol_info()
+
+        # 			#if player['qmprog'] == "molcas"
+
+        return position
+
+    def reset(self):
+
+        del self.nprocs
+        del self.altsteps
+        del self.molecule

diceplayer/DPpack/MolHandling.py

@@ -1,460 +0,0 @@
-from diceplayer.DPpack.PTable import *
-from diceplayer.DPpack.Misc import *
-
-from typing import IO, Tuple, List, TextIO, Union
-
-from numpy import linalg
-import numpy as np
-
-from copy import deepcopy
-import sys, math
-import textwrap
-import os, sys
-import shutil
-import math
-
-env = ["OMP_STACKSIZE"]
-
-bohr2ang = 0.52917721092
-ang2bohr = 1/bohr2ang
-
-class Atom:
-
-	def __init__(self, lbl: int, na: int, 
-				 rx: float,ry: float, rz: float,
-				 chg: float, eps: float, sig: float
-				 ) -> None:
-		
-		self.lbl = lbl                  # Integer
-		self.na = na                    # Integer
-		self.rx = rx                    # Double
-		self.ry = ry                    # Double
-		self.rz = rz                    # Double
-		self.chg = chg                  # Double
-		self.eps = eps                  # Double
-		self.sig = sig                  # Double
-		self.mass = atommass[self.na]   # Double
-
-# Classe que conterá toda informação e funções relacionadas a uma unica molecula
-
-class Molecule:
-
-	def __init__(self, molname: str ) -> None:
-
-		self.molname = molname
-
-		self.atom: List[Atom] = []					# Lista de instancias de Atom
-		self.position = None			# Array Numpy
-		self.energy = None				# Array Numpy
-		self.gradient = None			# Array Numpy
-		self.hessian = None				# Array Numpy
-
-		self.total_mass = 0
-		self.com = None
-		self.ghost_atoms: List[int] = [] # Stores the index of the ghost atoms in the atoms array
-		self.lp_atoms: List[int] = []
-
-	def add_atom(self, a: Atom) -> None:
-
-		self.atom.append(a)			# Inserção de um novo atomo
-		self.total_mass += a.mass
-
-		if (a.na == ghost_number):
-
-			self.ghost_atoms.append(self.atom.index(a))
-
-		self.center_of_mass()
-
-	def center_of_mass(self) -> None:
-	
-		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
-
-	def center_of_mass_to_origin(self)-> None:
-	
-		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]:
-	
-		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) -> np.ndarray:
-		
-		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) -> np.ndarray:
-		
-		self.center_of_mass()
-		Ixx = Ixy = Ixz = Iyy = Iyz = Izz = 0.0
-
-		for atom in self.atom:
-
-			####  Obtain the displacement from the center of mass
-			dx = atom.rx - self.com[0]
-			dy = atom.ry - self.com[1]
-			dz = atom.rz - self.com[2]
-			####  Update the diagonal components of the tensor
-			Ixx += atom.mass * (dy**2 + dz**2)
-			Iyy += atom.mass * (dz**2 + dx**2)
-			Izz += atom.mass * (dx**2 + dy**2)
-			####  Update the off-diagonal components of the tensor
-			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 eixos(self) -> np.ndarray:
-	
-		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]:
-		
-		try:
-			evals, evecs = linalg.eigh(self.inertia_tensor())
-		except:
-			sys.exit("Error: diagonalization of inertia tensor did not converge")
-	
-		return evals, evecs
-
-	def read_position(self) -> np.ndarray:
-	
-		position_list = []
-		for atom in self.atom:
-			position_list.extend([ atom.rx, atom.ry, atom.rz ])
-		position = np.array(position_list)
-		position *= ang2bohr
-		
-		return position
-
-	def update_hessian(self, step: np.ndarray,
-					   cur_gradient: np.ndarray,
-					   old_gradient: np.ndarray,
-					   hessian: np.ndarray
-					   ) -> None:	## According to the BFGS
-		
-		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]:
-	
-		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:
-		
-		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:
-			
-			sys.exit("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':
-	
-		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, fh: TextIO) -> None:
-
-		fh.write("    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()
-		
-		fh.write("    Moments of inertia =  {:>9E}  {:>9E}  {:>9E}\n".format(evals[0], 
-																		evals[1], evals[2]))
-		
-		fh.write("    Major principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
-														evecs[0,0], evecs[1,0], evecs[2,0]))
-		fh.write("    Inter principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
-														evecs[0,1], evecs[1,1], evecs[2,1]))
-		fh.write("    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()
-		fh.write("    Characteristic lengths = ( {:>6.2f} , {:>6.2f} , {:>6.2f} )\n".format(
-																sizes[0], sizes[1], sizes[2]))
-		fh.write("    Total mass = {:>8.2f} au\n".format(self.total_mass))
-		
-		chg_dip = self.charges_and_dipole()
-		fh.write("    Total charge = {:>8.4f} e\n".format(chg_dip[0]))
-		fh.write("    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') -> List[float]:
-	
-		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)
-
-# Usaremos uma nova classe que ira conter toda interação entre moleculas
-
-class System:
-
-	def __init__(self) -> None:
-
-		self.molecule: List[Molecule] = []
-		self.nmols: List[int] = []
-
-	def add_type(self, nmols: int, m: Molecule) -> None:
-		
-		self.molecule.append(m)
-		self.nmols.append(nmols)
-
-	# Função que calcula a distância entre dois centros de massa
-	# e por se tratar de uma função de dois atomos não deve ser
-	# inserida dentro de Molecule
-	def center_of_mass_distance(self, a: Molecule, b: Molecule) -> float:
-	
-		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):
-			sys.exit("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:
-			sys.exit("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:
-	
-		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":
-			sys.exit("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:
-	
-		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))
-
-
-	

diceplayer/DPpack/Molcas.py

@@ -1,348 +0,0 @@
-import os, sys
-import textwrap
-import subprocess
-
-from DPpack.PTable import *
-from DPpack.SetGlobals import *
-from DPpack.MolHandling import *
-from DPpack.Misc import *
-
-#######################################  functions  ######################################
-
-def read_forces_log(file, fh):
-	
-	forces = []
-	try:
-		with open(file) as tmpfh:
-			logfile = tmpfh.readlines()
-	except:
-		sys.exit("Error: cannot open file {}".format(file))
-	
-	start = logfile.pop(0).strip()
-	while  start.find("Molecular gradients") < 0:	##  expression not found
-		start = logfile.pop(0).strip()
-	
-	logfile = logfile[7:]		##  skip next 7 lines
-	
-	for i in range(len(molecules[0])):
-		values = logfile.pop(0).split()
-		values = values[1:]
-		forces.extend([ float(x) for x in values ])
-	
-	gradient = np.array(forces)
-	
-	fh.write("\nGradient read from file {}:\n".format(file))
-	fh.write("-----------------------------------------------------------------------\n"
-			 "Center     Atomic                     Forces (Hartree/Bohr)\n"
-			 "Number     Number              X                Y                Z\n"
-			 "-----------------------------------------------------------------------\n")
-	for i in range(len(molecules[0])):
-		fh.write("  {:>5d}     {:>3d}        {:>14.9f}   {:>14.9f}   {:>14.9f}\n".format(
-			   i + 1, molecules[0][i]['na'], forces.pop(0), forces.pop(0), forces.pop(0)))
-	
-	fh.write("-----------------------------------------------------------------------\n")
-	
-	force_max = np.amax(np.absolute(gradient))
-	force_rms = np.sqrt(np.mean(np.square(gradient)))
-	
-	fh.write("  Max Force = {:>14.9f}      RMS Force = {:>14.9f}\n\n".format(
-																	force_max, force_rms))
-	
-	return gradient
-
-
-
-def read_hessian_log(file):
-	
-	force_const = []
-	try:
-		with open(file) as tmpfh:
-			logfile = tmpfh.readlines()
-	except:
-		sys.exit("Error: cannot open file {}".format(file))
-	
-	start = logfile.pop(0).strip()
-	while  start.find("Force constant matrix") < 0:
-		start = logfile.pop(0).strip()
-	
-	logfile = logfile[1:]		## skip next 1 line
-	
-	degrees = 3 * len(molecules[0])
-	dim = degrees
-	last = round(dim * dim)
-	count = 0
-	while True:
-		values = logfile.pop(0).rstrip()
-		while len(values) != 0:
-			new_value = values[:16]
-			values = values[16:]
-			force_const.append(float(new_value))
-			count += 1
-		
-		if count >= last:
-			break
-	
-	hessian = np.array(force_const).reshape(dim, dim)
-	hessian = hessian[:degrees, :degrees]		## remove degrees related to ghost atoms
-	
-	for i in range(degrees):
-		for j in range(i + 1):
-			hessian[j,i] = hessian[i,j]		## force the hessian to be symmetric
-	
-	return hessian
-
-
-
-def print_grad_hessian(cycle, cur_gradient, hessian):
-	
-	try:
-		fh = open("grad_hessian.dat", "w")
-	except:
-		sys.exit("Error: cannot open file grad_hessian.dat")
-	
-	fh.write("Optimization cycle: {}\n".format(cycle))
-	fh.write("Cartesian Gradient\n")
-	degrees = 3 * len(molecules[0])
-	for i in range(degrees):
-		fh.write(" {:>11.8g}".format(cur_gradient[i]))
-		if (i + 1) % 5 == 0 or i == degrees - 1:
-			fh.write("\n")
-	
-	fh.write("Cartesian Force Constants\n")
-	last = degrees * (degrees + 1) / 2
-	count = 0
-	for i in range(degrees):
-		for j in range(i + 1):
-			count += 1
-			fh.write(" {:>11.8g}".format(hessian[i,j]))
-			if count % 5 == 0 or count == last:
-				fh.write("\n")
-	
-	fh.close()
-	
-	return
-
-
-
-def make_asec_file(cycle, asec_charges):
-	
-	path = "step{:02d}".format(cycle) + os.sep + "qm"
-	file = path + os.sep + "asec.xfield"
-	try:
-		fh = open(file, "w")
-	except:
-		sys.exit("Error: cannot open file {}".format(file))
-	
-	fh.write("{} Angstrom\n".format(len(asec_charges)))
-	
-	## Write the ASEC:
-	for charge in asec_charges:
-		fh.write("{:>10.5f}   {:>10.5f}   {:>10.5f}    {:>11.8f}    0.0  0.0  0.0\n".format(
-								 charge['rx'], charge['ry'], charge['rz'], charge['chg']))
-	
-	fh.write("End of input\n")
-	fh.close()
-	
-	return
-
-
-
-def make_force_input(cycle, asec_charges):
-	
-	path = "step{:02d}".format(cycle) + os.sep + "qm"
-	file = path + os.sep + "asec.input"
-	try:
-		fh = open(file, "w")
-	except:
-		sys.exit("Error: cannot open file {}".format(file))
-	
-	fh.write(" &Gateway\n")
-	fh.write("  Coord\n")
-	
-	nsites = len(molecules[0])
-	if cycle >= player['switchcyc']:
-		nsites += len(ghost_atoms) + len(lp_atoms)
-	
-	fh.write("{}\n".format(nsites))
-	fh.write("\nForce calculation - Cycle number {}\n".format(cycle))
-	
-	for atom in molecules[0]:
-		symbol = atomsymb[atom['na']]
-		fh.write("{:<2s}    {:>10.5f}   {:>10.5f}   {:>10.5f}\n".format(symbol, 
-													  atom['rx'], atom['ry'], atom['rz']))
-	
-	## If also performing charge fit in the same calculation
-	if cycle >= player['switchcyc']:
-		for ghost in ghost_atoms:
-			fh.write("Bq    {:>10.5f}   {:>10.5f}   {:>10.5f}\n".format(
-												   ghost['rx'], ghost['ry'], ghost['rz']))
-		
-		for lp in lp_atoms:
-			fh.write("Bq    {:>10.5f}   {:>10.5f}   {:>10.5f}\n".format(
-															lp['rx'], lp['ry'], lp['rz']))
-	
-	fh.write("basis = {}\n".format(molcas['basis']))
-	fh.write("group= nosym\n")
-	fh.write(" XFIELD\n")
-	fh.write(">>> Include asec.xfield\n")
-	
-	if not os.path.isfile(molcas['mbottom']):
-		sys.exit("Error: cannot find file {} in main directory".format(molcas['mbottom']))
-	try:
-		with open(molcas['mbottom']) as mbottomfile:
-			mbottom = mbottomfile.readlines()
-	except:
-		sys.exit("Error: cannot open file {}".format(molcas['mbottom']))
-	
-	for line in mbottom:
-		fh.write(line)
-	
-	fh.write(" &Alaska\nPNEW\n &SLAPAF\nCartesian\n")
-	fh.close()
-	
-	return
-
-
-
-def make_charge_input(cycle, asec_charges):
-	
-	path = "step{:02d}".format(cycle) + os.sep + "qm"
-	file = path + os.sep + "asec.input"
-	try:
-		fh = open(file, "w")
-	except:
-		sys.exit("Error: cannot open file {}".format(file))
-	
-	fh.write(" &Gateway\n")
-	fh.write("  Coord\n")
-	
-	nsites = len(molecules[0])
-	if cycle >= player['switchcyc']:
-		nsites += len(ghost_atoms) + len(lp_atoms)
-	
-	fh.write("{}\n".format(nsites))
-	fh.write("\nForce calculation - Cycle number {}\n".format(cycle))
-	
-	for atom in molecules[0]:
-		symbol = atomsymb[atom['na']]
-		fh.write("{:<2s}    {:>10.5f}   {:>10.5f}   {:>10.5f}\n".format(symbol, 
-													  atom['rx'], atom['ry'], atom['rz']))
-	
-	for ghost in ghost_atoms:
-		fh.write("Bq    {:>10.5f}   {:>10.5f}   {:>10.5f}\n".format(
-												   ghost['rx'], ghost['ry'], ghost['rz']))
-	
-	for lp in lp_atoms:
-		fh.write("Bq    {:>10.5f}   {:>10.5f}   {:>10.5f}\n".format(
-															lp['rx'], lp['ry'], lp['rz']))
-	
-	fh.write("basis = {}\n".format(molcas['basis']))
-	fh.write("group= nosym\n")
-	fh.write(" XFIELD\n")
-	fh.write(">>> Include asec.xfield\n")
-	
-	if not os.path.isfile(molcas['mbottom']):
-		sys.exit("Error: cannot find file {} in main directory".format(molcas['mbottom']))
-	try:
-		with open(molcas['mbottom']) as mbottomfile:
-			mbottom = mbottomfile.readlines()
-	except:
-		sys.exit("Error: cannot open file {}".format(molcas['mbottom']))
-	
-	for line in mbottom:
-		fh.write(line)
-	
-	fh.close()
-	
-	return
-
-
-
-def read_charges(file, fh):
-	
-	try:
-		with open(file) as tmpfh:
-			glogfile = tmpfh.readlines()
-	except:
-		sys.exit("Error: cannot open file {}".format(file))
-	
-	start = glogfile.pop(0).strip()
-	while  start != "Fitting point charges to electrostatic potential":
-		start = glogfile.pop(0).strip()
-	
-	glogfile = glogfile[3:]		## Consume 3 more lines
-	
-	fh.write("\nAtomic charges:\n")
-	fh.write("------------------------------------\n")
-	for atom in molecules[0]:
-		line = glogfile.pop(0).split()
-		atom_str = line[1]
-		charge = float(line[2])
-		atom['chg'] = charge
-		fh.write(" {:<2s}      {:>10.6f}\n".format(atom_str, charge))
-	
-	if gaussian['pop'] == "chelpg":
-		for ghost in ghost_atoms:
-			line = glogfile.pop(0).split()
-			atom_str = line[1]
-			charge = float(line[2])
-			ghost['chg'] = charge
-			fh.write(" {:<2s}      {:>10.6f}\n".format(atom_str, charge))
-		
-		for lp in lp_atoms:
-			line = glogfile.pop(0).split()
-			atom_str = line[1]
-			charge = float(line[2])
-			lp['chg'] = charge
-			fh.write(" {:<2s}      {:>10.6f}\n".format(atom_str, charge))
-	
-	fh.write("------------------------------------\n")
-		
-	return
-
-
-
-def run_gaussian(cycle, type, fh):
-	
-	path = "step{:02d}".format(cycle) + os.sep + "qm"
-	work_dir = os.getcwd()
-	os.chdir(path)
-	
-	if type == "force":
-		infile = "asec.gjf"
-	elif type == "charge":
-		infile = "asec2.gjf"
-	
-	fh.write("\nCalculation of {}s initiated with Gaussian on {}\n".format(type, date_time()))
-	
-	exit_status = subprocess.call([player['qmprog'], infile])
-	
-	if exit_status != 0:
-		sys.exit("Gaussian process did not exit properly")
-	
-	fh.write("Calculation of {}s finished on {}\n".format(type, date_time()))
-	
-	os.chdir(work_dir)
-	
-	return
-
-
-
-def run_formchk(cycle, fh):
-	
-	path = "step{:02d}".format(cycle) + os.sep + "qm"
-	work_dir = os.getcwd()
-	os.chdir(path)
-		
-	fh.write("Formatting the checkpoint file... ")
-	
-	exit_status = subprocess.call(["formchk", "asec.chk"])
-	
-	fh.write("Done\n")
-	
-	os.chdir(work_dir)
-	
-	return
-
-
-

diceplayer/DPpack/Optimization.py

@@ -1,266 +0,0 @@
-import sys, math
-from copy import deepcopy
-
-import numpy as np
-from numpy import linalg
-
-from DPpack.SetGlobals import *
-
-
-epsilon = 1e-8
-
-#######################################  functions  ######################################
-
-def best_previous_point():
-	
-	min_energy = 0
-	idx = 0
-	for energy in internal['energy'][:-1]:
-		if energy < min_energy or abs(energy - min_energy) < 1e-10:
-			min_energy = energy
-			min_idx = idx
-		idx += 1
-	
-	return min_idx
-
-
-
-def best_point():
-	
-	min_energy = 0
-	idx = 0
-	for energy in internal['energy']:
-		if energy < min_energy or abs(energy - min_energy) < 1e-10:
-			min_energy = energy
-			min_idx = idx
-		idx += 1
-	
-	return min_idx
-
-
-
-def line_search(fh):
-	
-	X1 = internal['position'][-1]		# numpy array
-	e1 = internal['energy'][-1]
-	G1 = internal['gradient'][-1]		# numpy array
-	
-	idx = best_previous_point()
-	X0 = internal['position'][idx]		# numpy array
-	e0 = internal['energy'][idx]
-	G0 = internal['gradient'][idx]		# numpy array
-	
-	# First try a quartic fit
-	fh.write("Attempting a quartic fit.\n")
-	success, y0 = quartic_fit(X0, X1, e0, e1, G0, G1, fh)
-	if success and y0 > 0:
-		if y0 < 1:
-			new_point = X0 + y0*(X1 - X0)
-			new_gradient = interpolate_gradient(G0, G1, y0)
-			new_gradient = perpendicular_projection(new_gradient, X1 - X0)
-			fh.write("Line search succeded.\n")
-			return True, new_point, new_gradient
-		else:
-			idx = best_point()
-			if idx == len(internal['energy']) - 1:
-				new_point = X0 + y0*(X1 - X0)
-				new_gradient = interpolate_gradient(G0, G1, y0)
-				new_gradient = perpendicular_projection(new_gradient, X1 - X0)
-				fh.write("Line search succeded.\n")
-				return True, new_point, new_gradient
-			else:
-				fh.write("Quartic step is not acceptable. ")
-	elif success:
-		fh.write("Quartic step is not acceptable. ")
-	
-	# If no condition is met, then y0 is unacceptable. Try the cubic fit next
-	fh.write("Attempting a cubic fit.\n")
-	success, y0 = cubic_fit(X0, X1, e0, e1, G0, G1, fh)
-	if success and y0 > 0:
-		if y0 < 1:
-			new_point = X0 + y0*(X1 - X0)
-			new_gradient = interpolate_gradient(G0, G1, y0)
-			new_gradient = perpendicular_projection(new_gradient, X1 - X0)
-			fh.write("Line search succeded.\n")
-			return True, new_point, new_gradient
-		else:
-			previous_step = X1 - internal['position'][-2]
-			previous_step_size = linalg.norm(previous_step)
-			new_point = X0 + y0*(X1 - X0)
-			step = new_point - X1
-			step_size = linalg.norm(step)
-			if step_size < previous_step_size:
-				new_gradient = interpolate_gradient(G0, G1, y0)
-				new_gradient = perpendicular_projection(new_gradient, X1 - X0)
-				fh.write("Line search succeded.\n")
-				return True, new_point, new_gradient
-			else:
-				fh.write("Cubic step is not acceptable. ")
-	elif success:
-		fh.write("Cubic step is not acceptable. ")
-	
-	# If no condition is met again, then all fits fail.
-	fh.write("All fits fail. ")
-	
-	# Then, if the latest point is not the best, use y0 = 0.5 (step to the midpoint)
-	idx = best_point()
-	if idx < len(internal['energy']) - 1:
-		y0 = 0.5
-		new_point = X0 + y0*(X1 - X0)
-		new_gradient = interpolate_gradient(G0, G1, y0)
-		new_gradient = perpendicular_projection(new_gradient, X1 - X0)
-		fh.write("Moving to the midpoint.\n")
-		return True, new_point, new_gradient
-	
-	# If the latest point is the best point, no linear search is done
-	fh.write("No linear search will be used in this step.\n")
-	
-	return False, None, None
-
-
-
-## For cubic and quartic fits, G0 and G1 are the gradient vectors
-
-def cubic_fit(X0, X1, e0, e1, G0, G1, fh):
-		
-	line = X1 - X0
-	line /= linalg.norm(line)
-	
-	g0 = np.dot(G0, line)
-	g1 = np.dot(G1, line)
-	
-	De = e1 - e0
-	
-	fh.write("De = {:<18.15e}      g0 = {:<12.8f}      g1 = {:<12.8f}\n".format(De, g0, g1))
-	
-	alpha = g1 + g0 - 2*De
-	if abs(alpha) < epsilon:
-		fh.write("Cubic fit failed: alpha too small\n")
-		return False, None
-	
-	beta = 3*De - 2*g0 - g1
-	discriminant = 4 * (beta**2 - 3*alpha*g0)
-	if discriminant < 0:
-		fh.write("Cubic fit failed: no minimum found (negative Delta)\n")
-		return False, None
-	if abs(discriminant) < epsilon:
-		fh.write("Cubic fit failed: no minimum found (null Delta)\n")
-		return False, None
-	
-	y0 = (-beta + math.sqrt(discriminant/4)) / (3*alpha)
-	fh.write("Minimum found with y0 = {:<8.4f}\n".format(y0))
-	
-	return True, y0
-
-
-
-def quartic_fit(X0, X1, e0, e1, G0, G1, fh):
-		
-	line = X1 - X0
-	line /= linalg.norm(line)
-	
-	g0 = np.dot(G0, line)
-	g1 = np.dot(G1, line)
-	
-	De = e1 - e0
-	Dg = g1 - g0
-	
-	fh.write("De = {:<18.15e}      g0 = {:<12.8f}      g1 = {:<12.8f}\n".format(De, g0, g1))
-	
-	if Dg < 0 or De - g0 < 0:
-		fh.write("Quartic fit failed: negative alpha\n")
-		return False, None
-	if abs(Dg) < epsilon or abs(De - g0) < epsilon:
-		fh.write("Quartic fit failed: alpha too small\n")
-		return False, None
-	
-	discriminant = 16 * (Dg**2 - 3*(g1 + g0 - 2*De)**2)
-	if discriminant < 0:
-		fh.write("Quartic fit failed: no minimum found (negative Delta)\n")
-		return False, None
-	
-	alpha1 = (Dg + math.sqrt(discriminant/16)) / 2
-	alpha2 = (Dg - math.sqrt(discriminant/16)) / 2
-	
-	fh.write("alpha1 = {:<7.4e}      alpha2 = {:<7.4e}\n".format(alpha1, alpha2))
-	
-	alpha = alpha1
-	beta = g1 + g0 - 2*De - 2*alpha
-	gamma = De - g0 - alpha - beta
-	
-	y0 = (-1/(2*alpha)) * ((beta**3 - 4*alpha*beta*gamma + 8*g0*alpha**2)/4)**(1/3)
-	fh.write("Minimum found with y0 = {:<8.4f}\n".format(y0))
-	
-	return True, y0
-
-
-
-def rfo_step(gradient, hessian, type):
-	
-	dim = len(gradient)
-	
-	aug_hessian = []
-	for i in range(dim):
-		aug_hessian.extend(hessian[i,:].tolist())
-		aug_hessian.append(gradient[i])
-	
-	aug_hessian.extend(gradient.tolist())
-	aug_hessian.append(0)
-	
-	aug_hessian = np.array(aug_hessian).reshape(dim + 1, dim + 1)
-	
-	evals, evecs = linalg.eigh(aug_hessian)
-	
-	if type == "min":
-		step = np.array(evecs[:-1,0])
-	elif type == "ts":
-		step = np.array(evecs[:-1,1])
-	
-	return step
-
-
-
-def update_trust_radius():
-	
-	if internal['trust_radius'] == None:
-		internal['trust_radius'] = player['maxstep']
-	elif len(internal['energy']) > 1:
-		X1 = internal['position'][-1]
-		X0 = internal['position'][-2]
-		Dx = X1 - X0
-		displace = linalg.norm(Dx)
-		e1 = internal['energy'][-1]
-		e0 = internal['energy'][-2]
-		De = e1 - e0
-		g0 = internal['gradient'][-2]
-		h0 = internal['hessian'][-2]
-		
-		rho = De / (np.dot(g0, Dx) + 0.5*np.dot(Dx, np.matmul(h0, Dx.T).T))
-		
-		if rho > 0.75 and displace > 0.8*internal['trust_radius']:
-			internal['trust_radius'] = 2*internal['trust_radius']
-		elif rho < 0.25:
-			internal['trust_radius'] = 0.25*displace
-	
-	return
-
-
-
-def interpolate_gradient(G0, G1, y0):
-	
-	DG = G1 - G0
-	gradient = G0 + y0*DG
-	
-	return gradient
-
-
-
-def perpendicular_projection(vector, line):
-	
-	direction = line / linalg.norm(line)
-	projection = np.dot(vector, direction) * direction
-	
-	return vector - projection
-
-
-

diceplayer/DPpack/Player.py

@@ -0,0 +1,830 @@
+from diceplayer.DPpack.Utils.Validations import NotNull
+from diceplayer.DPpack.Utils.PTable import *
+from diceplayer.DPpack.Utils.Misc import *
+
+from diceplayer.DPpack.External.Gaussian import Gaussian
+from diceplayer.DPpack.External.Dice import Dice
+
+from diceplayer.DPpack.Environment.System import System
+from diceplayer.DPpack.Environment.Molecule import Molecule
+from diceplayer.DPpack.Environment.Atom import Atom
+
+from typing import TextIO
+
+import os
+import sys
+import shutil
+import textwrap
+import types
+import yaml
+
+
+env = ["OMP_STACKSIZE"]
+
+
+class Player:
+    def __init__(self, infile: TextIO, outfile: TextIO) -> None:
+
+        self.infile = infile
+        self.outfile = outfile
+
+        self.continued: bool = None
+
+        self.system = System()
+
+        self.player = self.Player()
+        self.player_keywords = [
+            a
+            for a in dir(self.player)
+            if not a.startswith("__") and not callable(getattr(self.player, a))
+        ]
+
+        self.dice = Dice(infile, outfile)
+        self.dice_keywords = [
+            a
+            for a in dir(self.dice)
+            if not a.startswith("__") and not callable(getattr(self.dice, a))
+        ]
+
+        self.gaussian = Gaussian()
+        self.gaussian_keywords = [
+            a
+            for a in dir(self.gaussian)
+            if not a.startswith("__") and not callable(getattr(self.gaussian, a))
+        ]
+
+        self.TOL_RMS_FORCE = 3e-4
+        self.TOL_MAX_FORCE = 4.5e-4
+        self.TOL_RMS_STEP = 1.2e-3
+        self.TOL_MAX_SET = 1.8e-3
+        self.TRUST_RADIUS = None
+
+        self.combrule = None
+
+
+    @NotNull(requiredArgs=["maxcyc", "opt", "nprocs", "qmprog", "lps", "ghosts", "altsteps"])
+    def updateKeywords(self, **data):
+        self.__dict__.update(data)
+
+    def read_keywords(self) -> None:
+
+        with self.infile as f:
+            data = yaml.load(f, Loader=yaml.SafeLoader)
+
+        self.updateKeywords(data.get("diceplayer"))
+        self.dice.updateKeywords(data.get("dice"))
+        self.gaussian.updateKeywords(data.get("gaussian"))
+
+    def check_keywords(self) -> None:
+
+        min_steps = 20000
+
+        if self.dice.ljname == None:
+            sys.exit(
+                "Error: 'ljname' keyword not specified in file {}".format(self.infile)
+            )
+
+        if self.dice.outname == None:
+            sys.exit(
+                "Error: 'outname' keyword not specified in file {}".format(self.infile)
+            )
+
+        if self.dice.dens == None:
+            sys.exit(
+                "Error: 'dens' keyword not specified in file {}".format(self.infile)
+            )
+
+        if self.dice.nmol == 0:
+            sys.exit(
+                "Error: 'nmol' keyword not defined appropriately in file {}".format(
+                    self.infile
+                )
+            )
+
+        if self.dice.nstep == 0:
+            sys.exit(
+                "Error: 'nstep' keyword not defined appropriately in file {}".format(
+                    self.infile
+                )
+            )
+
+        # Check only if QM program is Gaussian:
+        if self.player.qmprog in ("g03", "g09", "g16"):
+
+            if self.gaussian.level == None:
+                sys.exit(
+                    "Error: 'level' keyword not specified in file {}".format(
+                        self.infile
+                    )
+                )
+
+            if self.gaussian.gmiddle != None:
+                if not os.path.isfile(self.gaussian.gmiddle):
+                    sys.exit("Error: file {} not found".format(self.gaussian.gmiddle))
+
+            if self.gaussian.gbottom != None:
+                if not os.path.isfile(self.gaussian.gbottom):
+                    sys.exit("Error: file {} not found".format(self.gaussian.gbottom))
+
+            if self.gaussian.pop != "chelpg" and (
+                self.player.ghosts == "yes" or self.player.lps == "yes"
+            ):
+                sys.exit(
+                    "Error: ghost atoms or lone pairs only available with 'pop = chelpg')"
+                )
+
+        # Check only if QM program is Molcas:
+        # if self.player.qmprog == "molcas":
+
+        # 	if self.molcas.mbottom == None:
+        # 		sys.exit("Error: 'mbottom' keyword not specified in file {}".format(self.infile))
+        # 	else:
+        # 		if not os.path.isfile(self.molcas.mbottom):
+        # 			sys.exit("Error: file {} not found".format(self.molcas.mbottom))
+
+        # 	if self.molcas.basis == None:
+        # 		sys.exit("Error: 'basis' keyword not specified in file {}".format(self.infile))
+
+        if self.player.altsteps != 0:
+
+            # Verifica se tem mais de 1 molecula QM
+            # (No futuro usar o RMSD fit para poder substituir todas as moleculas QM
+            # no arquivo outname.xy - Need to change the __make_init_file!!)
+            if self.dice.nmol[0] > 1:
+                sys.exit(
+                    "Error: altsteps > 0 only possible with 1 QM molecule (nmol = 1 n2 n3 n4)"
+                )
+
+            # if not zero, altsteps cannot be less than min_steps
+            self.player.altsteps = max(min_steps, self.player.altsteps)
+            # altsteps value is always the nearest multiple of 1000
+            self.player.altsteps = round(self.player.altsteps / 1000) * 1000
+
+        for i in range(len(self.dice.nstep)):
+            # nstep can never be less than min_steps
+            self.dice.nstep[i] = max(min_steps, self.dice.nstep[i])
+            # nstep values are always the nearest multiple of 1000
+            self.dice.nstep[i] = round(self.dice.nstep[i] / 1000) * 1000
+
+        # isave must be between 100 and 2000
+        self.dice.isave = max(100, self.dice.isave)
+        self.dice.isave = min(2000, self.dice.isave)
+        # isave value is always the nearest multiple of 100
+        self.dice.isave = round(self.dice.isave / 100) * 100
+
+    def print_keywords(self) -> None:
+
+        self.outfile.write(
+            "##########################################################################################\n"
+            "#############               Welcome to DICEPLAYER version 1.0                #############\n"
+            "##########################################################################################\n"
+            "\n"
+        )
+        self.outfile.write("Your python version is {}\n".format(sys.version))
+        self.outfile.write("\n")
+        self.outfile.write("Program started on {}\n".format(weekday_date_time()))
+        self.outfile.write("\n")
+        self.outfile.write("Environment variables:\n")
+        for var in env:
+            self.outfile.write(
+                "{} = {}\n".format(
+                    var, (os.environ[var] if var in os.environ else "Not set")
+                )
+            )
+
+        self.outfile.write(
+            "\n==========================================================================================\n"
+            "                         CONTROL variables being used in this run:\n"
+            "------------------------------------------------------------------------------------------\n"
+            "\n"
+        )
+
+        for key in sorted(self.player_keywords):
+            if getattr(self.player, key) != None:
+                if isinstance(getattr(self.player, key), list):
+                    string = " ".join(str(x) for x in getattr(self.player, key))
+                    self.outfile.write("{} = {}\n".format(key, string))
+                else:
+                    self.outfile.write(
+                        "{} = {}\n".format(key, getattr(self.player, key))
+                    )
+
+        self.outfile.write("\n")
+
+        self.outfile.write(
+            "------------------------------------------------------------------------------------------\n"
+            "                         DICE variables being used in this run:\n"
+            "------------------------------------------------------------------------------------------\n"
+            "\n"
+        )
+
+        for key in sorted(self.dice_keywords):
+            if getattr(self.dice, key) != None:
+                if isinstance(getattr(self.dice, key), list):
+                    string = " ".join(str(x) for x in getattr(self.dice, key))
+                    self.outfile.write("{} = {}\n".format(key, string))
+                else:
+                    self.outfile.write("{} = {}\n".format(key, getattr(self.dice, key)))
+
+        self.outfile.write("\n")
+
+        if self.player.qmprog in ("g03", "g09", "g16"):
+
+            self.outfile.write(
+                "------------------------------------------------------------------------------------------\n"
+                "                         GAUSSIAN variables being used in this run:\n"
+                "------------------------------------------------------------------------------------------\n"
+                "\n"
+            )
+
+            for key in sorted(self.gaussian_keywords):
+                if getattr(self.gaussian, key) != None:
+                    if isinstance(getattr(self.gaussian, key), list):
+                        string = " ".join(str(x) for x in getattr(self.gaussian, key))
+                        self.outfile.write("{} = {}\n".format(key, string))
+                    else:
+                        self.outfile.write(
+                            "{} = {}\n".format(key, getattr(self.gaussian, key))
+                        )
+
+            self.outfile.write("\n")
+
+        # elif self.player.qmprog == "molcas":
+
+        # 	self.outfile.write("------------------------------------------------------------------------------------------\n"
+        # 			"                         MOLCAS variables being used in this run:\n"
+        # 			"------------------------------------------------------------------------------------------\n"
+        # 			"\n")
+
+        # 	for key in sorted(molcas):
+        # 		if molcas[key] != None:
+        # 			if isinstance(molcas[key], list):
+        # 				string = " ".join(str(x) for x in molcas[key])
+        # 				self.outfile.write("{} = {}\n".format(key, string))
+        # 			else:
+        # 				self.outfile.write("{} = {}\n".format(key, molcas[key]))
+
+        # 	self.outfile.write("\n")
+
+    def read_potential(self) -> None:  # Deve ser atualizado para o uso de
+
+        try:
+            with open(self.dice.ljname) as file:
+                ljfile = file.readlines()
+        except EnvironmentError as err:
+            sys.exit(err)
+
+        combrule = ljfile.pop(0).split()[0]
+        if combrule not in ("*", "+"):
+            sys.exit(
+                "Error: expected a '*' or a '+' sign in 1st line of file {}".format(
+                    self.dice.ljname
+                )
+            )
+        self.dice.combrule = combrule
+
+        ntypes = ljfile.pop(0).split()[0]
+        if not ntypes.isdigit():
+            sys.exit(
+                "Error: expected an integer in the 2nd line of file {}".format(
+                    self.dice.ljname
+                )
+            )
+        ntypes = int(ntypes)
+
+        if ntypes != len(self.dice.nmol):
+            sys.exit(
+                "Error: number of molecule types in file {} must match that of 'nmol' keyword in file {}".format(
+                    self.dice.ljname, self.infile
+                )
+            )
+        line = 2
+        for i in range(ntypes):
+
+            line += 1
+            nsites, molname = ljfile.pop(0).split()[:2]
+
+            if not nsites.isdigit():
+                sys.exit(
+                    "Error: expected an integer in line {} of file {}".format(
+                        line, self.dice.ljname
+                    )
+                )
+
+            if molname is None:
+                sys.exit(
+                    "Error: expected a molecule name in line {} of file {}".format(
+                        line, self.dice.ljname
+                    )
+                )
+
+            nsites = int(nsites)
+
+            self.system.add_type(nsites, Molecule(molname))
+
+            for j in range(nsites):
+
+                line += 1
+                new_atom = ljfile.pop(0).split()
+
+                if len(new_atom) < 8:
+                    sys.exit(
+                        "Error: expected at least 8 fields in line {} of file {}".format(
+                            line, self.dice.ljname
+                        )
+                    )
+
+                if not new_atom[0].isdigit():
+                    sys.exit(
+                        "Error: expected an integer in field 1, line {} of file {}".format(
+                            line, self.dice.ljname
+                        )
+                    )
+                lbl = int(new_atom[0])
+
+                if not new_atom[1].isdigit():
+                    sys.exit(
+                        "Error: expected an integer in field 2, line {} of file {}".format(
+                            line, self.dice.ljname
+                        )
+                    )
+
+                atnumber = int(new_atom[1])
+                if (
+                    atnumber == ghost_number and i == 0
+                ):  # Ghost atom not allowed in the QM molecule
+                    sys.exit(
+                        "Error: found a ghost atom in line {} of file {}".format(
+                            line, self.dice.ljname
+                        )
+                    )
+                na = atnumber
+
+                try:
+                    rx = float(new_atom[2])
+                except:
+                    sys.exit(
+                        "Error: expected a float in field 3, line {} of file {}".format(
+                            line, self.dice.ljname
+                        )
+                    )
+
+                try:
+                    ry = float(new_atom[3])
+                except:
+                    sys.exit(
+                        "Error: expected a float in field 4, line {} of file {}".format(
+                            line, self.dice.ljname
+                        )
+                    )
+
+                try:
+                    rz = float(new_atom[4])
+                except:
+                    sys.exit(
+                        "Error: expected a float in field 5, line {} of file {}".format(
+                            line, self.dice.ljname
+                        )
+                    )
+
+                try:
+                    chg = float(new_atom[5])
+                except:
+                    sys.exit(
+                        "Error: expected a float in field 6, line {} of file {}".format(
+                            line, self.dice.ljname
+                        )
+                    )
+
+                try:
+                    eps = float(new_atom[6])
+                except:
+                    sys.exit(
+                        "Error: expected a float in field 7, line {} of file {}".format(
+                            line, self.dice.ljname
+                        )
+                    )
+
+                try:
+                    sig = float(new_atom[7])
+                except:
+                    sys.exit(
+                        "Error: expected a float in field 8, line {} of file {}".format(
+                            line, self.dice.ljname
+                        )
+                    )
+
+                mass = atommass[na]
+
+                if len(new_atom) > 8:
+                    masskey, mass = new_atom[8].partition("=")[::2]
+                    if masskey.lower() == "mass" and len(mass) != 0:
+                        try:
+                            new_mass = float(mass)
+                            if new_mass > 0:
+                                mass = new_mass
+                        except:
+                            sys.exit(
+                                "Error: expected a positive float after 'mass=' in field 9, line {} of file {}".format(
+                                    line, self.dice.ljname
+                                )
+                            )
+
+                self.system.molecule[i].add_atom(
+                    Atom(lbl, na, rx, ry, rz, chg, eps, sig)
+                )
+
+        to_delete = ["lbl", "na", "rx", "ry", "rz", "chg", "eps", "sig", "mass"]
+        for _var in to_delete:
+            if _var in locals() or _var in globals():
+                exec(f"del {_var}")
+
+    def print_potential(self) -> None:
+
+        formatstr = "{:<3d} {:>3d}  {:>10.5f} {:>10.5f} {:>10.5f}  {:>10.6f} {:>9.5f} {:>7.4f} {:>9.4f}\n"
+        self.outfile.write(
+            "\n"
+            "==========================================================================================\n"
+        )
+        self.outfile.write(
+            "                    Potential parameters from file {}:\n".format(
+                self.dice.ljname
+            )
+        )
+        self.outfile.write(
+            "------------------------------------------------------------------------------------------\n"
+            "\n"
+        )
+
+        self.outfile.write("Combination rule: {}\n".format(self.dice.combrule))
+        self.outfile.write(
+            "Types of molecules: {}\n\n".format(len(self.system.molecule))
+        )
+
+        i = 0
+        for mol in self.system.molecule:
+            i += 1
+            self.outfile.write(
+                "{} atoms in molecule type {}:\n".format(len(mol.atom), i)
+            )
+            self.outfile.write(
+                "---------------------------------------------------------------------------------\n"
+                "Lbl  AN       X          Y          Z         Charge    Epsilon   Sigma     Mass\n"
+            )
+            self.outfile.write(
+                "---------------------------------------------------------------------------------\n"
+            )
+
+            for atom in mol.atom:
+
+                self.outfile.write(
+                    formatstr.format(
+                        atom.lbl,
+                        atom.na,
+                        atom.rx,
+                        atom.ry,
+                        atom.rz,
+                        atom.chg,
+                        atom.eps,
+                        atom.sig,
+                        atom.mass,
+                    )
+                )
+
+            self.outfile.write("\n")
+
+        if self.player.ghosts == "yes" or self.player.lps == "yes":
+            self.outfile.write(
+                "\n"
+                "------------------------------------------------------------------------------------------\n"
+                "                    Aditional potential parameters:\n"
+                "------------------------------------------------------------------------------------------\n"
+            )
+
+        # if player['ghosts'] == "yes":
+
+        # 	self.outfile.write("\n")
+        # 	self.outfile.write("{} ghost atoms appended to molecule type 1 at:\n".format(len(ghost_types)))
+        # 	self.outfile.write("---------------------------------------------------------------------------------\n")
+
+        # 	atoms_string = ""
+        # 	for ghost in ghost_types:
+        # 		for atom in ghost['numbers']:
+        # 			atom_sym = atomsymb[ molecules[0][atom - 1]['na'] ].strip()
+        # 			atoms_string += "{}{} ".format(atom_sym,atom)
+
+        # 		if ghost['type'] == "g":
+        # 			self.outfile.write(textwrap.fill("* Geometric center of atoms {}".format(atoms_string), 80))
+        # 		elif ghost['type'] == "m":
+        # 			self.outfile.write(textwrap.fill("* Center of mass of atoms {}".format(atoms_string), 80))
+        # 		elif ghost['type'] == "z":
+        # 			self.outfile.write(textwrap.fill("* Center of atomic number of atoms {}".format(atoms_string), 80))
+
+        # 		self.outfile.write("\n")
+
+        # if player['lps'] == 'yes':
+
+        # 	self.outfile.write("\n")
+        # 	self.outfile.write("{} lone pairs appended to molecule type 1:\n".format(len(lp_types)))
+        # 	self.outfile.write("---------------------------------------------------------------------------------\n")
+
+        # 	for lp in lp_types:
+        # 		# LP type 1 or 2
+        # 		if lp['type'] in (1, 2):
+        # 			atom1_num = lp['numbers'][0]
+        # 			atom1_sym = atomsymb[ molecules[0][atom1_num - 1]['na'] ].strip()
+        # 			atom2_num = lp['numbers'][1]
+        # 			atom2_sym = atomsymb[ molecules[0][atom2_num - 1]['na'] ].strip()
+        # 			atom3_num = lp['numbers'][2]
+        # 			atom3_sym = atomsymb[ molecules[0][atom3_num - 1]['na'] ].strip()
+
+        # 			self.outfile.write(textwrap.fill(
+        # 			"* Type {} on atom {}{} with {}{} {}{}. Alpha = {:<5.1f} Deg and D = {:<4.2f} Angs".format(
+        # 				lp['type'], atom1_sym, atom1_num, atom2_sym, atom2_num, atom3_sym, atom3_num, lp['alpha'],
+        # 				lp['dist']), 86))
+        # 			self.outfile.write("\n")
+
+        # 		# Other LP types
+
+        self.outfile.write(
+            "\n"
+            "==========================================================================================\n"
+        )
+
+    def check_executables(self) -> None:
+
+        self.outfile.write("\n")
+        self.outfile.write(90 * "=")
+        self.outfile.write("\n\n")
+
+        dice_path = shutil.which(self.dice.progname)
+        if dice_path != None:
+            self.outfile.write(
+                "Program {} found at {}\n".format(self.dice.progname, dice_path)
+            )
+            self.dice.path = dice_path
+        else:
+            sys.exit("Error: cannot find dice executable")
+
+        qmprog_path = shutil.which(self.gaussian.qmprog)
+        if qmprog_path != None:
+            self.outfile.write(
+                "Program {} found at {}\n".format(self.gaussian.qmprog, qmprog_path)
+            )
+            self.gaussian.path = qmprog_path
+        else:
+            sys.exit("Error: cannot find {} executable".format(self.gaussian.qmprog))
+
+        if self.gaussian.qmprog in ("g03", "g09", "g16"):
+            formchk_path = shutil.which("formchk")
+            if formchk_path != None:
+                self.outfile.write("Program formchk found at {}\n".format(formchk_path))
+            else:
+                sys.exit("Error: cannot find formchk executable")
+
+    def dice_start(self, cycle: int):
+
+        self.dice.configure(
+            self.player.initcyc,
+            self.player.nprocs,
+            self.player.altsteps,
+            self.system.nmols,
+            self.system.molecule,
+        )
+
+        self.dice.start(cycle)
+
+        self.dice.reset()
+
+    def gaussian_start(self, cycle: int, geomsfh: TextIO):
+
+        self.gaussian.configure(
+            self.player.initcyc,
+            self.player.nprocs,
+            self.dice.ncores,
+            self.player.altsteps,
+            self.player.switchcyc,
+            self.player.opt,
+            self.system.nmols,
+            self.system.molecule,
+        )
+
+        position = self.gaussian.start(cycle, self.outfile, self.player.readhessian)
+
+        ## Update the geometry of the reference molecule
+        self.system.update_molecule(position, self.outfile)
+
+        ## Print new geometry in geoms.xyz
+        self.system.print_geom(cycle, geomsfh)
+
+        self.gaussian.reset()
+
+    # I still have to talk with Herbet about this function
+    def populate_asec_vdw(self, cycle):
+
+        # Both asec_charges and vdw_meanfield will utilize the Molecule() class and Atoms() with some None elements
+
+        asec_charges = Molecule(
+            "ASEC_CHARGES"
+        )  # (lbl=None, na=None, rx, ry, rz, chg, eps=None, sig=None)
+        # vdw_meanfield = (
+        #     Molecule()
+        # )  # (lbl=None, na=None, rx, ry, rz, chg=None, eps, sig)
+
+        if self.dice.nstep[-1] % self.dice.isave == 0:
+            nconfigs = round(self.dice.nstep[-1] / self.dice.isave)
+        else:
+            nconfigs = int(self.dice.nstep[-1] / self.dice.isave)
+
+        norm_factor = nconfigs * self.player.nprocs
+
+        nsitesref = len(self.system.molecule[0].atom)
+        # nsitesref = (
+        #     len(self.system.molecule[0].atom)
+        #     + len(self.system.molecule[0].ghost_atoms)
+        #     + len(self.system.molecule[0].lp_atoms)
+        # )
+
+        nsites_total = self.dice.nmol[0] * nsitesref
+        for i in range(1, len(self.dice.nmol)):
+            nsites_total += self.dice.nmol[i] * len(self.system.molecule[i].atom)
+
+        thickness = []
+        picked_mols = []
+
+        for proc in range(1, self.player.nprocs + 1):  # Run over folders
+
+            simdir = "simfiles"
+            path = (
+                simdir
+                + os.sep
+                + "step{:02d}".format(cycle)
+                + os.sep
+                + "p{:02d}".format(proc)
+            )
+            file = path + os.sep + self.dice.outname + ".xyz"
+            if not os.path.isfile(file):
+                sys.exit("Error: cannot find file {}".format(file))
+            try:
+                with open(file) as xyzfh:
+                    xyzfile = xyzfh.readlines()
+            except:
+                sys.exit("Error: cannot open file {}".format(file))
+
+            for config in range(nconfigs):  # Run over configs in a folder
+
+                if int(xyzfile.pop(0).split()[0]) != nsites_total:
+                    sys.exit("Error: wrong number of sites in file {}".format(file))
+
+                box = xyzfile.pop(0).split()[-3:]
+                box = [float(box[0]), float(box[1]), float(box[2])]
+                sizes = self.system.molecule[0].sizes_of_molecule()
+                thickness.append(
+                    min(
+                        [
+                            (box[0] - sizes[0]) / 2,
+                            (box[1] - sizes[1]) / 2,
+                            (box[2] - sizes[2]) / 2,
+                        ]
+                    )
+                )
+
+                # Skip the first (reference) molecule
+                xyzfile = xyzfile[nsitesref:]
+                mol_count = 0
+                for type in range(len(self.dice.nmol)):  # Run over types of molecules
+
+                    if type == 0:
+                        nmols = self.dice.nmol[0] - 1
+                    else:
+                        nmols = self.dice.nmol[type]
+
+                    for mol in range(nmols):  # Run over molecules of each type
+
+                        new_molecule = Molecule(self.system.molecule[type].molname)
+                        # Run over sites of each molecule
+                        for site in range(len(self.system.molecule[types].atom)):
+
+                            # new_molecule.append({})
+                            line = xyzfile.pop(0).split()
+
+                            if (
+                                line[0].title()
+                                != atomsymb[
+                                    self.system.molecule[type].atom[site].na.strip()
+                                ]
+                            ):
+                                sys.exit("Error reading file {}".format(file))
+
+                            new_molecule.add_atom(
+                                Atom(
+                                    self.system.molecule[type].atom[site].lbl,
+                                    self.system.molecule[type].atom[site].na,
+                                    self.system.molecule[type]
+                                    .atom[site]
+                                    .float(line[1]),
+                                    self.system.molecule[type]
+                                    .atom[site]
+                                    .float(line[2]),
+                                    self.system.molecule[type]
+                                    .atom[site]
+                                    .float(line[3]),
+                                    self.system.molecule[type].atom[site].chg,
+                                    self.system.molecule[type].atom[site].eps,
+                                    self.system.molecule[type].atom[site].sig,
+                                )
+                            )
+
+                        dist = self.system.molecule[0].minimum_distance(new_molecule)
+                        if dist < thickness[-1]:
+                            mol_count += 1
+                            for atom in new_molecule.atom:
+                                asec_charges.append({})
+                                # vdw_meanfield.append({})
+
+                                asec_charges[-1]["rx"] = atom.rx
+                                asec_charges[-1]["ry"] = atom.ry
+                                asec_charges[-1]["rz"] = atom.rz
+                                asec_charges[-1]["chg"] = atom.chg / norm_factor
+
+                                # if self.player.vdwforces == "yes":
+                                #     vdw_meanfield[-1]["rx"] = atom["rx"]
+                                #     vdw_meanfield[-1]["ry"] = atom["ry"]
+                                #     vdw_meanfield[-1]["rz"] = atom["rz"]
+                                #     vdw_meanfield[-1]["eps"] = atom["eps"]
+                                #     vdw_meanfield[-1]["sig"] = atom["sig"]
+
+                        # ####  Read lines with ghosts or lps in molecules of type 0 (reference)
+                        # ####  and, if dist < thickness, appends to asec
+                        # if type == 0:
+                        # 	for ghost in ghost_atoms:
+                        # 		line = xyzfile.pop(0).split()
+                        # 		if line[0] != dice_ghost_label:
+                        # 			sys.exit("Error reading file {}".format(file))
+                        # 		if dist < thickness[-1]:
+                        # 			asec_charges.append({})
+                        # 			asec_charges[-1]['rx'] = float(line[1])
+                        # 			asec_charges[-1]['ry'] = float(line[2])
+                        # 			asec_charges[-1]['rz'] = float(line[3])
+                        # 			asec_charges[-1]['chg'] = ghost['chg'] / norm_factor
+
+                        # 	for lp in lp_atoms:
+                        # 		line = xyzfile.pop(0).split()
+                        # 		if line[0] != dice_ghost_label:
+                        # 			sys.exit("Error reading file {}".format(file))
+                        # 		if dist < thickness[-1]:
+                        # 			asec_charges.append({})
+                        # 			asec_charges[-1]['rx'] = float(line[1])
+                        # 			asec_charges[-1]['ry'] = float(line[2])
+                        # 			asec_charges[-1]['rz'] = float(line[3])
+                        # 			asec_charges[-1]['chg'] = lp['chg'] / norm_factor
+
+                picked_mols.append(mol_count)
+
+        self.outfile.write("Done\n")
+
+        string = "In average, {:^7.2f} molecules ".format(
+            sum(picked_mols) / norm_factor
+        )
+        string += "were selected from each of the {} configurations ".format(
+            len(picked_mols)
+        )
+        string += (
+            "of the production simulations to form the ASEC, comprising a shell with "
+        )
+        string += "minimum thickness of {:>6.2f} Angstrom\n".format(
+            sum(thickness) / norm_factor
+        )
+
+        self.outfile.write(textwrap.fill(string, 86))
+        self.outfile.write("\n")
+
+        otherfh = open("ASEC.dat", "w")
+        for charge in asec_charges:
+            otherfh.write(
+                "{:>10.5f}   {:>10.5f}   {:>10.5f}     {:>11.8f}\n".format(
+                    charge["rx"], charge["ry"], charge["rz"], charge["chg"]
+                )
+            )
+        otherfh.close()
+
+        return asec_charges
+
+    class Player:
+        def __init__(self) -> None:
+
+            self.maxcyc = None
+            self.nprocs = 1
+            self.switchcyc = 3
+            self.altsteps = 20000
+            self.maxstep = 0.3
+            self.opt = "yes"
+            self.freq = "no"
+            self.readhessian = "no"
+            self.lps = "no"
+            self.ghosts = "no"
+            self.vdwforces = "no"
+            self.tol_factor = 1.2
+            self.qmprog = "g16"
+
+            self.initcyc = 1

diceplayer/DPpack/Utils/Optimization.py

@@ -0,0 +1,263 @@
+# import sys, math
+# from copy import deepcopy
+
+# import numpy as np
+# from numpy import linalg
+
+# from diceplayer.DPpack.SetGlobals import *
+
+
+# epsilon = 1e-8
+
+# #######################################  functions  ######################################
+
+
+# def best_previous_point():
+
+#     min_energy = 0
+#     idx = 0
+#     for energy in internal["energy"][:-1]:
+#         if energy < min_energy or abs(energy - min_energy) < 1e-10:
+#             min_energy = energy
+#             min_idx = idx
+#         idx += 1
+
+#     return min_idx
+
+
+# def best_point():
+
+#     min_energy = 0
+#     idx = 0
+#     for energy in internal["energy"]:
+#         if energy < min_energy or abs(energy - min_energy) < 1e-10:
+#             min_energy = energy
+#             min_idx = idx
+#         idx += 1
+
+#     return min_idx
+
+
+# def line_search(fh):
+
+#     X1 = internal["position"][-1]  # numpy array
+#     e1 = internal["energy"][-1]
+#     G1 = internal["gradient"][-1]  # numpy array
+
+#     idx = best_previous_point()
+#     X0 = internal["position"][idx]  # numpy array
+#     e0 = internal["energy"][idx]
+#     G0 = internal["gradient"][idx]  # numpy array
+
+#     # First try a quartic fit
+#     fh.write("Attempting a quartic fit.\n")
+#     success, y0 = quartic_fit(X0, X1, e0, e1, G0, G1, fh)
+#     if success and y0 > 0:
+#         if y0 < 1:
+#             new_point = X0 + y0 * (X1 - X0)
+#             new_gradient = interpolate_gradient(G0, G1, y0)
+#             new_gradient = perpendicular_projection(new_gradient, X1 - X0)
+#             fh.write("Line search succeded.\n")
+#             return True, new_point, new_gradient
+#         else:
+#             idx = best_point()
+#             if idx == len(internal["energy"]) - 1:
+#                 new_point = X0 + y0 * (X1 - X0)
+#                 new_gradient = interpolate_gradient(G0, G1, y0)
+#                 new_gradient = perpendicular_projection(new_gradient, X1 - X0)
+#                 fh.write("Line search succeded.\n")
+#                 return True, new_point, new_gradient
+#             else:
+#                 fh.write("Quartic step is not acceptable. ")
+#     elif success:
+#         fh.write("Quartic step is not acceptable. ")
+
+#     # If no condition is met, then y0 is unacceptable. Try the cubic fit next
+#     fh.write("Attempting a cubic fit.\n")
+#     success, y0 = cubic_fit(X0, X1, e0, e1, G0, G1, fh)
+#     if success and y0 > 0:
+#         if y0 < 1:
+#             new_point = X0 + y0 * (X1 - X0)
+#             new_gradient = interpolate_gradient(G0, G1, y0)
+#             new_gradient = perpendicular_projection(new_gradient, X1 - X0)
+#             fh.write("Line search succeded.\n")
+#             return True, new_point, new_gradient
+#         else:
+#             previous_step = X1 - internal["position"][-2]
+#             previous_step_size = linalg.norm(previous_step)
+#             new_point = X0 + y0 * (X1 - X0)
+#             step = new_point - X1
+#             step_size = linalg.norm(step)
+#             if step_size < previous_step_size:
+#                 new_gradient = interpolate_gradient(G0, G1, y0)
+#                 new_gradient = perpendicular_projection(new_gradient, X1 - X0)
+#                 fh.write("Line search succeded.\n")
+#                 return True, new_point, new_gradient
+#             else:
+#                 fh.write("Cubic step is not acceptable. ")
+#     elif success:
+#         fh.write("Cubic step is not acceptable. ")
+
+#     # If no condition is met again, then all fits fail.
+#     fh.write("All fits fail. ")
+
+#     # Then, if the latest point is not the best, use y0 = 0.5 (step to the midpoint)
+#     idx = best_point()
+#     if idx < len(internal["energy"]) - 1:
+#         y0 = 0.5
+#         new_point = X0 + y0 * (X1 - X0)
+#         new_gradient = interpolate_gradient(G0, G1, y0)
+#         new_gradient = perpendicular_projection(new_gradient, X1 - X0)
+#         fh.write("Moving to the midpoint.\n")
+#         return True, new_point, new_gradient
+
+#     # If the latest point is the best point, no linear search is done
+#     fh.write("No linear search will be used in this step.\n")
+
+#     return False, None, None
+
+
+# ## For cubic and quartic fits, G0 and G1 are the gradient vectors
+
+
+# def cubic_fit(X0, X1, e0, e1, G0, G1, fh):
+
+#     line = X1 - X0
+#     line /= linalg.norm(line)
+
+#     g0 = np.dot(G0, line)
+#     g1 = np.dot(G1, line)
+
+#     De = e1 - e0
+
+#     fh.write(
+#         "De = {:<18.15e}      g0 = {:<12.8f}      g1 = {:<12.8f}\n".format(De, g0, g1)
+#     )
+
+#     alpha = g1 + g0 - 2 * De
+#     if abs(alpha) < epsilon:
+#         fh.write("Cubic fit failed: alpha too small\n")
+#         return False, None
+
+#     beta = 3 * De - 2 * g0 - g1
+#     discriminant = 4 * (beta**2 - 3 * alpha * g0)
+#     if discriminant < 0:
+#         fh.write("Cubic fit failed: no minimum found (negative Delta)\n")
+#         return False, None
+#     if abs(discriminant) < epsilon:
+#         fh.write("Cubic fit failed: no minimum found (null Delta)\n")
+#         return False, None
+
+#     y0 = (-beta + math.sqrt(discriminant / 4)) / (3 * alpha)
+#     fh.write("Minimum found with y0 = {:<8.4f}\n".format(y0))
+
+#     return True, y0
+
+
+# def quartic_fit(X0, X1, e0, e1, G0, G1, fh):
+
+#     line = X1 - X0
+#     line /= linalg.norm(line)
+
+#     g0 = np.dot(G0, line)
+#     g1 = np.dot(G1, line)
+
+#     De = e1 - e0
+#     Dg = g1 - g0
+
+#     fh.write(
+#         "De = {:<18.15e}      g0 = {:<12.8f}      g1 = {:<12.8f}\n".format(De, g0, g1)
+#     )
+
+#     if Dg < 0 or De - g0 < 0:
+#         fh.write("Quartic fit failed: negative alpha\n")
+#         return False, None
+#     if abs(Dg) < epsilon or abs(De - g0) < epsilon:
+#         fh.write("Quartic fit failed: alpha too small\n")
+#         return False, None
+
+#     discriminant = 16 * (Dg**2 - 3 * (g1 + g0 - 2 * De) ** 2)
+#     if discriminant < 0:
+#         fh.write("Quartic fit failed: no minimum found (negative Delta)\n")
+#         return False, None
+
+#     alpha1 = (Dg + math.sqrt(discriminant / 16)) / 2
+#     alpha2 = (Dg - math.sqrt(discriminant / 16)) / 2
+
+#     fh.write("alpha1 = {:<7.4e}      alpha2 = {:<7.4e}\n".format(alpha1, alpha2))
+
+#     alpha = alpha1
+#     beta = g1 + g0 - 2 * De - 2 * alpha
+#     gamma = De - g0 - alpha - beta
+
+#     y0 = (-1 / (2 * alpha)) * (
+#         (beta**3 - 4 * alpha * beta * gamma + 8 * g0 * alpha**2) / 4
+#     ) ** (1 / 3)
+#     fh.write("Minimum found with y0 = {:<8.4f}\n".format(y0))
+
+#     return True, y0
+
+
+# def rfo_step(gradient, hessian, type):
+
+#     dim = len(gradient)
+
+#     aug_hessian = []
+#     for i in range(dim):
+#         aug_hessian.extend(hessian[i, :].tolist())
+#         aug_hessian.append(gradient[i])
+
+#     aug_hessian.extend(gradient.tolist())
+#     aug_hessian.append(0)
+
+#     aug_hessian = np.array(aug_hessian).reshape(dim + 1, dim + 1)
+
+#     evals, evecs = linalg.eigh(aug_hessian)
+
+#     if type == "min":
+#         step = np.array(evecs[:-1, 0])
+#     elif type == "ts":
+#         step = np.array(evecs[:-1, 1])
+
+#     return step
+
+
+# def update_trust_radius():
+
+#     if internal["trust_radius"] == None:
+#         internal["trust_radius"] = player["maxstep"]
+#     elif len(internal["energy"]) > 1:
+#         X1 = internal["position"][-1]
+#         X0 = internal["position"][-2]
+#         Dx = X1 - X0
+#         displace = linalg.norm(Dx)
+#         e1 = internal["energy"][-1]
+#         e0 = internal["energy"][-2]
+#         De = e1 - e0
+#         g0 = internal["gradient"][-2]
+#         h0 = internal["hessian"][-2]
+
+#         rho = De / (np.dot(g0, Dx) + 0.5 * np.dot(Dx, np.matmul(h0, Dx.T).T))
+
+#         if rho > 0.75 and displace > 0.8 * internal["trust_radius"]:
+#             internal["trust_radius"] = 2 * internal["trust_radius"]
+#         elif rho < 0.25:
+#             internal["trust_radius"] = 0.25 * displace
+
+#     return
+
+
+# def interpolate_gradient(G0, G1, y0):
+
+#     DG = G1 - G0
+#     gradient = G0 + y0 * DG
+
+#     return gradient
+
+
+# def perpendicular_projection(vector, line):
+
+#     direction = line / linalg.norm(line)
+#     projection = np.dot(vector, direction) * direction
+
+#     return vector - projection

diceplayer/DPpack/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

diceplayer/DPpack/external/Dice.py

@@ -1,662 +0,0 @@
-from diceplayer.DPpack.MolHandling import *
-from diceplayer.DPpack.PTable import *
-from diceplayer.DPpack.Misc import *
-
-
-from typing import IO, Tuple, List, TextIO, Union
-
-from numpy.core.numeric import partition
-from numpy import random
-
-from multiprocessing import Process, connection
-import setproctitle
-import subprocess
-import os
-import sys
-import shutil
-import textwrap
-import types
-
-
-dice_end_flag = "End of simulation"  # The normal end flag
-dice_flag_line = -2  # must be in the line before the last
-umaAng3_to_gcm3 = 1.6605  # Conversion between uma/Ang3 to g/cm3
-
-max_seed = 4294967295  # Maximum allowed value for a seed (numpy)
-
-
-class Dice:
-
-	def __init__(self, infile: TextIO, outfile: TextIO) -> None:
-
-		self.title = "Diceplayer run"
-		self.progname = "dice"
-
-		self.infile = infile
-		self.outfile = outfile
-
-		self.path = None
-		self.nprocs: int = None
-
-		self.randominit = 'first'
-		self.temp = 300.0
-		self.press = 1.0
-		self.isave = 1000         # ASEC construction will take this into account
-		self.ncores = 1
-
-		# Investigate the possibility of using 'box = Lx Ly Lz' instead.
-		self.dens = None
-		# self.box = None		# So 'geom' would be set by diceplayer and 'cutoff' would be
-		# switched off. One of them must be given.
-		self.combrule = "*"
-		self.ljname = None
-		self.outname = None
-		# Up to 4 integer values related to up to 4 molecule types
-		self.nmol: List[int] = []
-		# 2 or 3 integer values related to 2 or 3 simulations
-		self.nstep: List[int] = []
-		# (NVT th + NVT eq) or (NVT th + NPT th + NPT eq).
-		# This will control the 'nstep' keyword of Dice
-		self.upbuf = 360
-
-	def __new_density(self, cycle: int, proc: int) -> float:
-
-		sim_dir = "simfiles"
-		step_dir = "step{:02d}".format(cycle-1)
-		proc_dir = "p{:02d}".format(proc)
-		path = sim_dir + os.sep + step_dir + os.sep + proc_dir
-		file = path + os.sep + "last.xyz"
-
-		if not os.path.isfile(file):
-			sys.exit(
-				"Error: cannot find the xyz file {} in main directory".format(file))
-		try:
-			with open(file) as fh:
-				xyzfile = fh.readlines()
-		except:
-			sys.exit("Error: cannot open file {}".format(file))
-
-		box = xyzfile[1].split()
-		volume = float(box[-3]) * float(box[-2]) * float(box[-1])
-
-		total_mass = 0
-		for i in range(len(self.molecule)):
-
-			total_mass += self.molecule[i].total_mass * self.nmol[i]
-
-		density = (total_mass / volume) * umaAng3_to_gcm3
-
-		return density
-
-	def __print_last_config(self, cycle: int, proc: int) -> None:
-		
-		sim_dir = "simfiles"
-		step_dir = "step{:02d}".format(cycle)
-		proc_dir = "p{:02d}".format(proc)
-		path = sim_dir + os.sep + step_dir + os.sep + proc_dir
-		file = path + os.sep + "phb.xyz"
-		if not os.path.isfile(file):
-			sys.exit("Error: cannot find the xyz file {}".format(file))
-		try:
-			with open(file) as fh:
-				xyzfile = fh.readlines()
-		except:
-			sys.exit("Error: cannot open file {}".format(file))
-		
-		nsites = len(self.molecule[0].atom) * self.nmol[0]	
-		for i in range(1, len(self.nmol)):
-			nsites += self.nmol[i] * len(self.molecule[i].atom)
-		
-		nsites += 2  ## To include the comment line and the number of atoms (xyz file format)
-		
-		nsites *= -1  ## Become an index to count from the end of xyzfile (list)
-		xyzfile = xyzfile[nsites :]  ## Take the last configuration
-		
-		
-		file = path + os.sep + "last.xyz"
-		fh = open(file, "w")
-		for line in xyzfile:
-			fh.write(line)
-
-	def __make_dice_inputs(self, cycle: int, proc: int) -> None:
-
-		sim_dir = "simfiles"
-		step_dir = "step{:02d}".format(cycle)
-		proc_dir = "p{:02d}".format(proc)
-		path = sim_dir + os.sep + step_dir + os.sep + proc_dir
-
-		num = time.time()  # Take the decimal places 7 to 12 of the
-		num = (num - int(num)) * 1e6  # time in seconds as a floating point
-		# to make an integer in the range 1-1e6
-		num = int((num - int(num)) * 1e6)
-		random.seed((os.getpid() * num) % (max_seed + 1))
-
-		if (self.randominit == 'first' and cycle > self.initcyc):
-			last_step_dir = "step{:02d}".format(cycle-1)
-			last_path = sim_dir + os.sep + last_step_dir + os.sep + proc_dir
-			xyzfile = last_path + os.sep + "last.xyz"
-			self.__make_init_file(path, xyzfile)
-
-		if len(self.nstep) == 2:  # Means NVT simulation
-
-			self.__make_nvt_ter(cycle, path)
-			self.__make_nvt_eq(path)
-
-		elif len(self.nstep) == 3:  # Means NPT simulation
-
-			if (self.randominit == 'first' and cycle > self.initcyc):
-				self.dens = self.__new_density(cycle, proc)
-			else:
-				self.__make_nvt_ter(cycle, path)
-
-			self.__make_npt_ter(cycle, path)
-			self.__make_npt_eq(path)
-
-		else:
-			sys.exit("Error: bad number of entries for 'nstep'")
-
-		self.__make_potential(path)
-
-		# if (self.randominit == 'first' and cycle > self.initcyc):
-
-		# 	last_path = sim_dir + os.sep + "step{:02d}".format(cycle-1) + os.sep + proc_dir
-		# 	shutil.copyfile(last_path + os.sep + "phb.dat", path + os.sep + "phb.dat")
-
-	def __make_nvt_ter(self, cycle: int, path: str) -> None:
-
-		file = path + os.sep + "NVT.ter"
-		try:
-			fh = open(file, "w")
-		except:
-			sys.exit("Error: cannot open file {}".format(file))
-
-		fh.write("title = {} - NVT Thermalization\n".format(self.title))
-		fh.write("ncores = {}\n".format(self.ncores))
-		fh.write("ljname = {}\n".format(self.ljname))
-		fh.write("outname = {}\n".format(self.outname))
-
-		string = " ".join(str(x) for x in self.nmol)
-		fh.write("nmol = {}\n".format(string))
-
-		fh.write("dens = {}\n".format(self.dens))
-		fh.write("temp = {}\n".format(self.temp))
-
-		if (self.randominit == 'first' and cycle > self.initcyc):
-			fh.write("init = yesreadxyz\n")
-			fh.write("nstep = {}\n".format(self.altsteps))
-		else:
-			fh.write("init = yes\n")
-			fh.write("nstep = {}\n".format(self.nstep[0]))
-
-		fh.write("vstep = 0\n")
-		fh.write("mstop = 1\n")
-		fh.write("accum = no\n")
-		fh.write("iprint = 1\n")
-		fh.write("isave = 0\n")
-		fh.write("irdf = 0\n")
-
-		seed = int(1e6 * random.random())
-		fh.write("seed = {}\n".format(seed))
-		fh.write("upbuf = {}".format(self.upbuf))
-
-		fh.close()
-
-	def __make_nvt_eq(self, path: str) -> None:
-
-		file = path + os.sep + "NVT.eq"
-		try:
-			fh = open(file, "w")
-		except:
-			sys.exit("Error: cannot open file {}".format(file))
-
-		fh.write("title = {} - NVT Production\n".format(self.title))
-		fh.write("ncores = {}\n".format(self.ncores))
-		fh.write("ljname = {}\n".format(self.ljname))
-		fh.write("outname = {}\n".format(self.outname))
-
-		string = " ".join(str(x) for x in self.nmol)
-		fh.write("nmol = {}\n".format(string))
-
-		fh.write("dens = {}\n".format(self.dens))
-		fh.write("temp = {}\n".format(self.temp))
-		fh.write("init = no\n")
-		fh.write("nstep = {}\n".format(self.nstep[1]))
-		fh.write("vstep = 0\n")
-		fh.write("mstop = 1\n")
-		fh.write("accum = no\n")
-		fh.write("iprint = 1\n")
-		fh.write("isave = {}\n".format(self.isave))
-		fh.write("irdf = {}\n".format(10 * self.nprocs))
-
-		seed = int(1e6 * random.random())
-		fh.write("seed = {}\n".format(seed))
-
-		fh.close()
-
-	def __make_npt_ter(self, cycle: int, path: str) -> None:
-
-		file = path + os.sep + "NPT.ter"
-		try:
-			fh = open(file, "w")
-		except:
-			sys.exit("Error: cannot open file {}".format(file))
-
-		fh.write("title = {} - NPT Thermalization\n".format(self.title))
-		fh.write("ncores = {}\n".format(self.ncores))
-		fh.write("ljname = {}\n".format(self.ljname))
-		fh.write("outname = {}\n".format(self.outname))
-
-		string = " ".join(str(x) for x in self.nmol)
-		fh.write("nmol = {}\n".format(string))
-
-		fh.write("press = {}\n".format(self.press))
-		fh.write("temp = {}\n".format(self.temp))
-
-		if (self.randominit == 'first' and cycle > self.initcyc):
-			fh.write("init = yesreadxyz\n")
-			fh.write("dens = {:<8.4f}\n".format(self.dens))
-			fh.write("vstep = {}\n".format(int(self.altsteps / 5)))
-		else:
-			# Because there will be a previous NVT simulation
-			fh.write("init = no\n")
-			fh.write("vstep = {}\n".format(int(self.nstep[1] / 5)))
-
-		fh.write("nstep = 5\n")
-		fh.write("mstop = 1\n")
-		fh.write("accum = no\n")
-		fh.write("iprint = 1\n")
-		fh.write("isave = 0\n")
-		fh.write("irdf = 0\n")
-
-		seed = int(1e6 * random.random())
-		fh.write("seed = {}\n".format(seed))
-
-		fh.close()
-
-	def __make_npt_eq(self, path: str) -> None:
-
-		file = path + os.sep + "NPT.eq"
-		try:
-			fh = open(file, "w")
-		except:
-			sys.exit("Error: cannot open file {}".format(file))
-
-		fh.write("title = {} - NPT Production\n".format(self.title))
-		fh.write("ncores = {}\n".format(self.ncores))
-		fh.write("ljname = {}\n".format(self.ljname))
-		fh.write("outname = {}\n".format(self.outname))
-
-		string = " ".join(str(x) for x in self.nmol)
-		fh.write("nmol = {}\n".format(string))
-
-		fh.write("press = {}\n".format(self.press))
-		fh.write("temp = {}\n".format(self.temp))
-
-		fh.write("nstep = 5\n")
-
-		fh.write("vstep = {}\n".format(int(self.nstep[2] / 5)))
-		fh.write("init = no\n")
-		fh.write("mstop = 1\n")
-		fh.write("accum = no\n")
-		fh.write("iprint = 1\n")
-		fh.write("isave = {}\n".format(self.isave))
-		fh.write("irdf = {}\n".format(10 * self.nprocs))
-
-		seed = int(1e6 * random.random())
-		fh.write("seed = {}\n".format(seed))
-
-		fh.close()
-
-	def __make_init_file(self, path: str, file: TextIO) -> None:
-
-		if not os.path.isfile(file):
-			sys.exit(
-				"Error: cannot find the xyz file {} in main directory".format(file))
-		try:
-			with open(file) as fh:
-				xyzfile = fh.readlines()
-		except:
-			sys.exit("Error: cannot open file {}".format(file))
-
-		nsites_mm = 0
-		for i in range(1, len(self.nmol)):
-			nsites_mm += self.nmol[i] * len(self.molecule[i].atom)
-
-		# Become an index to count from the end of xyzfile (list)
-		nsites_mm *= -1
-		# Only the MM atoms of the last configuration remains
-		xyzfile = xyzfile[nsites_mm:]
-
-		file = path + os.sep + self.outname + ".xy"
-
-		try:
-			fh = open(file, "w", 1)
-		except:
-			sys.exit("Error: cannot open file {}".format(file))
-
-		for atom in self.molecule[0].atom:
-			fh.write("{:>10.6f}  {:>10.6f}  {:>10.6f}\n".format(
-				atom.rx, atom.ry, atom.rz))
-
-		# for i in self.molecule[0].ghost_atoms:
-		# 	with self.molecule[0].atom[i] as ghost:
-		# 		fh.write("{:>10.6f}  {:>10.6f}  {:>10.6f}\n".format(ghost.rx, ghost.ry, ghost.rz))
-
-		# for i in self.molecule[0].lp_atoms:
-		# 	with self.molecule[0].atom[i] as lp:
-		# 		fh.write("{:>10.6f}  {:>10.6f}  {:>10.6f}\n".format(lp.rx, lp.ry, lp.rz))
-
-		for line in xyzfile:
-			atom = line.split()
-			rx = float(atom[1])
-			ry = float(atom[2])
-			rz = float(atom[3])
-			fh.write("{:>10.6f}  {:>10.6f}  {:>10.6f}\n".format(rx, ry, rz))
-
-		fh.write("$end")
-
-		fh.close()
-
-	def __make_potential(self, path: str) -> None:
-
-		fstr = "{:<3d} {:>3d}  {:>10.5f} {:>10.5f} {:>10.5f}  {:>10.6f} {:>9.5f} {:>7.4f}\n"
-
-		file = path + os.sep + self.ljname
-		try:
-			fh = open(file, "w")
-		except:
-			sys.exit("Error: cannot open file {}".format(file))
-
-		fh.write("{}\n".format(self.combrule))
-		fh.write("{}\n".format(len(self.nmol)))
-
-		nsites_qm = len(self.molecule[0].atom) + len(
-			self.molecule[0].ghost_atoms) + len(self.molecule[0].lp_atoms)
-
-		# Print the sites of the QM self.molecule
-		fh.write("{} {}\n".format(nsites_qm, self.molecule[0].molname))
-		for atom in self.molecule[0].atom:
-			fh.write(fstr.format(atom.lbl, atom.na, atom.rx, atom.ry, atom.rz,
-								 atom.chg, atom.eps, atom.sig))
-
-		ghost_label = self.molecule[0].atom[-1].lbl + 1
-		for i in self.molecule[0].ghost_atoms:
-			fh.write(fstr.format(ghost_label, ghost_number, self.molecule[0].atom[i].rx, self.molecule[0].atom[i].ry,
-								 self.molecule[0].atom[i].rz, self.molecule[0].atom[i].chg, 0, 0))
-
-		ghost_label += 1
-		for lp in self.molecule[0].lp_atoms:
-			fh.write(fstr.format(ghost_label, ghost_number, lp['rx'], lp['ry'], lp['rz'],
-								 lp['chg'], 0, 0))
-
-		# Print the sites of the other self.molecules
-		for mol in self.molecule[1:]:
-			fh.write("{} {}\n".format(len(mol.atom), mol.molname))
-			for atom in mol.atom:
-				fh.write(fstr.format(atom.lbl, atom.na, atom.rx, atom.ry,
-									 atom.rz, atom.chg, atom.eps, atom.sig))
-
-	def __make_proc_dir(self, cycle: int, proc: int) -> None:
-
-		sim_dir = "simfiles"
-		step_dir = "step{:02d}".format(cycle)
-		proc_dir = "p{:02d}".format(proc)
-		path = sim_dir + os.sep + step_dir + os.sep + proc_dir
-		try:
-			os.makedirs(path)
-		except:
-			sys.exit("Error: cannot make directory {}".format(path))
-
-	def __run_dice(self, cycle: int, proc: int, fh: str) -> None:
-
-		sim_dir = "simfiles"
-		step_dir = "step{:02d}".format(cycle)
-		proc_dir = "p{:02d}".format(proc)
-
-		try:
-			fh.write("Simulation process {} initiated with pid {}\n".format(
-				sim_dir + os.sep + step_dir + os.sep + proc_dir, os.getpid()))
-
-		except Exception as err:
-			print("I/O error({0}): {1}".format(err))
-
-		path = sim_dir + os.sep + step_dir + os.sep + proc_dir
-		working_dir = os.getcwd()
-		os.chdir(path)
-
-		if len(self.nstep) == 2:  # Means NVT simulation
-
-			if (self.randominit == 'first' and cycle > self.initcyc):
-				string_tmp = 'previous'
-			else:
-				string_tmp = 'random'
-
-			# NVT thermalization
-			string = "(from " + string_tmp + " configuration)"
-			fh.write("p{:02d}> NVT thermalization finished {} on {}\n".format(proc, string,
-																			  date_time()))
-
-			infh = open("NVT.ter")
-			outfh = open("NVT.ter.out", "w")
-
-			if shutil.which("bash") != None:
-				exit_status = subprocess.call(
-					["bash", "-c", "exec -a dice-step{}-p{} {} < {} > {}".format(cycle, proc, self.progname, infh.name, outfh.name)])
-			else:
-				exit_status = subprocess.call(
-					self.progname, stin=infh.name, stout=outfh.name)
-
-			infh.close()
-			outfh.close()
-
-			if os.getppid() == 1:  # Parent process is dead
-				sys.exit()
-
-			if exit_status != 0:
-				sys.exit(
-					"Dice process step{:02d}-p{:02d} did not exit properly".format(cycle, proc))
-			else:
-				# Open again to seek the normal end flag
-				outfh = open("NVT.ter.out")
-				flag = outfh.readlines()[dice_flag_line].strip()
-				outfh.close()
-				if flag != dice_end_flag:
-					sys.exit(
-						"Dice process step{:02d}-p{:02d} did not exit properly".format(cycle, proc))
-
-			# NVT production
-			fh.write("p{:02d}> NVT production initiated on {}\n".format(
-				proc, date_time()))
-
-			infh = open("NVT.eq")
-			outfh = open("NVT.eq.out", "w")
-
-			if shutil.which("bash") != None:
-				exit_status = subprocess.call(
-					["bash", "-c", "exec -a dice-step{}-p{} {} < {} > {}".format(cycle, proc, self.progname, infh.name, outfh.name)])
-			else:
-				exit_status = subprocess.call(
-					self.progname, stin=infh.name, stout=outfh.name)
-
-			infh.close()
-			outfh.close()
-
-			if os.getppid() == 1:  # Parent process is dead
-				sys.exit()
-
-			if exit_status != 0:
-				sys.exit(
-					"Dice process step{:02d}-p{:02d} did not exit properly".format(cycle, proc))
-			else:
-				# Open again to seek the normal end flag
-				outfh = open("NVT.eq.out")
-				flag = outfh.readlines()[dice_flag_line].strip()
-				outfh.close()
-				if flag != dice_end_flag:
-					sys.exit(
-						"Dice process step{:02d}-p{:02d} did not exit properly".format(cycle, proc))
-
-			fh.write("p{:02d}> ----- NVT production finished on {}\n".format(proc,
-																			 date_time()))
-
-		elif len(self.nstep) == 3:  # Means NPT simulation
-
-			# NVT thermalization if randominit
-			if self.randominit == 'always' or (self.randominit == 'first' and cycle == 1) or self.continued:
-				string = "(from random configuration)"
-				fh.write("p{:02d}> NVT thermalization initiated {} on {}\n".format(proc,
-																				   string, date_time()))
-				infh = open("NVT.ter")
-				outfh = open("NVT.ter.out", "w")
-
-				if shutil.which("bash") != None:
-					exit_status = subprocess.call(
-						["bash", "-c", "exec -a dice-step{}-p{} {} < {} > {}".format(cycle, proc, self.progname, infh.name, outfh.name)])
-				else:
-					exit_status = subprocess.call(
-						self.progname, stin=infh.name, stout=outfh.name)
-
-				infh.close()
-				outfh.close()
-
-				if os.getppid() == 1:  # Parent process is dead
-					sys.exit()
-
-				if exit_status != 0:
-					sys.exit(
-						"Dice process step{:02d}-p{:02d} did not exit properly".format(cycle, proc))
-				else:
-					# Open again to seek the normal end flag
-					outfh = open("NVT.ter.out")
-					flag = outfh.readlines()[dice_flag_line].strip()
-					outfh.close()
-					if flag != dice_end_flag:
-						sys.exit(
-							"Dice process step{:02d}-p{:02d} did not exit properly".format(cycle, proc))
-
-			# NPT thermalization
-			if not self.randominit == 'always' or ((self.randominit == 'first' and cycle > self.initcyc)):
-				string = " (from previous configuration) "
-			else:
-				string = " "
-			fh.write("p{:02d}> NPT thermalization finished {} on {}\n".format(proc, string,
-																			  date_time()))
-
-			infh = open("NPT.ter")
-			outfh = open("NPT.ter.out", "w")
-
-			if shutil.which("bash") != None:
-				exit_status = subprocess.call(
-					["bash", "-c", "exec -a dice-step{}-p{} {} < {} > {}".format(cycle, proc, self.progname, infh.name, outfh.name)])
-			else:
-				exit_status = subprocess.call(
-					self.progname, stin=infh.name, stout=outfh.name)
-
-			infh.close()
-			outfh.close()
-
-			if os.getppid() == 1:  # Parent process is dead
-				sys.exit()
-
-			if exit_status != 0:
-				sys.exit(
-					"Dice process step{:02d}-p{:02d} did not exit properly".format(cycle, proc))
-			else:
-				# Open again to seek the normal end flag
-				outfh = open("NPT.ter.out")
-				flag = outfh.readlines()[dice_flag_line].strip()
-				outfh.close()
-				if flag != dice_end_flag:
-					sys.exit(
-						"Dice process step{:02d}-p{:02d} did not exit properly".format(cycle, proc))
-
-			# NPT production
-			fh.write("p{:02d}> NPT production initiated on {}\n".format(
-				proc, date_time()))
-
-			infh = open("NPT.eq")
-			outfh = open("NPT.eq.out", "w")
-
-			if shutil.which("bash") != None:
-				exit_status = subprocess.call(
-					["bash", "-c", "exec -a dice-step{}-p{} {} < {} > {}".format(cycle, proc, self.progname, infh.name, outfh.name)])
-			else:
-				exit_status = subprocess.call(
-					self.progname, stin=infh.name, stout=outfh.name)
-
-			infh.close()
-			outfh.close()
-
-			if os.getppid() == 1:  # Parent process is dead
-				sys.exit()
-
-			if exit_status != 0:
-				sys.exit(
-					"Dice process step{:02d}-p{:02d} did not exit properly".format(cycle, proc))
-			else:
-				# Open again to seek the normal end flag
-				outfh = open("NPT.eq.out")
-				flag = outfh.readlines()[dice_flag_line].strip()
-				outfh.close()
-				if flag != dice_end_flag:
-					sys.exit(
-						"Dice process step{:02d}-p{:02d} did not exit properly".format(cycle, proc))
-
-			fh.write("p{:02d}> ----- NPT production finished on {}\n".format(proc,
-																			 date_time()))
-
-		os.chdir(working_dir)
-
-	def __simulation_process(self, cycle: int, proc: int):
-		setproctitle.setproctitle(
-			"diceplayer-step{:0d}-p{:0d}".format(cycle, proc))
-
-		try:
-			self.__make_proc_dir(cycle, proc)
-			self.__make_dice_inputs(cycle, proc)
-			self.__run_dice(cycle, proc, self.outfile)
-		except Exception as err:
-			sys.exit(err)
-
-	def configure(self, initcyc: int, nprocs: int, altsteps: int, nmol: List[int], molecule: List[Molecule]):
-
-		self.initcyc = initcyc
-
-		self.nprocs = nprocs
-		self.altsteps = altsteps
-
-		self.molecule = molecule
-
-	def start(self, cycle: int) -> None:
-
-		procs = []
-		sentinels = []
-
-		for proc in range(1, self.nprocs + 1):
-
-			p = Process(target=self.__simulation_process, args=(cycle, proc))
-			p.start()
-
-			procs.append(p)
-			sentinels.append(p.sentinel)
-
-		while procs:
-			finished = connection.wait(sentinels)
-			for proc_sentinel in finished:
-				i = sentinels.index(proc_sentinel)
-				status = procs[i].exitcode
-				procs.pop(i)
-				sentinels.pop(i)
-				if status != 0:
-					for p in procs:
-						p.terminate()
-					sys.exit(status)
-
-		for proc in range(1, self.nprocs + 1):
-			self.__print_last_config(cycle, proc)
-
-	def reset(self):
-
-		del self.nprocs
-		del self.altsteps
-		del self.molecule

diceplayer/__main__.py

@@ -8,48 +8,64 @@ import argparse
 import shutil
 import pickle
 
-from diceplayer.DPpack.PTable import *
-from diceplayer.DPpack.SetGlobals import *
-from diceplayer.DPpack.MolHandling import *
-from diceplayer.DPpack.Misc import *
+from diceplayer.DPpack.Utils.Misc import *
 
-__version = 'dev'
+from diceplayer.DPpack.Player import Player
+
+from diceplayer.DPpack.Environment.Molecule import Molecule
+from diceplayer.DPpack.Environment.Atom import Atom
+
+__version = "dev"
 setproctitle.setproctitle("diceplayer-{}".format(__version))
 
-if __name__ == '__main__':
-####  Read and store the arguments passed to the program  ####
-####  and set the usage and help messages                 ####
+if __name__ == "__main__":
+    ####  Read and store the arguments passed to the program  ####
+    ####  and set the usage and help messages                 ####
 
-	parser = argparse.ArgumentParser(prog='Diceplayer')
-	parser.add_argument('--continue', dest='opt_continue' , default=False, action='store_true')
-	parser.add_argument('--version', action='version', version= "diceplayer-"+__version)
-	parser.add_argument('-i', dest='infile', default='control.in', metavar='INFILE', 
-								   help='input file of diceplayer [default = control.in]')
-	parser.add_argument('-o', dest='outfile', default='run.log', metavar='OUTFILE', 
-									 help='output file of diceplayer [default = run.log]')
+    parser = argparse.ArgumentParser(prog="Diceplayer")
+    parser.add_argument(
+        "--continue", dest="opt_continue", default=False, action="store_true"
+    )
+    parser.add_argument(
+        "--version", action="version", version="diceplayer-" + __version
+    )
+    parser.add_argument(
+        "-i",
+        dest="infile",
+        default="control.in",
+        metavar="INFILE",
+        help="input file of diceplayer [default = control.in]",
+    )
+    parser.add_argument(
+        "-o",
+        dest="outfile",
+        default="run.log",
+        metavar="OUTFILE",
+        help="output file of diceplayer [default = run.log]",
+    )
     ## Study the option of a parameter for continuing the last process via data from control.in and run.log files
 
     args = parser.parse_args()
 
-####  Open OUTFILE for writing and print keywords and initial info
+    ####  Open OUTFILE for writing and print keywords and initial info
 
     try:
 
         if args.opt_continue and os.path.exists(args.outfile):
 
-			save = pickle.load(open("latest-step.pkl","rb"))
+            save = pickle.load(open("latest-step.pkl", "rb"))
 
-			if os.path.isfile(args.outfile+".backup"):
-				os.remove(args.outfile+".backup")
+            if os.path.isfile(args.outfile + ".backup"):
+                os.remove(args.outfile + ".backup")
 
-			os.rename(args.outfile,args.outfile+".backup")
-			outfile = open(args.outfile,'w',1)
+            os.rename(args.outfile, args.outfile + ".backup")
+            outfile = open(args.outfile, "w", 1)
 
         elif os.path.exists(args.outfile):
-			os.rename(args.outfile, args.outfile+".backup")
-			outfile = open(args.outfile,'w',1)
+            os.rename(args.outfile, args.outfile + ".backup")
+            outfile = open(args.outfile, "w", 1)
         else:
-			outfile = open(args.outfile,"w",1)
+            outfile = open(args.outfile, "w", 1)
 
     except Exception as err:
         sys.exit(err)
@@ -57,65 +73,71 @@ if __name__ == '__main__':
     try:
 
         if os.path.exists(args.infile):
-			infile = open(args.infile,"r")
+            infile = open(args.infile, "r")
 
     except Exception as err:
         sys.exit(err)
 
-####  Read and check the keywords in INFILE
+    ####  Read and check the keywords in INFILE
 
-	internal = Internal(infile, outfile)
+    player = Player(infile, outfile)
 
-	internal.read_keywords()
+    player.read_keywords()
 
-	internal.check_keywords()
-	internal.print_keywords()
+    player.check_keywords()
+    player.print_keywords()
 
     if args.opt_continue:
-		internal.player.initcyc = save[0] + 1
-		internal.system = save[1]
+        player.player.initcyc = save[0] + 1
+        player.system = save[1]
     else:
-		internal.read_potential()
+        player.read_potential()
 
-####  Check whether the executables are in the path
-####		and print potential to Log File
+    ####  Check whether the executables are in the path
+    ####		and print potential to Log File
 
-	internal.check_executables()
+    player.check_executables()
 
-	internal.print_potential()
+    player.print_potential()
 
-####  Bring the molecules to standard orientation and prints info about them
+    ####  Bring the molecules to standard orientation and prints info about them
 
-	for i in range(len(internal.system.molecule)):
+    for i in range(len(player.system.molecule)):
 
-		internal.outfile.write("\nMolecule type {} - {}:\n\n".format(i + 1, internal.system.molecule[i].molname))
-		internal.system.molecule[i].print_mol_info(internal.outfile)
-		internal.outfile.write("    Translating and rotating molecule to standard orientation...")
-		internal.system.molecule[i].standard_orientation()
-		internal.outfile.write(" Done\n\n    New values:\n")
-		internal.system.molecule[i].print_mol_info(internal.outfile)
+        player.outfile.write(
+            "\nMolecule type {} - {}:\n\n".format(
+                i + 1, player.system.molecule[i].molname
+            )
+        )
+        player.system.molecule[i].print_mol_info(player.outfile)
+        player.outfile.write(
+            "    Translating and rotating molecule to standard orientation..."
+        )
+        player.system.molecule[i].standard_orientation()
+        player.outfile.write(" Done\n\n    New values:\n")
+        player.system.molecule[i].print_mol_info(player.outfile)
 
-	internal.outfile.write(90 * "=")
-	internal.outfile.write("\n")
+    player.outfile.write(90 * "=")
+    player.outfile.write("\n")
 
     if not args.opt_continue:
         make_simulation_dir()
     else:
         simdir = "simfiles"
-		stepdir = "step{:02d}".format(internal.player.initcyc)
-		if os.path.exists(simdir+os.sep+stepdir):
-			shutil.rmtree(simdir+os.sep+stepdir)
+        stepdir = "step{:02d}".format(player.player.initcyc)
+        if os.path.exists(simdir + os.sep + stepdir):
+            shutil.rmtree(simdir + os.sep + stepdir)
 
-####  Open the geoms.xyz file and prints the initial geometry if starting from zero
+    ####  Open the geoms.xyz file and prints the initial geometry if starting from zero
 
-	if internal.player.initcyc == 1:
+    if player.player.initcyc == 1:
         try:
             path = "geoms.xyz"
             geomsfh = open(path, "w", 1)
         except EnvironmentError as err:
             sys.exit(err)
-		internal.system.print_geom(0, geomsfh)
-		geomsfh.write(40 * "-"+"\n")
+        player.system.print_geom(0, geomsfh)
+        geomsfh.write(40 * "-" + "\n")
     else:
         try:
             path = "geoms.xyz"
@@ -123,31 +145,33 @@ if __name__ == '__main__':
         except EnvironmentError as err:
             sys.exit(err)
 
-	internal.outfile.write("\nStarting the iterative process.\n")
+    player.outfile.write("\nStarting the iterative process.\n")
 
     ## Initial position (in Bohr)
-	position = internal.system.molecule[0].read_position()
+    position = player.system.molecule[0].read_position()
 
     ## If restarting, read the last gradient and hessian
-	# if internal.player.initcyc > 1:
-	# 	if internal.player.qmprog in ("g03", "g09", "g16"):
+    # if player.player.initcyc > 1:
+    # 	if player.player.qmprog in ("g03", "g09", "g16"):
     # 		Gaussian.read_forces("grad_hessian.dat")
     # 		Gaussian.read_hessian_fchk("grad_hessian.dat")
 
-		#if player['qmprog'] == "molcas":
-			#Molcas.read_forces("grad_hessian.dat")
-			#Molcas.read_hessian("grad_hessian.dat")
+    # if player['qmprog'] == "molcas":
+    # Molcas.read_forces("grad_hessian.dat")
+    # Molcas.read_hessian("grad_hessian.dat")
 
     ###
     ###  Start the iterative process
     ###
 
-	internal.outfile.write("\n" + 90 * "-" + "\n")
+    player.outfile.write("\n" + 90 * "-" + "\n")
 
-	for cycle in range(internal.player.initcyc, internal.player.initcyc + internal.player.maxcyc):
+    for cycle in range(
+        player.player.initcyc, player.player.initcyc + player.player.maxcyc
+    ):
 
-		internal.outfile.write("{} Step # {}\n".format(40 * " ", cycle))
-		internal.outfile.write(90 * "-" + "\n\n")
+        player.outfile.write("{} Step # {}\n".format(40 * " ", cycle))
+        player.outfile.write(90 * "-" + "\n\n")
 
         make_step_dir(cycle)
 
@@ -155,160 +179,42 @@ if __name__ == '__main__':
         ####  Start block of parallel simulations
         ####
 
-		internal.dice_start(cycle)
+        player.dice_start(cycle)
 
         ###
         ###  End of parallel simulations block
         ###
 
-		## Make ASEC
-		# internal.outfile.write("\nBuilding the ASEC and vdW meanfields... ")
-		# asec_charges = internal.populate_asec_vdw(cycle)
+        # Make ASEC
+        player.outfile.write("\nBuilding the ASEC and vdW meanfields... ")
+        asec_charges = player.populate_asec_vdw(cycle)
 
-		# ## After ASEC is built, compress files bigger than 1MB
-		# for proc in range(1, internal.player.nprocs + 1):
-		# 	path = "step{:02d}".format(cycle) + os.sep + "p{:02d}".format(proc)
-		# 	compress_files_1mb(path)
+        ## After ASEC is built, compress files bigger than 1MB
+        for proc in range(1, player.player.nprocs + 1):
+            path = "step{:02d}".format(cycle) + os.sep + "p{:02d}".format(proc)
+            compress_files_1mb(path)
 
         ###
         ###  Start QM calculation
         ###
 
-		# make_qm_dir(cycle)
+        player.gaussian_start(cycle, geomsfh)
 
-		# if internal.player.qmprog in ("g03", "g09", "g16"):
+        player.system.print_geom(cycle, geomsfh)
+        geomsfh.write(40 * "-" + "\n")
 
-		# 	if cycle > 1:
+        player.outfile.write("\n+" + 88 * "-" + "+\n")
 
-		# 		src = "simfiles" + os.sep + "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk"
-		# 		dst = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk"
-		# 		shutil.copyfile(src, dst)
-			
-		# 	internal.make_gaussian_input(cycle)
-		# 	internal.gaussian.run_gaussian(cycle, "force", internal.outfile)
-		# 	internal.gaussian.run_formchk(cycle, internal.outfile)
-				
-		# 	## Read the gradient
-		# 	file = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.fchk"
-
-		# 	try: 
-		# 		gradient
-		# 		old_gradient = gradient
-		# 	except:
-		# 		pass
-
-		# 	gradient = internal.read_forces_fchk(file, internal.outfile)
-
-		# 	# If 1st step, read the hessian
-		# 	if cycle == 1:
-
-		# 		if internal.player.readhessian == "yes":
-
-		# 			file = "grad_hessian.dat"
-		# 			internal.outfile.write("\nReading the hessian matrix from file {}\n".format(file))
-		# 			hessian = internal.read_hessian_log(file)
-
-		# 		else:
-
-		# 			file = "simfiles" + os.sep + "step01" + os.sep + "qm" + os.sep + "asec.fchk"
-		# 			internal.outfile.write("\nReading the hessian matrix from file {}\n".format(file))
-		# 			hessian = internal.read_hessian_fchk(file)
-
-		# 	# From 2nd step on, update the hessian
-		# 	else:
-		# 		internal.outfile.write("\nUpdating the hessian matrix using the BFGS method... ")
-		# 		hessian = internal.system.molecule[0].update_hessian(step, gradient, old_gradient, hessian)
-		# 		internal.outfile.write("Done\n")
-				
-		# 	# Save gradient and hessian
-		# 	internal.print_grad_hessian(cycle, gradient, hessian)
-
-		# 	# Calculate the step and update the position
-		# 	step = internal.calculate_step(cycle, gradient, hessian)
-
-		# 	position += step
-				
-		# 	# ## Update the geometry of the reference molecule
-		# 	internal.system.update_molecule(position, internal.outfile)
-				
-# 				## If needed, calculate the charges
-# 				if cycle < internal.player.switchcyc:
-					
-# 					# internal.gaussian.make_charge_input(cycle, asec_charges)
-# 					internal.gaussian.run_gaussian(cycle, "charge", internal.outfile)
-				
-# 				## Read the new charges and update molecules[0]
-# 				if cycle < internal.player.switchcyc:
-# 					file = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log"
-# 					internal.gaussian.read_charges(file, internal.outfile)
-# 				else:
-# 					file = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.log"
-# 					internal.gaussian.read_charges(file, internal.outfile)
-				
-# 				## Print new info for molecule[0]
-# 				internal.outfile.write("\nNew values for molecule type 1:\n\n")
-# 				internal.system.molecule[0].print_mol_info()
-				
-# 				## Print new geometry in geoms.xyz
-# 				internal.system.molecule[0].print_geom(cycle, geomsfh)
-				
-# 			##
-# 			##  Molcas block
-# 			##
-# 			#if player['qmprog'] == "molcas":
-		
-		
-# 		#elif player['opt'] == "ts":
-			
-# 			##
-# 			##  Gaussian block
-# 			##
-# 			#if player['qmprog'] in ("g03", "g09", "g16"):
-				
-				
-				
-# 			##
-# 			##  Molcas block
-# 			##
-# 			#if player['qmprog'] == "molcas":
-				
-		
-# 		else:    ## Only relax the charge distribution
-			
-# 			if internal.player.qmprog in ("g03", "g09", "g16"):
-				
-# 				if cycle > 1:
-# 					src = "simfiles" + os.sep + "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk"
-# 					dst = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk"
-# 					shutil.copyfile(src, dst)
-				
-# 				# internal.gaussian.make_charge_input(cycle, asec_charges)
-# 				internal.gaussian.run_gaussian(cycle, "charge", internal.outfile)
-				
-# 				file = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log"
-# 				internal.read_charges(file)
-				
-# 				## Print new info for molecule[0]
-# 				internal.outfile.write("\nNew values for molecule type 1:\n\n")
-# 				internal.system.molecule[0].print_mol_info()
-				
-# 			#if player['qmprog'] == "molcas":
-				
-		internal.system.print_geom(cycle, geomsfh)
-		geomsfh.write(40 * "-"+"\n")
-
-		internal.outfile.write("\n+" + 88 * "-" + "+\n")
-
-		pickle.dump([cycle,internal.system], open("latest-step.pkl", "wb"))
+        pickle.dump([cycle, player.system], open("latest-step.pkl", "wb"))
     ####
     ####  End of the iterative process
     ####
 
-## imprimir ultimas mensagens, criar um arquivo de potencial para ser usado em eventual
-## continuacao, fechar arquivos (geoms.xyz, run.log, ...)
+    ## imprimir ultimas mensagens, criar um arquivo de potencial para ser usado em eventual
+    ## continuacao, fechar arquivos (geoms.xyz, run.log, ...)
 
-	internal.outfile.write("\nDiceplayer finished normally!\n")
-	internal.outfile.close()
+    player.outfile.write("\nDiceplayer finished normally!\n")
+    player.outfile.close()
 ####
 ####  End of the program
 ####

run.log

@@ -1,169 +0,0 @@
-##########################################################################################
-#############               Welcome to DICEPLAYER version 1.0                #############
-##########################################################################################
-
-Your python version is 3.8.12 (default, Oct 12 2021, 13:58:32) 
-[GCC 7.5.0]
-
-Program started on Tuesday, 22 Feb 2022 at 02:37:19
-
-Environment variables:
-OMP_STACKSIZE = Not set
-
-==========================================================================================
-                         CONTROL variables being used in this run:
-------------------------------------------------------------------------------------------
-
-altsteps = 20000
-freq = no
-ghosts = no
-initcyc = 1
-lps = no
-maxcyc = 3
-maxstep = 0.3
-nprocs = 4
-opt = no
-qmprog = g16
-readhessian = no
-switchcyc = 3
-tol_factor = 1.2
-vdwforces = no
-
-------------------------------------------------------------------------------------------
-                         DICE variables being used in this run:
-------------------------------------------------------------------------------------------
-
-combrule = *
-dens = 0.75
-infile = <_io.TextIOWrapper name='control.in' mode='r' encoding='UTF-8'>
-isave = 1000
-ljname = phb.ljc
-ncores = 3
-nmol = 1 50
-nstep = 20000 20000
-outfile = <_io.TextIOWrapper name='run.log' mode='w' encoding='UTF-8'>
-outname = phb
-press = 1.0
-progname = dice
-randominit = first
-temp = 300.0
-title = Diceplayer run
-upbuf = 360
-
-------------------------------------------------------------------------------------------
-                         GAUSSIAN variables being used in this run:
-------------------------------------------------------------------------------------------
-
-chgmult = 0 1
-keywords = freq=intmodes
-level = MP2/aug-cc-pVDZ
-pop = chelpg
-qmprog = g16
-
-
-==========================================================================================
-
-Program dice found at /usr/local/bin/dice
-Program g16 found at /usr/local/g16/g16
-Program formchk found at /usr/local/g16/formchk
-
-==========================================================================================
-                    Potential parameters from file phb.ljc:
-------------------------------------------------------------------------------------------
-
-Combination rule: *
-Types of molecules: 2
-
-12 atoms in molecule type 1:
----------------------------------------------------------------------------------
-Lbl  AN       X          Y          Z         Charge    Epsilon   Sigma     Mass
----------------------------------------------------------------------------------
-1     6     0.00000    1.40000    0.00000    0.000000   0.11000  3.7500   12.0110
-1     6     1.21244    0.70000    0.00000    0.000000   0.11000  3.7500   12.0110
-1     6     1.21244   -0.70000    0.00000    0.000000   0.11000  3.7500   12.0110
-1     6     0.00000   -1.40000    0.00000    0.000000   0.11000  3.7500   12.0110
-1     6    -1.21244   -0.70000    0.00000    0.000000   0.11000  3.7500   12.0110
-1     6    -1.21244    0.70000    0.00000    0.000000   0.11000  3.7500   12.0110
-2     1     0.00000    2.48810    0.00000    0.000000   0.00000  0.0000    1.0079
-2     1     2.15467    1.24400    0.00000    0.000000   0.00000  0.0000    1.0079
-2     1     2.15467   -1.24400    0.00000    0.000000   0.00000  0.0000    1.0079
-2     1     0.00000   -2.48810    0.00000    0.000000   0.00000  0.0000    1.0079
-2     1    -2.15467   -1.24400    0.00000    0.000000   0.00000  0.0000    1.0079
-2     1    -2.15467    1.24400    0.00000    0.000000   0.00000  0.0000    1.0079
-
-16 atoms in molecule type 2:
----------------------------------------------------------------------------------
-Lbl  AN       X          Y          Z         Charge    Epsilon   Sigma     Mass
----------------------------------------------------------------------------------
-1     6     0.67203   -2.82345    0.00263   -0.115000   0.07000  3.5500   12.0110
-1     6     2.07203   -2.82345    0.00263   -0.115000   0.07000  3.5500   12.0110
-1     6     2.76823   -1.61764    0.00263   -0.115000   0.07000  3.5500   12.0110
-1     6     2.06824   -0.40521    0.00264   -0.115000   0.07000  3.5500   12.0110
-1    14     0.67589   -0.40522    0.00264    0.150000   0.07000  3.5500   28.0860
-1    10    -0.02420   -1.61760    0.00263   -0.115000   0.07000  3.5500   20.1800
-2     1     0.13203   -3.75875    0.00263    0.115000   0.03000  2.4200    1.0079
-2     1     2.61203   -3.75875    0.00263    0.115000   0.03000  2.4200    1.0079
-2     1     2.60824    0.53010    0.00264    0.115000   0.03000  2.4200    1.0079
-2     1    -1.10420   -1.61760    0.00263    0.115000   0.03000  2.4200    1.0079
-3     8    -0.00411    0.77257    0.00264   -0.585000   0.17000  3.0700   15.9990
-3     1     0.61978    1.50220    0.00264    0.435000   0.00000  0.0000    1.0079
-4     6     4.27823   -1.61764    0.00263    0.115000   0.17000  3.8000   12.0110
-4     1     4.63490   -0.74399    0.50704    0.000000   0.00000  0.0000    1.0079
-4     1     4.63490   -2.49130    0.50704    0.000000   0.00000  0.0000    1.0079
-4     1     4.63490   -1.61764   -1.00617    0.000000   0.00000  0.0000    1.0079
-
-
-==========================================================================================
-
-Molecule type 1 - BENZENO:
-
-    Center of mass = (     0.0000 ,    -0.0000 ,     0.0000 )
-    Moments of inertia =  8.934187E+01  8.934282E+01  1.786847E+02
-    Major principal axis = (   0.000000 ,   1.000000 ,   0.000000 )
-    Inter principal axis = (   1.000000 ,   0.000000 ,   0.000000 )
-    Minor principal axis = (   0.000000 ,   0.000000 ,   1.000000 )
-    Characteristic lengths = (   4.31 ,   4.98 ,   0.00 )
-    Total mass =    78.11 au
-    Total charge =   0.0000 e
-    Dipole moment = (    0.0000 ,    0.0000 ,    0.0000 )     Total =    0.0000 Debye
-
-    Translating and rotating molecule to standard orientation... Done
-
-    New values:
-    Center of mass = (     0.0000 ,    -0.0000 ,     0.0000 )
-    Moments of inertia =  8.934187E+01  8.934282E+01  1.786847E+02
-    Major principal axis = (   1.000000 ,   0.000000 ,   0.000000 )
-    Inter principal axis = (   0.000000 ,   1.000000 ,   0.000000 )
-    Minor principal axis = (   0.000000 ,   0.000000 ,   1.000000 )
-    Characteristic lengths = (   4.98 ,   4.31 ,   0.00 )
-    Total mass =    78.11 au
-    Total charge =   0.0000 e
-    Dipole moment = (    0.0000 ,    0.0000 ,    0.0000 )     Total =    0.0000 Debye
-
-
-Molecule type 2 - PLACEHOLDER:
-
-    Center of mass = (     1.3581 ,    -1.1728 ,     0.0026 )
-    Moments of inertia =  1.506334E+02  3.060213E+02  4.535775E+02
-    Major principal axis = (   0.879376 ,  -0.476127 ,  -0.000001 )
-    Inter principal axis = (   0.476127 ,   0.879376 ,   0.000001 )
-    Minor principal axis = (   0.000000 ,  -0.000002 ,   1.000000 )
-    Characteristic lengths = (   5.74 ,   5.26 ,   1.51 )
-    Total mass =   132.38 au
-    Total charge =  -0.0000 e
-    Dipole moment = (    2.3293 ,    0.1593 ,   -0.0000 )     Total =    2.3347 Debye
-
-    Translating and rotating molecule to standard orientation... Done
-
-    New values:
-    Center of mass = (     0.0000 ,    -0.0000 ,     0.0000 )
-    Moments of inertia =  1.506334E+02  3.060213E+02  4.535775E+02
-    Major principal axis = (   1.000000 ,   0.000000 ,   0.000000 )
-    Inter principal axis = (   0.000000 ,   1.000000 ,   0.000000 )
-    Minor principal axis = (   0.000000 ,   0.000000 ,   1.000000 )
-    Characteristic lengths = (   5.63 ,   4.95 ,   1.51 )
-    Total mass =   132.38 au
-    Total charge =  -0.0000 e
-    Dipole moment = (    1.9725 ,    1.2491 ,   -0.0000 )     Total =    2.3347 Debye
-
-==========================================================================================

run.log.backup

@@ -1,169 +0,0 @@
-##########################################################################################
-#############               Welcome to DICEPLAYER version 1.0                #############
-##########################################################################################
-
-Your python version is 3.8.12 | packaged by conda-forge | (default, Oct 12 2021, 21:57:06) 
-[GCC 9.4.0]
-
-Program started on Tuesday, 22 Feb 2022 at 02:04:46
-
-Environment variables:
-OMP_STACKSIZE = Not set
-
-==========================================================================================
-                         CONTROL variables being used in this run:
-------------------------------------------------------------------------------------------
-
-altsteps = 20000
-freq = no
-ghosts = no
-initcyc = 1
-lps = no
-maxcyc = 3
-maxstep = 0.3
-nprocs = 4
-opt = no
-qmprog = g16
-readhessian = no
-switchcyc = 3
-tol_factor = 1.2
-vdwforces = no
-
-------------------------------------------------------------------------------------------
-                         DICE variables being used in this run:
-------------------------------------------------------------------------------------------
-
-combrule = *
-dens = 0.75
-infile = <_io.TextIOWrapper name='control.in' mode='r' encoding='UTF-8'>
-isave = 1000
-ljname = phb.ljc
-ncores = 3
-nmol = 1 50
-nstep = 20000 20000
-outfile = <_io.TextIOWrapper name='run.log' mode='w' encoding='UTF-8'>
-outname = phb
-press = 1.0
-progname = dice
-randominit = first
-temp = 300.0
-title = Diceplayer run
-upbuf = 360
-
-------------------------------------------------------------------------------------------
-                         GAUSSIAN variables being used in this run:
-------------------------------------------------------------------------------------------
-
-chgmult = 0 1
-keywords = freq=intmodes
-level = MP2/aug-cc-pVDZ
-pop = chelpg
-qmprog = g16
-
-
-==========================================================================================
-
-Program dice found at /usr/local/bin/dice
-Program g16 found at /usr/local/g16/g16
-Program formchk found at /usr/local/g16/formchk
-
-==========================================================================================
-                    Potential parameters from file phb.ljc:
-------------------------------------------------------------------------------------------
-
-Combination rule: *
-Types of molecules: 2
-
-12 atoms in molecule type 1:
----------------------------------------------------------------------------------
-Lbl  AN       X          Y          Z         Charge    Epsilon   Sigma     Mass
----------------------------------------------------------------------------------
-1     6     0.00000    1.40000    0.00000    0.000000   0.11000  3.7500   12.0110
-1     6     1.21244    0.70000    0.00000    0.000000   0.11000  3.7500   12.0110
-1     6     1.21244   -0.70000    0.00000    0.000000   0.11000  3.7500   12.0110
-1     6     0.00000   -1.40000    0.00000    0.000000   0.11000  3.7500   12.0110
-1     6    -1.21244   -0.70000    0.00000    0.000000   0.11000  3.7500   12.0110
-1     6    -1.21244    0.70000    0.00000    0.000000   0.11000  3.7500   12.0110
-2     1     0.00000    2.48810    0.00000    0.000000   0.00000  0.0000    1.0079
-2     1     2.15467    1.24400    0.00000    0.000000   0.00000  0.0000    1.0079
-2     1     2.15467   -1.24400    0.00000    0.000000   0.00000  0.0000    1.0079
-2     1     0.00000   -2.48810    0.00000    0.000000   0.00000  0.0000    1.0079
-2     1    -2.15467   -1.24400    0.00000    0.000000   0.00000  0.0000    1.0079
-2     1    -2.15467    1.24400    0.00000    0.000000   0.00000  0.0000    1.0079
-
-16 atoms in molecule type 2:
----------------------------------------------------------------------------------
-Lbl  AN       X          Y          Z         Charge    Epsilon   Sigma     Mass
----------------------------------------------------------------------------------
-1     6     0.67203   -2.82345    0.00263   -0.115000   0.07000  3.5500   12.0110
-1     6     2.07203   -2.82345    0.00263   -0.115000   0.07000  3.5500   12.0110
-1     6     2.76823   -1.61764    0.00263   -0.115000   0.07000  3.5500   12.0110
-1     6     2.06824   -0.40521    0.00264   -0.115000   0.07000  3.5500   12.0110
-1    14     0.67589   -0.40522    0.00264    0.150000   0.07000  3.5500   28.0860
-1    10    -0.02420   -1.61760    0.00263   -0.115000   0.07000  3.5500   20.1800
-2     1     0.13203   -3.75875    0.00263    0.115000   0.03000  2.4200    1.0079
-2     1     2.61203   -3.75875    0.00263    0.115000   0.03000  2.4200    1.0079
-2     1     2.60824    0.53010    0.00264    0.115000   0.03000  2.4200    1.0079
-2     1    -1.10420   -1.61760    0.00263    0.115000   0.03000  2.4200    1.0079
-3     8    -0.00411    0.77257    0.00264   -0.585000   0.17000  3.0700   15.9990
-3     1     0.61978    1.50220    0.00264    0.435000   0.00000  0.0000    1.0079
-4     6     4.27823   -1.61764    0.00263    0.115000   0.17000  3.8000   12.0110
-4     1     4.63490   -0.74399    0.50704    0.000000   0.00000  0.0000    1.0079
-4     1     4.63490   -2.49130    0.50704    0.000000   0.00000  0.0000    1.0079
-4     1     4.63490   -1.61764   -1.00617    0.000000   0.00000  0.0000    1.0079
-
-
-==========================================================================================
-
-Molecule type 1 - BENZENO:
-
-    Center of mass = (     0.0000 ,    -0.0000 ,     0.0000 )
-    Moments of inertia =  8.934187E+01  8.934282E+01  1.786847E+02
-    Major principal axis = (   0.000000 ,   1.000000 ,   0.000000 )
-    Inter principal axis = (   1.000000 ,   0.000000 ,   0.000000 )
-    Minor principal axis = (   0.000000 ,   0.000000 ,   1.000000 )
-    Characteristic lengths = (   4.31 ,   4.98 ,   0.00 )
-    Total mass =    78.11 au
-    Total charge =   0.0000 e
-    Dipole moment = (    0.0000 ,    0.0000 ,    0.0000 )     Total =    0.0000 Debye
-
-    Translating and rotating molecule to standard orientation... Done
-
-    New values:
-    Center of mass = (     0.0000 ,    -0.0000 ,     0.0000 )
-    Moments of inertia =  8.934187E+01  8.934282E+01  1.786847E+02
-    Major principal axis = (   1.000000 ,   0.000000 ,   0.000000 )
-    Inter principal axis = (   0.000000 ,   1.000000 ,   0.000000 )
-    Minor principal axis = (   0.000000 ,   0.000000 ,   1.000000 )
-    Characteristic lengths = (   4.98 ,   4.31 ,   0.00 )
-    Total mass =    78.11 au
-    Total charge =   0.0000 e
-    Dipole moment = (    0.0000 ,    0.0000 ,    0.0000 )     Total =    0.0000 Debye
-
-
-Molecule type 2 - PLACEHOLDER:
-
-    Center of mass = (     1.3581 ,    -1.1728 ,     0.0026 )
-    Moments of inertia =  1.506334E+02  3.060213E+02  4.535775E+02
-    Major principal axis = (   0.879376 ,  -0.476127 ,  -0.000001 )
-    Inter principal axis = (   0.476127 ,   0.879376 ,   0.000001 )
-    Minor principal axis = (   0.000000 ,  -0.000002 ,   1.000000 )
-    Characteristic lengths = (   5.74 ,   5.26 ,   1.51 )
-    Total mass =   132.38 au
-    Total charge =  -0.0000 e
-    Dipole moment = (    2.3293 ,    0.1593 ,   -0.0000 )     Total =    2.3347 Debye
-
-    Translating and rotating molecule to standard orientation... Done
-
-    New values:
-    Center of mass = (     0.0000 ,    -0.0000 ,     0.0000 )
-    Moments of inertia =  1.506334E+02  3.060213E+02  4.535775E+02
-    Major principal axis = (   1.000000 ,   0.000000 ,   0.000000 )
-    Inter principal axis = (   0.000000 ,   1.000000 ,   0.000000 )
-    Minor principal axis = (   0.000000 ,   0.000000 ,   1.000000 )
-    Characteristic lengths = (   5.63 ,   4.95 ,   1.51 )
-    Total mass =   132.38 au
-    Total charge =  -0.0000 e
-    Dipole moment = (    1.9725 ,    1.2491 ,   -0.0000 )     Total =    2.3347 Debye
-
-==========================================================================================