HideyoshiNakazone/DicePlayer
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'devel' into main

Browse Source
This commit is contained in:
Vitor Hideyoshi
2022-06-01 20:07:01 -03:00
parent 4d122a43aa e54b3e55fa
commit e253c4fe3e
24 changed files with 3403 additions and 3132 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
.vscode/settings.json vendored
View File

@@ -1,3 +0,0 @@
{
"python.pythonPath": "/home/hideyoshi/Python/3.8/bin/python3.8"
}

22
control.in
View File

@@ -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

23
control.yml.exemple Normal file
View File

@@ -0,0 +1,23 @@
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: [2000, 3000]
isave: 1000
ljname: 'phb.ljc'
outname: 'phb'
randominit: 'always'
gaussian:
qmprog: 'g16'
level: 'MP2/aug-cc-pVDZ'
keywords: 'freq'

52
diceplayer/DPpack/Environment/Atom.py Normal file
View File

@@ -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]

421
diceplayer/DPpack/Environment/Molecule.py Normal file
View File

@@ -0,0 +1,421 @@
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
total_mass (int): The total mass of the molecule
com (NDArray[Any, Any]): The center of mass of the molecule
"""
def __init__(self, molname: str) -> None:
"""
The constructor function __init__ is used to create new instances of the Molecule class.
Args:
molname (str): Molecule name
"""
self.molname: str = molname
self.atom: List[Atom] = []
self.position: NDArray[Any, Any]
self.energy: NDArray[Any, Any]
self.gradient: NDArray[Any, Any]
self.hessian: NDArray[Any, Any]
self.ghost_atoms: List[Atom] = []
self.lp_atoms: List[Atom] = []
self.total_mass: int = 0
self.com: 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)

224
diceplayer/DPpack/Environment/System.py Normal file
View File

@@ -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
)
)

0
diceplayer/DPpack/Environment/__init__.py Normal file
View File

753
diceplayer/DPpack/External/Dice.py vendored Normal file
View File

@@ -0,0 +1,753 @@
import os
import shutil
import subprocess
import sys
from multiprocessing import Process, connection
from typing import Final, List, TextIO
import setproctitle
from diceplayer.DPpack.Utils.Misc import *
from diceplayer.DPpack.Utils.PTable import *
from diceplayer.DPpack.Utils.StepDTO import StepDTO
from diceplayer.DPpack.Utils.Validations import NotNull
from numpy import random
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"
nprocs: int = None
randominit = "first"
combrule = "*"
temp = 300.0
press = 1.0
isave = 1000
dens = None
ljname = None
outname = None
nmol: List[int] = None
nstep: List[int] = None
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.step.molecule)):
total_mass += self.step.molecule[i].total_mass * self.step.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.step.molecule[0].atom) * self.step.nmol[0]
for i in range(1, len(self.step.nmol)):
nsites += self.step.nmol[i] * len(self.step.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.step.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.step.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.step.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.step.initcyc:
fh.write("init = yesreadxyz\n")
fh.write("nstep = {}\n".format(self.step.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.step.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.step.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.step.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.step.initcyc:
fh.write("init = yesreadxyz\n")
fh.write("dens = {:<8.4f}\n".format(self.dens))
fh.write("vstep = {}\n".format(int(self.step.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.step.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.step.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.step.nmol)):
nsites_mm += self.step.nmol[i] * len(self.step.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.step.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.step.nmol)))
nsites_qm = (
len(self.step.molecule[0].atom)
+ len(self.step.molecule[0].ghost_atoms)
+ len(self.step.molecule[0].lp_atoms)
)
fh.write("{} {}\n".format(nsites_qm, self.step.molecule[0].molname))
for atom in self.step.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.step.molecule[0].atom[-1].lbl + 1
for i in self.step.molecule[0].ghost_atoms:
fh.write(
fstr.format(
ghost_label,
ghost_number,
self.step.molecule[0].atom[i].rx,
self.step.molecule[0].atom[i].ry,
self.step.molecule[0].atom[i].rz,
self.step.molecule[0].atom[i].chg,
0,
0,
)
)
ghost_label += 1
for lp in self.step.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.step.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.step.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.step.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, step: StepDTO):
self.step = step
def start(self, cycle: int) -> None:
procs = []
sentinels = []
for proc in range(1, self.step.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.step.nprocs + 1):
self.__print_last_config(cycle, proc)
def reset(self):
del self.step

654
diceplayer/DPpack/External/Gaussian.py vendored Normal file
View File

@@ -0,0 +1,654 @@
import os
import shutil
import subprocess
import sys
import textwrap
from typing import TextIO
import numpy as np
from diceplayer.DPpack.Environment.Atom import Atom
from diceplayer.DPpack.Environment.Molecule import Molecule
from diceplayer.DPpack.Utils.Misc import *
from diceplayer.DPpack.Utils.PTable import *
from diceplayer.DPpack.Utils.StepDTO import StepDTO
from diceplayer.DPpack.Utils.Validations import NotNull
class Gaussian:
mem = 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"
def __init__(self) -> None:
pass
@NotNull(requiredArgs=["qmprog","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.step.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.step.molecule[0].atom)):
fh.write(
" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format(
i + 1,
self.step.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.step.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.step.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.step.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.step.nprocs * self.step.ncores))
kword_line = "#P " + str(self.level)
if self.keywords != None:
kword_line += " " + self.keywords
if self.step.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.step.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.step.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, step: StepDTO):
self.step = step
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.step.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.step.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.step.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.step

0
diceplayer/DPpack/External/__init__.py vendored Normal file
View File

460
diceplayer/DPpack/MolHandling.py
View File

@@ -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))

348
diceplayer/DPpack/Molcas.py
View File

@@ -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

266
diceplayer/DPpack/Optimization.py
View File

@@ -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

804
diceplayer/DPpack/Player.py Normal file
View File

@@ -0,0 +1,804 @@
import os
import shutil
import sys
import textwrap
import types
from typing import TextIO
import yaml
from diceplayer.DPpack.Environment.Atom import Atom
from diceplayer.DPpack.Environment.Molecule import Molecule
from diceplayer.DPpack.Environment.System import System
from diceplayer.DPpack.External.Dice import Dice
from diceplayer.DPpack.External.Gaussian import Gaussian
from diceplayer.DPpack.Utils.Misc import *
from diceplayer.DPpack.Utils.PTable import *
from diceplayer.DPpack.Utils.StepDTO import StepDTO
from diceplayer.DPpack.Utils.Validations import NotNull
env = ["OMP_STACKSIZE"]
class Player:
maxcyc = None
opt = None
nprocs = None
qmprog = None
lps = None
ghosts = None
altsteps = None
combrule = None
TOL_RMS_FORCE = 3e-4
TOL_MAX_FORCE = 4.5e-4
TOL_RMS_STEP = 1.2e-3
TOL_MAX_SET = 1.8e-3
TRUST_RADIUS = None
continued: bool = False
def __init__(self, infile: TextIO, outfile: TextIO) -> None:
self.infile = infile
self.outfile = outfile
self.system = System()
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))
]
@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.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.ghosts == "yes" or self.lps == "yes"
):
sys.exit(
"Error: ghost atoms or lone pairs only available with 'pop = chelpg')"
)
# Check only if QM program is Molcas:
# if self.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.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.altsteps = max(min_steps, self.altsteps)
# altsteps value is always the nearest multiple of 1000
self.altsteps = round(self.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"
)
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.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.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.ghosts == "yes" or self.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(
StepDTO(
initcyc=self.initcyc,
nprocs=self.nprocs,
altsteps=self.altsteps,
nmol=self.system.nmols,
molecule=self.system.molecule,
)
)
self.dice.start(cycle)
self.dice.reset()
def gaussian_start(self, cycle: int, geomsfh: TextIO):
self.gaussian.configure(
self.initcyc,
self.nprocs,
self.dice.ncores,
self.altsteps,
self.switchcyc,
self.opt,
self.system.nmols,
self.system.molecule,
)
position = self.gaussian.start(cycle, self.outfile, self.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.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.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.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

1722
diceplayer/DPpack/SetGlobals.py
View File

File diff suppressed because it is too large Load Diff

0
diceplayer/DPpack/Misc.py → diceplayer/DPpack/Utils/Misc.py
View File

263
diceplayer/DPpack/Utils/Optimization.py Normal file
View File

@@ -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

0
diceplayer/DPpack/PTable.py → diceplayer/DPpack/Utils/PTable.py
View File

18
diceplayer/DPpack/Utils/StepDTO.py Normal file
View File

@@ -0,0 +1,18 @@
from dataclasses import dataclass
from typing import List
from diceplayer.DPpack.Environment.Molecule import Molecule
@dataclass
class StepDTO:
cycle: int = None
initcyc: int = None
nprocs: int = None
ncores: int = None
altsteps: int = None
switchcyc: int = None
opt: str = None
nmol: List[int] = None
molecule: List[Molecule] = None

15
diceplayer/DPpack/Utils/Validations.py Normal file
View File

@@ -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

0
diceplayer/DPpack/Utils/__init__.py Normal file
View File

366
diceplayer/__main__.py
View File

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

24
phb.ljc → phb.ljc.exemple
View File

@@ -1,18 +1,18 @@
*
2
12 BENZENO
1 6 0.000000 1.400000 0.000000 0.00000000 0.1100 3.7500
1 6 1.212436 0.700000 0.000000 0.00000000 0.1100 3.7500
1 6 1.212436 -0.700000 0.000000 0.00000000 0.1100 3.7500
1 6 0.000000 -1.400000 0.000000 0.00000000 0.1100 3.7500
1 6 -1.212436 -0.700000 0.000000 0.00000000 0.1100 3.7500
1 6 -1.212436 0.700000 0.000000 0.00000000 0.1100 3.7500
2 1 0.000000 2.488100 0.000000 0.00000000 0.0000 0.0000
2 1 2.154671 1.244000 0.000000 0.00000000 0.0000 0.0000
2 1 2.154671 -1.244000 0.000000 0.00000000 0.0000 0.0000
2 1 0.000000 -2.488100 0.000000 0.00000000 0.0000 0.0000
2 1 -2.154671 -1.244000 0.000000 0.00000000 0.0000 0.0000
2 1 -2.154671 1.244000 0.000000 0.00000000 0.0000 0.0000
1 6 1.407784 0.000000 0.000000 0.00000000 0.1100 3.7500
1 6 0.796966 1.274231 -0.338230 0.00000000 0.1100 3.7500
1 6 -0.703892 1.219173 -0.000000 0.00000000 0.1100 3.7500
1 6 -1.389932 -0.018988 0.227698 0.00000000 0.1100 3.7500
1 6 -0.703892 -1.219173 0.000000 0.00000000 0.1100 3.7500
1 6 0.703892 -1.219173 0.000000 0.00000000 0.1100 3.7500
2 1 2.501890 -0.000000 0.000000 0.00000000 0.0000 0.0000
2 1 1.415017 2.294961 -0.608843 0.00000000 0.0000 0.0000
2 1 -1.250943 2.166698 -0.000000 0.00000000 0.0000 0.0000
2 1 -2.501890 0.000000 0.000000 0.00000000 0.0000 0.0000
2 1 -0.519043 -2.235211 0.017991 0.00000000 0.0000 0.0000
2 1 1.250943 -2.166698 -0.000000 0.00000000 0.0000 0.0000
16 PLACEHOLDER
1 6 0.672026 -2.823446 0.002631 -0.11500000 0.0700 3.5500
1 6 2.072026 -2.823446 0.002631 -0.11500000 0.0700 3.5500

9
setup.py Normal file
View File

@@ -0,0 +1,9 @@
import PyInstaller.__main__
name = 'diceplayer'
PyInstaller.__main__.run([
'diceplayer/__main__.py',
'--onefile',
'-n{}'.format(name)
])