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

Implements Refactoring in Player Class and Initial Working Version

Browse Source
This commit is contained in:
Vitor Hideyoshi
2023-06-02 20:20:38 -03:00
parent b440a0f05d
commit 33612f2d7b
21 changed files with 1193 additions and 983 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
diceplayer/shared/config/dice_dto.py → diceplayer/shared/config/dice_config.py
View File

@@ -6,7 +6,7 @@ from typing import List
@dataclass
class DiceDTO(Dataclass):
class DiceConfig(Dataclass):
"""
Data Transfer Object for the Dice configuration.
"""
@@ -53,4 +53,4 @@ class DiceDTO(Dataclass):
@classmethod
def from_dict(cls, param: dict):
return from_dict(DiceDTO, param)
return from_dict(DiceConfig, param)

3
diceplayer/shared/config/gaussian_dto.py → diceplayer/shared/config/gaussian_config.py
View File

@@ -11,10 +11,11 @@ class GaussianDTO(Dataclass):
"""
level: str
qmprog: str
keywords: str
chgmult = [0, 1]
pop: str = 'chelpg'
chg_tol: float = 0.01
keywords: str = None
def __post_init__(self):
if self.qmprog not in ("g03", "g09", "g16"):

47
diceplayer/shared/config/player_config.py Normal file
View File

@@ -0,0 +1,47 @@
from diceplayer.shared.utils.dataclass_protocol import Dataclass
from diceplayer.shared.config.gaussian_config import GaussianDTO
from diceplayer.shared.config.dice_config import DiceConfig
from dataclasses import dataclass
from dacite import from_dict
@dataclass
class PlayerConfig(Dataclass):
"""
Data Transfer Object for the player configuration.
"""
opt: bool
maxcyc: int
nprocs: int
ncores: int
dice: DiceConfig
gaussian: GaussianDTO
mem: int = None
switchcyc: int = 3
qmprog: str = 'g16'
altsteps: int = 20000
simulation_dir = 'simfiles'
def __post_init__(self):
MIN_STEP = 20000
# altsteps value is always the nearest multiple of 1000
self.altsteps = round(max(MIN_STEP, self.altsteps) / 1000) * 1000
@classmethod
def from_dict(cls, param: dict):
if param['dice'] is None:
raise ValueError(
"Error: 'dice' keyword not specified in config file."
)
param['dice'] = DiceConfig.from_dict(param['dice'])
if param['gaussian'] is None:
raise ValueError(
"Error: 'gaussian' keyword not specified in config file."
)
param['gaussian'] = GaussianDTO.from_dict(param['gaussian'])
return from_dict(PlayerConfig, param)

28
diceplayer/shared/config/player_dto.py
View File

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

22
diceplayer/shared/config/step_dto.py
View File

@@ -1,22 +0,0 @@
from diceplayer.shared.environment.molecule import Molecule
from diceplayer.shared.config.player_dto import PlayerDTO
from dataclasses import dataclass
from typing import List
@dataclass
class StepDTO:
"""
Data Transfer Object for the step configuration.
"""
ncores: int
nprocs: int
simulation_dir: str
altsteps: int
nmol: List[int] = None
molecule: List[Molecule] = None
charges: List[float] = None
position: List[float] = None

64
diceplayer/shared/environment/molecule.py
View File

@@ -1,15 +1,18 @@
import logging
import math
from copy import deepcopy
from typing import List, Any, Tuple, Final, Union
import numpy as np
from nptyping import NDArray, Shape, Float
from numpy.linalg import linalg
from __future__ import annotations
from diceplayer.shared.utils.ptable import ghost_number
from diceplayer.shared.environment.atom import Atom
from diceplayer.shared.utils.misc import BOHR2ANG
from diceplayer.shared.utils.ptable import ghost_number
from diceplayer import logger
from nptyping import NDArray, Shape, Float
from numpy.linalg import linalg
import numpy as np
from typing import List, Any, Tuple, Union
from copy import deepcopy
import logging
import math
class Molecule:
@@ -185,11 +188,18 @@ class Molecule:
return position
def update_charges(self, charges: List[float]) -> None:
def update_charges(self, charges: NDArray) -> int:
"""
Updates the charges of the atoms of the molecule and
returns the max difference between the new and old charges
"""
diff = 0
for i, atom in enumerate(self.atom):
diff = max(diff, abs(atom.chg - charges[i]))
atom.chg = charges[i]
return diff
# @staticmethod
# def update_hessian(
# step: np.ndarray,
@@ -299,48 +309,48 @@ class Molecule:
Prints the Molecule information into a Output File
"""
logging.info(
" Center of mass = ( {:>10.4f} , {:>10.4f} , {:>10.4f} )\n".format(
logger.info(
" Center of mass = ( {:>10.4f} , {:>10.4f} , {:>10.4f} )".format(
self.com[0], self.com[1], self.com[2]
)
)
inertia = self.inertia_tensor()
evals, evecs = self.principal_axes()
logging.info(
" Moments of inertia = {:>9E} {:>9E} {:>9E}\n".format(
logger.info(
" Moments of inertia = {:>9E} {:>9E} {:>9E}".format(
evals[0], evals[1], evals[2]
)
)
logging.info(
" Major principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
logger.info(
" Major principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )".format(
evecs[0, 0], evecs[1, 0], evecs[2, 0]
)
)
logging.info(
" Inter principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
logger.info(
" Inter principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )".format(
evecs[0, 1], evecs[1, 1], evecs[2, 1]
)
)
logging.info(
" Minor principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
logger.info(
" Minor principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )".format(
evecs[0, 2], evecs[1, 2], evecs[2, 2]
)
)
sizes = self.sizes_of_molecule()
logging.info(
" Characteristic lengths = ( {:>6.2f} , {:>6.2f} , {:>6.2f} )\n".format(
logger.info(
" Characteristic lengths = ( {:>6.2f} , {:>6.2f} , {:>6.2f} )".format(
sizes[0], sizes[1], sizes[2]
)
)
logging.info(" Total mass = {:>8.2f} au\n".format(self.total_mass))
logger.info(" Total mass = {:>8.2f} au".format(self.total_mass))
chg_dip = self.charges_and_dipole()
logging.info(" Total charge = {:>8.4f} e\n".format(chg_dip[0]))
logging.info(
" Dipole moment = ( {:>9.4f} , {:>9.4f} , {:>9.4f} ) Total = {:>9.4f} Debye\n\n".format(
logger.info(" Total charge = {:>8.4f} e".format(chg_dip[0]))
logger.info(
" Dipole moment = ( {:>9.4f} , {:>9.4f} , {:>9.4f} ) Total = {:>9.4f} Debye".format(
chg_dip[1], chg_dip[2], chg_dip[3], chg_dip[4]
)
)

92
diceplayer/shared/environment/system.py
View File

@@ -1,6 +1,7 @@
from diceplayer.shared.environment.molecule import Molecule
from diceplayer.shared.utils.ptable import atomsymb
from diceplayer.shared.utils.misc import BOHR2ANG
from diceplayer import logger
from typing import List, Tuple, TextIO
from copy import deepcopy
@@ -11,50 +12,42 @@ import math
class System:
"""
System class declaration. This class is used throughout the DicePlayer program to represent the system containing the molecules.
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
"""
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:
Initializes an empty system object that will be populated afterwards
"""
Adds a new molecule type to the system
self.nmols: List[int] = []
self.molecule: List[Molecule] = []
Args:
nmols (int): Number of molecules of the new type in the system
m (Molecule): The instance of the new type of molecule
"""
def add_type(self, m: Molecule) -> None:
"""
Adds a new molecule type to the system
Args:
m (Molecule): The instance of the new type of molecule
"""
if isinstance(m, Molecule) is False:
raise TypeError("Error: molecule is not a Molecule instance")
self.molecule.append(m)
if isinstance(nmols, int) is False:
raise TypeError("Error: nmols is not an integer")
self.nmols.append(nmols)
def update_molecule(self, position: np.ndarray, fh: TextIO) -> None:
def update_molecule(self, position: np.ndarray) -> None:
"""Updates the position of the molecule in the Output file
Args:
position (np.ndarray): numpy position vector
fh (TextIO): Output file
"""
Args:
position (np.ndarray): numpy position vector
"""
position_in_ang = (position * BOHR2ANG).tolist()
self.add_type(self.nmols[0], deepcopy(self.molecule[0]))
self.add_type(deepcopy(self.molecule[0]))
for atom in self.molecule[-1].atom:
atom.rx = position_in_ang.pop(0)
atom.ry = position_in_ang.pop(0)
atom.rz = position_in_ang.pop(0)
@@ -62,8 +55,8 @@ class System:
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))
logger.info("Projected new conformation of reference molecule with RMSD fit")
logger.info(f"RMSD = {rmsd:>8.5f} Angstrom")
def rmsd_fit(self, p_index: int, r_index: int) -> Tuple[float, Molecule]:
@@ -200,7 +193,6 @@ class System:
#
# return min_dist, nearestmol
# def print_geom(self, cycle: int, fh: TextIO) -> None:
# """
# Print the geometry of the molecule in the Output file
@@ -220,22 +212,22 @@ class System:
# )
# )
#
# def printChargesAndDipole(self, cycle: int, fh: TextIO) -> None:
# """
# Print the charges and dipole of the molecule in the Output file
#
# Args:
# cycle (int): Number of the cycle
# fh (TextIO): Output file
# """
#
# fh.write("Cycle # {}\n".format(cycle))
# fh.write("Number of site: {}\n".format(len(self.molecule[0].atom)))
#
# chargesAndDipole = self.molecule[0].charges_and_dipole()
#
# fh.write(
# "{:>10.6f} {:>10.6f} {:>10.6f} {:>10.6f} {:>10.6f}\n".format(
# chargesAndDipole[0], chargesAndDipole[1], chargesAndDipole[2], chargesAndDipole[3], chargesAndDipole[4]
# )
# )
def print_charges_and_dipole(self, cycle: int) -> None:
"""
Print the charges and dipole of the molecule in the Output file
Args:
cycle (int): Number of the cycle
fh (TextIO): Output file
"""
logger.info("Cycle # {}\n".format(cycle))
logger.info("Number of site: {}\n".format(len(self.molecule[0].atom)))
chargesAndDipole = self.molecule[0].charges_and_dipole()
logger.info(
"{:>10.6f} {:>10.6f} {:>10.6f} {:>10.6f} {:>10.6f}\n".format(
chargesAndDipole[0], chargesAndDipole[1], chargesAndDipole[2], chargesAndDipole[3], chargesAndDipole[4]
)
)

17
diceplayer/shared/interface/__interface.py
View File

@@ -1,20 +1,23 @@
from diceplayer.shared.utils.dataclass_protocol import Dataclass
from __future__ import annotations
from diceplayer.shared.config.player_config import PlayerConfig
from diceplayer.shared.environment.system import System
from abc import ABC, abstractmethod
class Interface(ABC):
__slots__ = [
'config'
'step',
'system'
]
@abstractmethod
def __init__(self, data: dict):
pass
def __init__(self):
self.system: System | None = None
self.step: PlayerConfig | None = None
@staticmethod
@abstractmethod
def set_config(data: dict) -> Dataclass:
def configure(self, step: PlayerConfig, system: System):
pass
@abstractmethod

134
diceplayer/shared/interface/dice_interface.py
View File

@@ -1,7 +1,7 @@
from __future__ import annotations
from diceplayer.shared.config.dice_dto import DiceDTO
from diceplayer.shared.config.step_dto import StepDTO
from diceplayer.shared.config.player_config import PlayerConfig
from diceplayer.shared.environment.system import System
from diceplayer.shared.interface import Interface
from diceplayer import logger
@@ -26,23 +26,14 @@ MAX_SEED: Final[int] = 4294967295
class DiceInterface(Interface):
title = "Diceplayer run"
def __init__(self, data: dict):
self.config: DiceDTO = self.set_config(data)
self.step: StepDTO | None = None
@staticmethod
def set_config(data: dict) -> DiceDTO:
return DiceDTO.from_dict(data)
def configure(self, step: any):
def configure(self, step: PlayerConfig, system: System):
self.step = step
self.system = system
def start(self, cycle: int):
procs = []
sentinels = []
logger.info(f"---------------------- DICE - CYCLE {cycle} --------------------------\n")
for proc in range(1, self.step.nprocs + 1):
p = Process(target=self._simulation_process, args=(cycle, proc))
p.start()
@@ -66,6 +57,7 @@ class DiceInterface(Interface):
def reset(self):
del self.step
del self.system
def _simulation_process(self, cycle: int, proc: int):
setproctitle(f"diceplayer-step{cycle:0d}-p{proc:0d}")
@@ -102,7 +94,7 @@ class DiceInterface(Interface):
# This is logic is used to make the initial configuration file
# for the next cycle using the last.xyz file from the previous cycle.
if self.config.randominit == 'first' and cycle > 1:
if self.step.dice.randominit == 'first' and cycle > 1:
last_xyz = Path(
self.step.simulation_dir,
f"step{(cycle - 1):02d}",
@@ -115,15 +107,15 @@ class DiceInterface(Interface):
with open(last_xyz, 'r') as last_xyz_file:
self._make_init_file(proc_dir, last_xyz_file)
last_xyz_file.seek(0)
self.config.dens = self._new_density(last_xyz_file)
self.step.dice.dens = self._new_density(last_xyz_file)
else:
self._make_nvt_ter(cycle, proc_dir)
if len(self.config.nstep) == 2:
self._make_nvt_eq(proc_dir)
if len(self.step.dice.nstep) == 2:
self._make_nvt_eq(cycle, proc_dir)
elif len(self.config.nstep) == 3:
elif len(self.step.dice.nstep) == 3:
self._make_npt_ter(cycle, proc_dir)
self._make_npt_eq(proc_dir)
@@ -142,13 +134,13 @@ class DiceInterface(Interface):
os.chdir(proc_dir)
if not (self.config.randominit == 'first' and cycle > 1):
if not (self.step.dice.randominit == 'first' and cycle > 1):
self.run_dice_file(cycle, proc, "NVT.ter")
if len(self.config.nstep) == 2:
if len(self.step.dice.nstep) == 2:
self.run_dice_file(cycle, proc, "NVT.eq")
elif len(self.config.nstep) == 3:
elif len(self.step.dice.nstep) == 3:
self.run_dice_file(cycle, proc, "NPT.ter")
self.run_dice_file(cycle, proc, "NPT.eq")
@@ -175,15 +167,15 @@ class DiceInterface(Interface):
xyz_lines = last_xyz_file.readlines()
nsites_mm = 0
for i in range(1, len(self.step.nmol)):
nsites_mm += self.step.nmol[i] * len(self.step.molecule[i].atom)
for i in range(1, len(self.step.dice.nmol)):
nsites_mm += self.step.dice.nmol[i] * len(self.system.molecule[i].atom)
xyz_lines = xyz_lines[-nsites_mm:]
input_file = Path(proc_dir, self.config.outname + ".xy")
input_file = Path(proc_dir, self.step.dice.outname + ".xy")
with open(input_file, 'w') as f:
for atom in self.step.molecule[0].atom:
for atom in self.system.molecule[0].atom:
f.write(
f"{atom.rx:>10.6f} {atom.ry:>10.6f} {atom.rz:>10.6f}\n"
)
@@ -204,8 +196,8 @@ class DiceInterface(Interface):
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]
for i in range(len(self.system.molecule)):
total_mass += self.system.molecule[i].total_mass * self.step.dice.nmol[i]
density = (total_mass / volume) * UMAANG3_TO_GCM3
@@ -216,21 +208,21 @@ class DiceInterface(Interface):
with open(file, 'w') as f:
f.write(f"title = {self.title} - NVT Thermalization\n")
f.write(f"ncores = {self.step.ncores}\n")
f.write(f"ljname = {self.config.ljname}\n")
f.write(f"outname = {self.config.outname}\n")
f.write(f"ljname = {self.step.dice.ljname}\n")
f.write(f"outname = {self.step.dice.outname}\n")
mol_string = " ".join(str(x) for x in self.config.nmol)
mol_string = " ".join(str(x) for x in self.step.dice.nmol)
f.write(f"nmol = {mol_string}\n")
f.write(f"dens = {self.config.dens}\n")
f.write(f"temp = {self.config.temp}\n")
f.write(f"dens = {self.step.dice.dens}\n")
f.write(f"temp = {self.step.dice.temp}\n")
if self.config.randominit == "first" and cycle > 1:
if self.step.dice.randominit == "first" and cycle > 1:
f.write(f"init = yesreadxyz\n")
f.write(f"nstep = {self.step.altsteps}\n")
else:
f.write(f"init = yes\n")
f.write(f"nstep = {self.config.nstep[0]}\n")
f.write(f"nstep = {self.step.dice.nstep[0]}\n")
f.write("vstep = 0\n")
f.write("mstop = 1\n")
@@ -241,30 +233,36 @@ class DiceInterface(Interface):
seed = int(1e6 * random.random())
f.write(f"seed = {seed}\n")
f.write(f"upbuf = {self.config.upbuf}")
f.write(f"upbuf = {self.step.dice.upbuf}")
def _make_nvt_eq(self, proc_dir):
def _make_nvt_eq(self, cycle, proc_dir):
file = Path(proc_dir, "NVT.eq")
with open(file, 'w') as f:
f.write(f"title = {self.title} - NVT Production\n")
f.write(f"ncores = {self.step.ncores}\n")
f.write(f"ljname = {self.config.ljname}\n")
f.write(f"outname = {self.config.outname}\n")
f.write(f"ljname = {self.step.dice.ljname}\n")
f.write(f"outname = {self.step.dice.outname}\n")
mol_string = " ".join(str(x) for x in self.config.nmol)
mol_string = " ".join(str(x) for x in self.step.dice.nmol)
f.write(f"nmol = {mol_string}\n")
f.write(f"dens = {self.config.dens}\n")
f.write(f"temp = {self.config.temp}\n")
f.write("init = no\n")
f.write(f"nstep = {self.config.nstep[1]}\n")
f.write(f"dens = {self.step.dice.dens}\n")
f.write(f"temp = {self.step.dice.temp}\n")
if self.step.dice.randominit == "first" and cycle > 1:
f.write("init = yesreadxyz\n")
else:
f.write("init = no\n")
f.write(f"nstep = {self.step.dice.nstep[1]}\n")
f.write("vstep = 0\n")
f.write("mstop = 1\n")
f.write("accum = no\n")
f.write("iprint = 1\n")
f.write(f"isave = {self.config.isave}\n")
f.write(f"isave = {self.step.dice.isave}\n")
f.write(f"irdf = {10 * self.step.nprocs}\n")
seed = int(1e6 * random.random())
@@ -276,22 +274,22 @@ class DiceInterface(Interface):
with open(file, 'w') as f:
f.write(f"title = {self.title} - NPT Thermalization\n")
f.write(f"ncores = {self.step.ncores}\n")
f.write(f"ljname = {self.config.ljname}\n")
f.write(f"outname = {self.config.outname}\n")
f.write(f"ljname = {self.step.dice.ljname}\n")
f.write(f"outname = {self.step.dice.outname}\n")
mol_string = " ".join(str(x) for x in self.config.nmol)
mol_string = " ".join(str(x) for x in self.step.dice.nmol)
f.write(f"nmol = {mol_string}\n")
f.write(f"press = {self.config.press}\n")
f.write(f"temp = {self.config.temp}\n")
f.write(f"press = {self.step.dice.press}\n")
f.write(f"temp = {self.step.dice.temp}\n")
if self.config.randominit == "first" and cycle > 1:
if self.step.dice.randominit == "first" and cycle > 1:
f.write("init = yesreadxyz\n")
f.write(f"dens = {self.config.dens:<8.4f}\n")
f.write(f"dens = {self.step.dice.dens:<8.4f}\n")
f.write(f"vstep = {int(self.step.altsteps / 5)}\n")
else:
f.write("init = no\n")
f.write(f"vstep = {int(self.config.nstep[1] / 5)}\n")
f.write(f"vstep = {int(self.step.dice.nstep[1] / 5)}\n")
f.write("nstep = 5\n")
f.write("mstop = 1\n")
@@ -308,23 +306,23 @@ class DiceInterface(Interface):
with open(file, 'w') as f:
f.write(f"title = {self.title} - NPT Production\n")
f.write(f"ncores = {self.step.ncores}\n")
f.write(f"ljname = {self.config.ljname}\n")
f.write(f"outname = {self.config.outname}\n")
f.write(f"ljname = {self.step.dice.ljname}\n")
f.write(f"outname = {self.step.dice.outname}\n")
mol_string = " ".join(str(x) for x in self.config.nmol)
mol_string = " ".join(str(x) for x in self.step.dice.nmol)
f.write(f"nmol = {mol_string}\n")
f.write(f"press = {self.config.press}\n")
f.write(f"temp = {self.config.temp}\n")
f.write(f"press = {self.step.dice.press}\n")
f.write(f"temp = {self.step.dice.temp}\n")
f.write(f"nstep = 5\n")
f.write(f"vstep = {int(self.config.nstep[2] / 5)}\n")
f.write(f"vstep = {int(self.step.dice.nstep[2] / 5)}\n")
f.write("init = no\n")
f.write("mstop = 1\n")
f.write("accum = no\n")
f.write("iprint = 1\n")
f.write(f"isave = {self.config.isave}\n")
f.write(f"isave = {self.step.dice.isave}\n")
f.write(f"irdf = {10 * self.step.nprocs}\n")
seed = int(1e6 * random.random())
@@ -333,15 +331,15 @@ class DiceInterface(Interface):
def _make_potentials(self, proc_dir):
fstr = "{:<3d} {:>3d} {:>10.5f} {:>10.5f} {:>10.5f} {:>10.6f} {:>9.5f} {:>7.4f}\n"
file = Path(proc_dir, self.config.ljname)
file = Path(proc_dir, self.step.dice.ljname)
with open(file, 'w') as f:
f.write(f"{self.config.combrule}\n")
f.write(f"{len(self.step.nmol)}\n")
f.write(f"{self.step.dice.combrule}\n")
f.write(f"{len(self.step.dice.nmol)}\n")
nsites_qm = len(self.step.molecule[0].atom)
f.write(f"{nsites_qm} {self.step.molecule[0].molname}\n")
nsites_qm = len(self.system.molecule[0].atom)
f.write(f"{nsites_qm} {self.system.molecule[0].molname}\n")
for atom in self.step.molecule[0].atom:
for atom in self.system.molecule[0].atom:
f.write(
fstr.format(
atom.lbl,
@@ -355,7 +353,7 @@ class DiceInterface(Interface):
)
)
for mol in self.step.molecule[1:]:
for mol in self.system.molecule[1:]:
f.write(f"{len(mol.atom)} {mol.molname}\n")
for atom in mol.atom:
f.write(
@@ -378,12 +376,12 @@ class DiceInterface(Interface):
[
"bash",
"-c",
f"exec -a dice-step{cycle}-p{proc} {self.config.progname} < {infile.name} > {outfile.name}",
f"exec -a dice-step{cycle}-p{proc} {self.step.dice.progname} < {infile.name} > {outfile.name}",
]
)
else:
exit_status = subprocess.call(
self.config.progname, stdin=infile, stdout=outfile
self.step.dice.progname, stdin=infile, stdout=outfile
)
if exit_status != 0:

380
diceplayer/shared/interface/gaussian_interface.py
View File

@@ -1,22 +1,380 @@
from diceplayer.shared.config.gaussian_dto import GaussianDTO
from __future__ import annotations
from diceplayer.shared.config.player_config import PlayerConfig
from diceplayer.shared.environment.molecule import Molecule
from diceplayer.shared.environment.system import System
from diceplayer.shared.environment.atom import Atom
from diceplayer.shared.utils.misc import date_time
from diceplayer.shared.utils.ptable import atomsymb
from diceplayer.shared.interface import Interface
from diceplayer import logger
from typing import Tuple, List, Dict, Any
from nptyping import NDArray
import numpy as np
from pathlib import Path
import subprocess
import textwrap
import shutil
import os
class GaussianInterface(Interface):
def configure(self, step_dto: PlayerConfig, system: System):
self.system = system
self.step = step_dto
def __init__(self, data: dict):
self.config: GaussianDTO = self.set_config(data)
def start(self, cycle: int) -> Dict[str, NDArray]:
self._make_qm_dir(cycle)
@staticmethod
def set_config(data: dict) -> GaussianDTO:
return GaussianDTO.from_dict(data)
if cycle > 1:
self._copy_chk_file_from_previous_step(cycle)
def configure(self):
pass
asec_charges = self.populate_asec_vdw(cycle)
self._make_gaussian_input_file(
cycle,
asec_charges
)
def start(self, cycle: int):
pass
self._run_gaussian(cycle)
self._run_formchk(cycle)
return_value = {}
if self.step.opt:
# return_value['position'] = np.array(
# self._run_optimization(cycle)
# )
raise NotImplementedError("Optimization not implemented yet.")
else:
return_value['charges'] = np.array(
self._read_charges_from_fchk(cycle)
)
return return_value
def reset(self):
pass
del self.step
del self.system
def _make_qm_dir(self, cycle: int):
qm_dir_path = Path(
self.step.simulation_dir,
f"step{cycle:02d}",
"qm"
)
if not qm_dir_path.exists():
qm_dir_path.mkdir()
def _copy_chk_file_from_previous_step(self, cycle: int):
current_chk_file_path = Path(
self.step.simulation_dir,
f"step{cycle:02d}",
"qm",
f"asec.chk"
)
if current_chk_file_path.exists():
raise FileExistsError(
f"File {current_chk_file_path} already exists."
)
previous_chk_file_path = Path(
self.step.simulation_dir,
f"step{(cycle - 1):02d}",
"qm",
f"asec.chk"
)
if not previous_chk_file_path.exists():
raise FileNotFoundError(
f"File {previous_chk_file_path} does not exist."
)
shutil.copy(previous_chk_file_path, current_chk_file_path)
def populate_asec_vdw(self, cycle: int) -> list[dict]:
norm_factor = self._calculate_norm_factor()
nsitesref = len(self.system.molecule[0].atom)
nsites_total = self._calculate_total_number_of_sites(nsitesref)
proc_charges = []
for proc in range(1, self.step.nprocs + 1):
proc_charges.append(self._read_charges_from_last_step(cycle, proc))
asec_charges, thickness, picked_mols = \
self._evaluate_proc_charges(nsites_total, proc_charges)
logger.info(f"In average, {(sum(picked_mols) / norm_factor):^7.2f} molecules\n"
f"were selected from each of the {len(picked_mols)} configurations\n"
f"of the production simulations to form the ASEC, comprising a shell with\n"
f"minimum thickness of {(sum(thickness) / norm_factor):>6.2f} Angstrom\n"
)
for charge in asec_charges:
charge['chg'] = charge['chg'] / norm_factor
return asec_charges
def _calculate_norm_factor(self) -> int:
if self.step.dice.nstep[-1] % self.step.dice.isave == 0:
nconfigs = round(self.step.dice.nstep[-1] / self.step.dice.isave)
else:
nconfigs = int(self.step.dice.nstep[-1] / self.step.dice.isave)
return nconfigs * self.step.nprocs
def _calculate_total_number_of_sites(self, nsitesref) -> int:
nsites_total = self.step.dice.nmol[0] * nsitesref
for i in range(1, len(self.step.dice.nmol)):
nsites_total += self.step.dice.nmol[i] * len(self.system.molecule[i].atom)
return nsites_total
def _read_charges_from_last_step(self, cycle: int, proc: int) -> list[str]:
last_xyz_file_path = Path(
self.step.simulation_dir,
f"step{cycle:02d}",
f"p{proc:02d}",
"last.xyz"
)
if not last_xyz_file_path.exists():
raise FileNotFoundError(
f"File {last_xyz_file_path} does not exist."
)
with open(last_xyz_file_path, 'r') as last_xyz_file:
lines = last_xyz_file.readlines()
return lines
def _evaluate_proc_charges(self, total_nsites: int, proc_charges: list[list[str]]) -> Tuple[
List[Dict[str, float | Any]], List[float], List[int]]:
asec_charges = []
thickness = []
picked_mols = []
for charges in proc_charges:
charges_nsites = int(charges.pop(0))
if int(charges_nsites) != total_nsites:
raise ValueError(
f"Number of sites does not match total number of sites."
)
thickness.append(
self._calculate_proc_thickness(charges)
)
nsites_ref_mol = len(self.system.molecule[0].atom)
charges = charges[nsites_ref_mol:]
mol_count = 0
for type in range(len(self.step.dice.nmol)):
if type == 0:
# Reference Molecule must be ignored from type 0
nmols = self.step.dice.nmol[type] - 1
else:
nmols = self.step.dice.nmol[type]
for mol in range(nmols):
new_molecule = Molecule("ASEC TMP MOLECULE")
for site in range(len(self.system.molecule[type].atom)):
line = charges.pop(0).split()
if line[0].title() != atomsymb[self.system.molecule[type].atom[site].na].strip():
raise SyntaxError(
f"Error: Invalid Dice Output. Atom type does not match."
)
new_molecule.add_atom(
Atom(
self.system.molecule[type].atom[site].lbl,
self.system.molecule[type].atom[site].na,
float(line[1]),
float(line[2]),
float(line[3]),
self.system.molecule[type].atom[site].chg,
self.system.molecule[type].atom[site].eps,
self.system.molecule[type].atom[site].sig,
)
)
distance = self.system.molecule[0] \
.minimum_distance(new_molecule)
if distance < thickness[-1]:
for atom in new_molecule.atom:
asec_charges.append(
{"lbl": atomsymb[atom.na], "rx": atom.rx, "ry": atom.ry, "rz": atom.rz, "chg": atom.chg}
)
mol_count += 1
picked_mols.append(mol_count)
return asec_charges, thickness, picked_mols
def _calculate_proc_thickness(self, charges: list[str]) -> float:
box = charges.pop(0).split()[-3:]
box = [float(box[0]), float(box[1]), float(box[2])]
sizes = self.system.molecule[0].sizes_of_molecule()
return min(
[
(box[0] - sizes[0]) / 2,
(box[1] - sizes[1]) / 2,
(box[2] - sizes[2]) / 2,
]
)
def _make_gaussian_input_file(self, cycle: int, asec_charges: list[dict]) -> None:
gaussian_input_file_path = Path(
self.step.simulation_dir,
f"step{cycle:02d}",
"qm",
f"asec.gjf"
)
with open(gaussian_input_file_path, 'w') as gaussian_input_file:
gaussian_input_file.writelines(
self._generate_gaussian_input(cycle, asec_charges)
)
def _generate_gaussian_input(self, cycle: int, asec_charges: list[dict]) -> list[str]:
gaussian_input = ["%Chk=asec.chk\n"]
if self.step.mem is not None:
gaussian_input.append(f"%Mem={self.step.mem}GB\n")
gaussian_input.append(f"%Nprocs={self.step.nprocs * self.step.ncores}\n")
kwords_line = f"#P {self.step.gaussian.level}"
if self.step.gaussian.keywords:
kwords_line += " " + self.step.gaussian.keywords
if self.step.opt == "yes":
kwords_line += " Force"
kwords_line += " NoSymm"
kwords_line += f" Pop={self.step.gaussian.pop} Density=Current"
if cycle > 1:
kwords_line += " Guess=Read"
gaussian_input.append(textwrap.fill(kwords_line, 90))
gaussian_input.append("\n")
gaussian_input.append("\nForce calculation - Cycle number {}\n".format(cycle))
gaussian_input.append("\n")
gaussian_input.append(f"{self.step.gaussian.chgmult[0]},{self.step.gaussian.chgmult[1]}\n")
for atom in self.system.molecule[0].atom:
symbol = atomsymb[atom.na]
gaussian_input.append(
"{:<2s} {:>10.5f} {:>10.5f} {:>10.5f}\n".format(
symbol, atom.rx, atom.ry, atom.rz
)
)
gaussian_input.append("\n")
for charge in asec_charges:
gaussian_input.append(
"{:>10.5f} {:>10.5f} {:>10.5f} {:>11.8f}\n".format(
charge['rx'], charge['ry'], charge['rz'], charge['chg']
)
)
gaussian_input.append("\n")
return gaussian_input
def _run_gaussian(self, cycle: int) -> None:
qm_dir = Path(
self.step.simulation_dir,
f"step{(cycle):02d}",
"qm"
)
working_dir = os.getcwd()
os.chdir(qm_dir)
infile = "asec.gjf"
operation = None
if self.step.opt:
operation = "forces"
else:
operation = "charges"
logger.info(
f"Calculation of {operation} initiated with Gaussian on {date_time()}\n"
)
if shutil.which("bash") is not None:
exit_status = subprocess.call(
[
"bash",
"-c",
"exec -a {}-step{} {} {}".format(
self.step.gaussian.qmprog, cycle, self.step.gaussian.qmprog, infile
),
]
)
else:
exit_status = subprocess.call([self.step.gaussian.qmprog, infile])
if exit_status != 0:
raise SystemError("Gaussian process did not exit properly")
logger.info(f"Calculation of {operation} finished on {date_time()}")
os.chdir(working_dir)
def _run_formchk(self, cycle: int):
qm_dir = Path(
self.step.simulation_dir,
f"step{(cycle):02d}",
"qm"
)
work_dir = os.getcwd()
os.chdir(qm_dir)
logger.info("Formatting the checkpoint file... \n")
exit_status = subprocess.call(["formchk", "asec.chk"], stdout=subprocess.DEVNULL)
if exit_status != 0:
raise SystemError("Formchk process did not exit properly")
logger.info("Done\n")
os.chdir(work_dir)
def _read_charges_from_fchk(self, cycle: int):
fchk_file_path = Path(
"simfiles",
f"step{cycle:02d}",
"qm",
"asec.fchk"
)
with open(fchk_file_path) as fchk:
fchkfile = fchk.readlines()
if self.step.gaussian.pop in ["chelpg", "mk"]:
CHARGE_FLAG = "ESP Charges"
else:
CHARGE_FLAG = "ESP Charges"
start = fchkfile.pop(0).strip()
while start.find(CHARGE_FLAG) != 0: # expression in begining of line
start = fchkfile.pop(0).strip()
charges: List[float] = []
while len(charges) < len(self.system.molecule[0].atom):
charges.extend([float(x) for x in fchkfile.pop(0).split()])
return charges