From 4ae8385918c93a6eeb3fe7e3281a420a9dfccacb Mon Sep 17 00:00:00 2001 From: Vitor Hideyoshi Date: Tue, 31 May 2022 09:13:08 -0300 Subject: [PATCH] Initial Folder Rework Implementation Adds the Environment, External, Utils folder inside de DPpack. All classes are going to be implemented there --- control.in | 22 - control.yml | 22 + diceplayer/DPpack/Environment/Atom.py | 52 + diceplayer/DPpack/Environment/Molecule.py | 418 +++++++ diceplayer/DPpack/Environment/System.py | 224 ++++ diceplayer/DPpack/Environment/__init__.py | 0 diceplayer/DPpack/External/Dice.py | 786 +++++++++++++ diceplayer/DPpack/External/Gaussian.py | 688 +++++++++++ diceplayer/DPpack/External/__init__.py | 0 diceplayer/DPpack/MolHandling.py | 460 -------- diceplayer/DPpack/Molcas.py | 348 ------ diceplayer/DPpack/Optimization.py | 266 ----- diceplayer/DPpack/Player.py | 830 ++++++++++++++ diceplayer/DPpack/SetGlobals.py | 1261 --------------------- diceplayer/DPpack/{ => Utils}/Misc.py | 0 diceplayer/DPpack/Utils/Optimization.py | 263 +++++ diceplayer/DPpack/{ => Utils}/PTable.py | 0 diceplayer/DPpack/Utils/Validations.py | 15 + diceplayer/DPpack/Utils/__init__.py | 0 diceplayer/DPpack/external/Dice.py | 662 ----------- diceplayer/__main__.py | 414 +++---- run.log | 169 --- run.log.backup | 169 --- 23 files changed, 3458 insertions(+), 3611 deletions(-) delete mode 100644 control.in create mode 100644 control.yml create mode 100644 diceplayer/DPpack/Environment/Atom.py create mode 100644 diceplayer/DPpack/Environment/Molecule.py create mode 100644 diceplayer/DPpack/Environment/System.py create mode 100644 diceplayer/DPpack/Environment/__init__.py create mode 100644 diceplayer/DPpack/External/Dice.py create mode 100644 diceplayer/DPpack/External/Gaussian.py create mode 100644 diceplayer/DPpack/External/__init__.py delete mode 100644 diceplayer/DPpack/MolHandling.py delete mode 100644 diceplayer/DPpack/Molcas.py delete mode 100644 diceplayer/DPpack/Optimization.py create mode 100644 diceplayer/DPpack/Player.py delete mode 100644 diceplayer/DPpack/SetGlobals.py rename diceplayer/DPpack/{ => Utils}/Misc.py (100%) create mode 100644 diceplayer/DPpack/Utils/Optimization.py rename diceplayer/DPpack/{ => Utils}/PTable.py (100%) create mode 100644 diceplayer/DPpack/Utils/Validations.py create mode 100644 diceplayer/DPpack/Utils/__init__.py delete mode 100644 diceplayer/DPpack/external/Dice.py delete mode 100644 run.log delete mode 100644 run.log.backup diff --git a/control.in b/control.in deleted file mode 100644 index 28820ac..0000000 --- a/control.in +++ /dev/null @@ -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 \ No newline at end of file diff --git a/control.yml b/control.yml new file mode 100644 index 0000000..5d14715 --- /dev/null +++ b/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' \ No newline at end of file diff --git a/diceplayer/DPpack/Environment/Atom.py b/diceplayer/DPpack/Environment/Atom.py new file mode 100644 index 0000000..1da6b06 --- /dev/null +++ b/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] diff --git a/diceplayer/DPpack/Environment/Molecule.py b/diceplayer/DPpack/Environment/Molecule.py new file mode 100644 index 0000000..9f5c3bb --- /dev/null +++ b/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) diff --git a/diceplayer/DPpack/Environment/System.py b/diceplayer/DPpack/Environment/System.py new file mode 100644 index 0000000..9e5d52a --- /dev/null +++ b/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 + ) + ) diff --git a/diceplayer/DPpack/Environment/__init__.py b/diceplayer/DPpack/Environment/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/diceplayer/DPpack/External/Dice.py b/diceplayer/DPpack/External/Dice.py new file mode 100644 index 0000000..db6e259 --- /dev/null +++ b/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 diff --git a/diceplayer/DPpack/External/Gaussian.py b/diceplayer/DPpack/External/Gaussian.py new file mode 100644 index 0000000..84e11e6 --- /dev/null +++ b/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 diff --git a/diceplayer/DPpack/External/__init__.py b/diceplayer/DPpack/External/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/diceplayer/DPpack/MolHandling.py b/diceplayer/DPpack/MolHandling.py deleted file mode 100644 index 6c4509f..0000000 --- a/diceplayer/DPpack/MolHandling.py +++ /dev/null @@ -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)) - - - diff --git a/diceplayer/DPpack/Molcas.py b/diceplayer/DPpack/Molcas.py deleted file mode 100644 index 4d8a39a..0000000 --- a/diceplayer/DPpack/Molcas.py +++ /dev/null @@ -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 - - - diff --git a/diceplayer/DPpack/Optimization.py b/diceplayer/DPpack/Optimization.py deleted file mode 100644 index 26c0de4..0000000 --- a/diceplayer/DPpack/Optimization.py +++ /dev/null @@ -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 - - - diff --git a/diceplayer/DPpack/Player.py b/diceplayer/DPpack/Player.py new file mode 100644 index 0000000..9aeb8d0 --- /dev/null +++ b/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 diff --git a/diceplayer/DPpack/SetGlobals.py b/diceplayer/DPpack/SetGlobals.py deleted file mode 100644 index 1495eb9..0000000 --- a/diceplayer/DPpack/SetGlobals.py +++ /dev/null @@ -1,1261 +0,0 @@ -from diceplayer.DPpack.MolHandling import * -from diceplayer.DPpack.PTable import * -from diceplayer.DPpack.Misc import * - -from diceplayer.DPpack.external.Dice import * - -from typing import IO, Tuple, List, TextIO, Union - -from numpy.core.numeric import partition -from numpy import random - -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 Internal: - - 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 = self.Gaussian() - self.gaussian_keywords = [a for a in dir(self.gaussian) if not a.startswith( - '__') and not callable(getattr(self.gaussian, a))] - - # self.molcas = self.Molcas() - # self.molcas_keywords = [a for a in dir(self.molcas) if not a.startswith('__') and not callable(getattr(self.molcas, a))] - - # Constanst that shall be set for global use - - self.tol_rms_force = 3e-4 # Hartree/Bohr - self.tol_max_force = 4.5e-4 # Hartree/Bohr - self.tol_rms_step = 1.2e-3 # Bohr - self.tol_max_step = 1.8e-3 # Bohr - self.trust_radius = None - - # Dice: - self.combrule = None - - def read_keywords(self) -> None: - - try: - controlfile = self.infile.readlines() - except EnvironmentError: - sys.exit("Error: cannot read file {}".format(self.infile)) - - for line in controlfile: - - key, value = line.partition("=")[::2] # Discards the '=' - key = key.strip().lower() - if key in ('title', 'keywords'): - value = value.strip() - else: - value = value.split() - - # Read the Diceplayer related keywords - # 'value' is not empty! - if key in self.player_keywords and len(value) != 0: - - if key == 'qmprog' and value[0].lower() in ("g03", "g09", "g16", "molcas"): - - setattr(self.player, key, value[0].lower()) - - if self.player.qmprog in ("g03", "g09", "g16"): - - self.gaussian.qmprog = self.player.qmprog - - # if self.player.qmprog == "molcas": - - # pass - - elif key == 'opt' and value[0].lower() in ("yes", "no", "ts"): - setattr(self.player, key, value[0].lower()) - - # elif key == 'zipprog' and value[0].lower() in ("zip", "gzip", "bzip"): - #player[key] = value[0].lower() - - elif key in ('lps', 'ghosts') and value[0].lower() in ("yes", "no"): - setattr(self.player, key, value[0].lower()) - - elif key in ('readhessian', 'vdwforces') and value[0].lower() in ("yes", "no"): - setattr(self.player, key, value[0].lower()) - - elif key in ('maxcyc', 'nprocs', 'altsteps', 'switchcyc'): - err = "Error: expected a positive integer for keyword {} in file {}".format( - key, self.infile) - try: - new_value = int(value[0]) - if new_value >= 1: - setattr(self.player, key, new_value) - elif key == 'altsteps' and new_value == 0: - setattr(self.player, key, 0) - except ValueError: - sys.exit(err) - - elif key == 'maxstep': # Cannot be less than 0.01 - err = "Error: expected a float greater than 0.01 for keyword {} in file {}".format( - key, self.infile) - try: - new_value = float(value[0]) - if new_value < 0.01: - sys.exit(err) - else: - setattr(self.player, key, new_value) - except ValueError: - sys.exit(err) - - # Read the Dice related keywords - # 'value' is not empty! - elif key in self.dice_keywords and len(value) != 0: - - if key == 'title': - setattr(self.dice, key, value) - - elif key in ('ljname', 'outname', 'progname'): - setattr(self.dice, key, value[0]) - - elif key == 'randominit': - if value in ('always', 'first'): - setattr(self.dice, key, value[0]) - - elif key in ('ncores', 'isave'): - err = "Error: expected a positive integer for keyword {} in file {}".format( - key, self.infile) - if not value[0].isdigit(): - sys.exit(err) - new_value = int(value[0]) - if new_value >= 1: - setattr(self.dice, key, new_value) - - elif key in ('temp', 'press', 'dens'): # Cannot be less than 1e-10 - err = "Error: expected a positive float for keyword {} in file {}".format( - key, self.infile) - try: - new_value = float(value[0]) - if new_value < 1e-10: - sys.exit(err) - else: - setattr(self.dice, key, new_value) - except ValueError: - sys.exit(err) - - # If defined, must be well defined (only positive integer values) - elif key == 'nmol': - err = "Error: expected 1 to 4 positive integers for keyword {} in file {}".format( - key, self.infile) - args = min(4, len(value)) - for i in range(args): - if value[i].isdigit(): - new_value = int(value[i]) - if new_value < 1: - sys.exit(err) - else: - getattr(self.dice, key).append(new_value) - elif i == 0: - sys.exit(err) - else: - break - - # If defined, must be well defined (only positive integer values) - elif key == 'nstep': - err = "Error: expected 2 or 3 positive integers for keyword {} in file {}".format( - key, self.infile) - if len(value) < 2: - sys.exit(err) - args = min(3, len(value)) - for i in range(args): - if value[i].isdigit(): - new_value = int(value[i]) - if new_value < 1: - sys.exit(err) - else: - getattr(self.dice, key).append(new_value) - elif i < 2: - sys.exit(err) - else: - break - - # Read the Gaussian related keywords - # 'value' is not empty! - elif key in self.gaussian_keywords and len(value) != 0: - - if key == 'mem': # Memory in MB (minimum of 100) - err = "Error: expected a positive integer for keyword {} in file {}".format( - key, self.infile) - if not value[0].isdigit(): - sys.exit(err) - new_value = int(value[0]) - if new_value >= 100: - setattr(self.gaussian, key, new_value) - - elif key == 'keywords': - setattr(self.gaussian, key, value) - - # If defined, must be well defined (2 integer values) - elif key == 'chgmult': - err = "Error: expected 2 integers for keyword {} in file {}".format( - key, self.infile) - if len(value) < 2: - sys.exit(err) - for i in range(2): - try: - getattr(self.gaussian, key)[i] = int(value[i]) - except ValueError: - sys.exit(err) - - elif key == 'level': - setattr(self.gaussian, key, value[0]) - - elif key in ('gmiddle', 'gbottom'): - setattr(self.gaussian, key, value[0]) - - elif key == 'pop' and value[0].lower() in ("chelpg", "mk", "nbo"): - setattr(self.gaussian, key, value[0].lower()) - - # #### Read the Molcas related keywords - # elif key in self.molcas_keywords and len(value) != 0: ## 'value' is not empty! - - # if key == 'root': # If defined, must be well defined (only positive integer values) - # err = "Error: expected a positive integer for keyword {} in file {}".format(key, self.infile) - # if not value[0].isdigit(): - # sys.exit(err) - # new_value = int(value[0]) - # if new_value >= 1: - # setattr(self.molcas, key, new_value) - - # elif key in ('mbottom', 'orbfile'): - # setattr(self.molcas, key, value[0]) - - # elif key == 'basis': - # setattr(self.molcas ,key, value[0]) - - # #### End - - 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 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 - - # 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() # (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) + 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].molnale) - # 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: - 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 - - # ---------------------------------------------------------------------------------------------------------------------------------------- - # Gaussian related methods - # ---------------------------------------------------------------------------------------------------------------------------------------- - 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.system.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.system.molecule[0].atom)): - fh.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format( - i + 1, self.system.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.system.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.system.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.system.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.gaussian.mem != None: - fh.write("%Mem={}MB\n".format(self.gaussian.mem)) - fh.write("%Nprocs={}\n".format(self.player.nprocs * self.dice.ncores)) - - kword_line = "#P " + str(self.gaussian.level) - - if self.gaussian.keywords != None: - kword_line += " " + self.gaussian.keywords - - if self.player.opt == 'yes': - kword_line += " Force" - - # kword_line += " Charge" - kword_line += " NoSymm" - kword_line += " Pop={} Density=Current".format(self.gaussian.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.gaussian.chgmult[0], self.gaussian.chgmult[1])) - - for atom in self.system.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.gaussian.gmiddle != None: - # if not os.path.isfile(self.gaussian.gmiddle): - # sys.exit("Error: cannot find file {} in main directory".format( - # self.gaussian.gmiddle)) - # try: - # with open(self.gaussian.gmiddle) as gmiddlefile: - # gmiddle = gmiddlefile.readlines() - # except: - # sys.exit("Error: cannot open file {}".format(self.gaussian.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.gaussian.gbottom != None: - # if not os.path.isfile(self.gaussian.gbottom): - # sys.exit("Error: cannot find file {} in main directory".format( - # self.gaussian.gbottom)) - # try: - # with open(self.gaussian.gbottom) as gbottomfile: - # gbottom = gbottomfile.readlines() - # except: - # sys.exit("Error: cannot open file {}".format(self.gaussian.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.system.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.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") - - class Player: - - def __init__(self) -> None: - - self.maxcyc = None - self.nprocs = 1 - self.switchcyc = 3 - self.altsteps = 20000 - self.maxstep = .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 - - class Gaussian: - - def __init__(self) -> None: - - self.qmprog = "g09" - self.path = None - - self.mem = None - self.keywords = None - self.chgmult = [0, 1] - self.gmiddle = None # In each case, if a filename is given, its content will be - self.gbottom = None # inserted in the gaussian input - self.pop = "chelpg" - self.level = None - - 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 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) - - # class Molcas: - - # def __init(self): - - # self.orbfile = "input.exporb" - # self.root = 1 - - # self.mbottom = None - # self.basis = None diff --git a/diceplayer/DPpack/Misc.py b/diceplayer/DPpack/Utils/Misc.py similarity index 100% rename from diceplayer/DPpack/Misc.py rename to diceplayer/DPpack/Utils/Misc.py diff --git a/diceplayer/DPpack/Utils/Optimization.py b/diceplayer/DPpack/Utils/Optimization.py new file mode 100644 index 0000000..cc8cf91 --- /dev/null +++ b/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 diff --git a/diceplayer/DPpack/PTable.py b/diceplayer/DPpack/Utils/PTable.py similarity index 100% rename from diceplayer/DPpack/PTable.py rename to diceplayer/DPpack/Utils/PTable.py diff --git a/diceplayer/DPpack/Utils/Validations.py b/diceplayer/DPpack/Utils/Validations.py new file mode 100644 index 0000000..5381863 --- /dev/null +++ b/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 diff --git a/diceplayer/DPpack/Utils/__init__.py b/diceplayer/DPpack/Utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/diceplayer/DPpack/external/Dice.py b/diceplayer/DPpack/external/Dice.py deleted file mode 100644 index 7e1c367..0000000 --- a/diceplayer/DPpack/external/Dice.py +++ /dev/null @@ -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 diff --git a/diceplayer/__main__.py b/diceplayer/__main__.py index 2fb788e..ae0e82b 100644 --- a/diceplayer/__main__.py +++ b/diceplayer/__main__.py @@ -8,307 +8,213 @@ 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]') - ## Study the option of a parameter for continuing the last process via data from control.in and run.log files + 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() + 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: + try: - if args.opt_continue and os.path.exists(args.outfile): - - save = pickle.load(open("latest-step.pkl","rb")) + if args.opt_continue and os.path.exists(args.outfile): - if os.path.isfile(args.outfile+".backup"): - os.remove(args.outfile+".backup") + save = pickle.load(open("latest-step.pkl", "rb")) - os.rename(args.outfile,args.outfile+".backup") - outfile = open(args.outfile,'w',1) + if os.path.isfile(args.outfile + ".backup"): + os.remove(args.outfile + ".backup") - elif os.path.exists(args.outfile): - os.rename(args.outfile, args.outfile+".backup") - outfile = open(args.outfile,'w',1) - else: - outfile = open(args.outfile,"w",1) + os.rename(args.outfile, args.outfile + ".backup") + outfile = open(args.outfile, "w", 1) - except Exception as err: - sys.exit(err) + elif os.path.exists(args.outfile): + os.rename(args.outfile, args.outfile + ".backup") + outfile = open(args.outfile, "w", 1) + else: + outfile = open(args.outfile, "w", 1) - try: + except Exception as err: + sys.exit(err) - if os.path.exists(args.infile): - infile = open(args.infile,"r") + try: - except Exception as err: - sys.exit(err) + if os.path.exists(args.infile): + infile = open(args.infile, "r") -#### Read and check the keywords in INFILE + except Exception as err: + sys.exit(err) - internal = Internal(infile, outfile) + #### Read and check the keywords in INFILE - internal.read_keywords() + player = Player(infile, outfile) - internal.check_keywords() - internal.print_keywords() + player.read_keywords() - if args.opt_continue: - internal.player.initcyc = save[0] + 1 - internal.system = save[1] - else: - internal.read_potential() + player.check_keywords() + player.print_keywords() -#### Check whether the executables are in the path -#### and print potential to Log File + if args.opt_continue: + player.player.initcyc = save[0] + 1 + player.system = save[1] + else: + player.read_potential() - internal.check_executables() + #### Check whether the executables are in the path + #### and print potential to Log File - internal.print_potential() + player.check_executables() -#### Bring the molecules to standard orientation and prints info about them + player.print_potential() - for i in range(len(internal.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) - - internal.outfile.write(90 * "=") - internal.outfile.write("\n") + #### Bring the molecules to standard orientation and prints info about them - 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) + for i in range(len(player.system.molecule)): -#### Open the geoms.xyz file and prints the initial geometry if starting from zero + 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) - if internal.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") - else: - try: - path = "geoms.xyz" - geomsfh = open(path, "a", 1) - except EnvironmentError as err: - sys.exit(err) + player.outfile.write(90 * "=") + player.outfile.write("\n") - internal.outfile.write("\nStarting the iterative process.\n") - - ## Initial position (in Bohr) - position = internal.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"): - # 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") - - ### - ### Start the iterative process - ### + if not args.opt_continue: + make_simulation_dir() + else: + simdir = "simfiles" + stepdir = "step{:02d}".format(player.player.initcyc) + if os.path.exists(simdir + os.sep + stepdir): + shutil.rmtree(simdir + os.sep + stepdir) - internal.outfile.write("\n" + 90 * "-" + "\n") - - for cycle in range(internal.player.initcyc, internal.player.initcyc + internal.player.maxcyc): - - internal.outfile.write("{} Step # {}\n".format(40 * " ", cycle)) - internal.outfile.write(90 * "-" + "\n\n") + #### Open the geoms.xyz file and prints the initial geometry if starting from zero - make_step_dir(cycle) - - #### - #### Start block of parallel simulations - #### - - internal.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) - - # ## 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) - - ### - ### Start QM calculation - ### - - # make_qm_dir(cycle) + if player.player.initcyc == 1: + try: + path = "geoms.xyz" + geomsfh = open(path, "w", 1) + except EnvironmentError as err: + sys.exit(err) + player.system.print_geom(0, geomsfh) + geomsfh.write(40 * "-" + "\n") + else: + try: + path = "geoms.xyz" + geomsfh = open(path, "a", 1) + except EnvironmentError as err: + sys.exit(err) - # if internal.player.qmprog in ("g03", "g09", "g16"): - - # if cycle > 1: + player.outfile.write("\nStarting the iterative process.\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" + ## Initial position (in Bohr) + position = player.system.molecule[0].read_position() - # try: - # gradient - # old_gradient = gradient - # except: - # pass + ## If restarting, read the last gradient and hessian + # 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") - # gradient = internal.read_forces_fchk(file, internal.outfile) + # if player['qmprog'] == "molcas": + # Molcas.read_forces("grad_hessian.dat") + # Molcas.read_hessian("grad_hessian.dat") - # # If 1st step, read the hessian - # if cycle == 1: + ### + ### Start the iterative process + ### - # if internal.player.readhessian == "yes": + player.outfile.write("\n" + 90 * "-" + "\n") - # file = "grad_hessian.dat" - # internal.outfile.write("\nReading the hessian matrix from file {}\n".format(file)) - # hessian = internal.read_hessian_log(file) + for cycle in range( + player.player.initcyc, player.player.initcyc + player.player.maxcyc + ): - # else: + player.outfile.write("{} Step # {}\n".format(40 * " ", cycle)) + player.outfile.write(90 * "-" + "\n\n") - # 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) + make_step_dir(cycle) - # # 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) + #### + #### Start block of parallel simulations + #### - # # Calculate the step and update the position - # step = internal.calculate_step(cycle, gradient, hessian) + player.dice_start(cycle) - # 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") + ### + ### End of parallel simulations block + ### - internal.outfile.write("\n+" + 88 * "-" + "+\n") + # Make ASEC + player.outfile.write("\nBuilding the ASEC and vdW meanfields... ") + asec_charges = player.populate_asec_vdw(cycle) - pickle.dump([cycle,internal.system], open("latest-step.pkl", "wb")) - #### - #### End of the iterative process - #### + ## 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) -## imprimir ultimas mensagens, criar um arquivo de potencial para ser usado em eventual -## continuacao, fechar arquivos (geoms.xyz, run.log, ...) + ### + ### Start QM calculation + ### - internal.outfile.write("\nDiceplayer finished normally!\n") - internal.outfile.close() + player.gaussian_start(cycle, geomsfh) + + player.system.print_geom(cycle, geomsfh) + geomsfh.write(40 * "-" + "\n") + + player.outfile.write("\n+" + 88 * "-" + "+\n") + + 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, ...) + + player.outfile.write("\nDiceplayer finished normally!\n") + player.outfile.close() #### #### End of the program #### diff --git a/run.log b/run.log deleted file mode 100644 index 40c6461..0000000 --- a/run.log +++ /dev/null @@ -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 - -========================================================================================== diff --git a/run.log.backup b/run.log.backup deleted file mode 100644 index 51ae859..0000000 --- a/run.log.backup +++ /dev/null @@ -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 - -==========================================================================================