From b9b6c373b1780af14271a377bf9cf9a7ef885541 Mon Sep 17 00:00:00 2001 From: Vitor Hideyoshi Date: Thu, 23 Feb 2023 04:07:40 -0300 Subject: [PATCH] Initial Polarization Process --- config.example.yml | 2 +- crystalpol/gaussian.py | 44 +++++++++++++++++++++------- crystalpol/polarization.py | 36 +++++++++++++++++++---- crystalpol/shared/config/config.py | 1 + crystalpol/shared/system/crystal.py | 3 ++ crystalpol/shared/system/molecule.py | 10 +++---- 6 files changed, 73 insertions(+), 23 deletions(-) diff --git a/config.example.yml b/config.example.yml index 37a8f6a..c768993 100644 --- a/config.example.yml +++ b/config.example.yml @@ -1,5 +1,5 @@ crystal_pol: - mem: 28 + mem: 24 n_procs: 20 level: "b3lyp/aug-cc-pVDZ" pop: "chelpg" diff --git a/crystalpol/gaussian.py b/crystalpol/gaussian.py index 9cb10a2..e995e5e 100644 --- a/crystalpol/gaussian.py +++ b/crystalpol/gaussian.py @@ -3,10 +3,11 @@ from crystalpol.shared.system.crystal import Crystal from crystalpol.shared.config import Config from pathlib import Path, PosixPath -from typing import TextIO, Union +from typing import TextIO, Union, List import subprocess import textwrap import shutil +import sys import os @@ -23,9 +24,7 @@ class Gaussian: if self.config.pop not in ["chelpg", "mk", "nbo"]: self.config.pop = "chelpg" - def run(self, cycle: int, crystal: Crystal) -> None: - - self.create_simulation_dir() + def run(self, cycle: int, crystal: Crystal) -> List[float]: file = Path( "simfiles", @@ -51,11 +50,14 @@ class Gaussian: if exit_status != 0: raise RuntimeError("Gaussian process did not exit properly") - return self.read_charges_from_gaussian_output() + return self.read_charges_from_gaussian_output( + cycle, + crystal.get_number_of_charges() + ) def create_step_dir(self, cycle): step_dir = Path( - "simfiles", + self.config.simulation_dir, f"crystal-{str(cycle).zfill(2)}" ) if not os.path.exists(step_dir): @@ -79,6 +81,9 @@ class Gaussian: with open(file, 'w+') as fh: + chk_path = file.with_suffix('.chk') + fh.write(f"%Chk={chk_path}\n") + fh.write(f"%Mem={self.config.mem}Gb\n") fh.write(f"%Nprocs={self.config.n_procs}\n") @@ -123,13 +128,32 @@ class Gaussian: for atom in molecule: symbol = atom_symbol[atom.na] fh.write( - f"{symbol:<2s} " f"{float(atom.rx):>10.5f} " f"{float(atom.ry):>10.5f} " - f"{float(atom.rz):>10.5f}\n" + f"{float(atom.rz):>10.5f} " + f"{float(atom.chg):>10.5f}\n" ) fh.write("\n") - def read_charges_from_gaussian_output(self) -> None: - pass + def read_charges_from_gaussian_output(self, cycle, number_of_charges: int) -> List[float]: + step_dir = Path( + self.config.simulation_dir, + f"crystal-{str(cycle).zfill(2)}" + ) + filename = f"crystal-{str(cycle).zfill(2)}.log" + + file = Path(step_dir, filename) + try: + with open(file) as data: + lines = data.readlines() + except FileNotFoundError: + sys.exit("Error: cannot open file {}".format(file)) + + start = lines.pop(0).strip() + while start != "Fitting point charges to electrostatic potential": + start = lines.pop(0).strip() + + lines = lines[3:] # Consume 3 more lines + + return list(map(lambda x: float(x.split()[2]), lines[:number_of_charges])) \ No newline at end of file diff --git a/crystalpol/polarization.py b/crystalpol/polarization.py index 04564dc..85ab3b2 100644 --- a/crystalpol/polarization.py +++ b/crystalpol/polarization.py @@ -8,6 +8,8 @@ from crystalpol.gaussian import Gaussian from typing import List +import sys + class Polarization: __slots__ = ('geom_file', 'outfile', 'config', 'gaussian', 'crystal') @@ -23,9 +25,18 @@ class Polarization: def run(self): + self.gaussian.create_simulation_dir() + self.read_crystal() - self.gaussian.run(1, self.crystal) + cycle = 1 + charge_diff = sys.float_info.max + while charge_diff > self.config.charge_tolerance: + + charge_diff = self.update_crystal_charges( + self.gaussian.run(cycle, self.crystal), + ) + cycle += 1 def read_crystal(self) -> None: with open(self.geom_file, 'r') as geom_file: @@ -38,11 +49,24 @@ class Polarization: for molecule in molecules: self.crystal.add_cell([molecule]) + def update_crystal_charges(self, charges: List[float]) -> float: + + charge_diff = [] + + for cell in self.crystal: + for molecule in cell: + for index, atom in enumerate(molecule): + if atom.chg: + charge_diff.append(abs(atom.chg - charges[index])) + atom.chg = charges[index] + + return abs(max(charge_diff, key=abs)) if charge_diff else sys.float_info.max + def _get_molecules_from_lines(self, lines: List[str]) -> List[Molecule]: if (len(lines) % self.config.n_atoms) == 0: molecules: List[Molecule] = [] - for index, molecule in enumerate(split(lines, self.config.n_atoms)): + for index, molecule in enumerate(self.split(lines, self.config.n_atoms)): mol = Molecule(f"Molecule-{index}") for atom_line in molecule: symbol, rx, ry, rz = tuple(atom_line.split()) @@ -68,7 +92,7 @@ class Polarization: structure.append(atom.symbol) return structure - -def split(array: List, partitions: int): - for i in range(0, len(array), partitions): - yield array[i: i + partitions] + @staticmethod + def split(array: List, partitions: int): + for i in range(0, len(array), partitions): + yield array[i: i + partitions] diff --git a/crystalpol/shared/config/config.py b/crystalpol/shared/config/config.py index 781b923..0583c2d 100644 --- a/crystalpol/shared/config/config.py +++ b/crystalpol/shared/config/config.py @@ -10,6 +10,7 @@ class Config: n_procs: int = 1 pop: str = "chelpg" + charge_tolerance = 0.02 comment: str = "crystalpol" simulation_dir = "simfiles" mult: list = \ diff --git a/crystalpol/shared/system/crystal.py b/crystalpol/shared/system/crystal.py index 32c0c5f..56ada82 100644 --- a/crystalpol/shared/system/crystal.py +++ b/crystalpol/shared/system/crystal.py @@ -34,6 +34,9 @@ class Crystal: else: self.cells.append(cell) + def get_number_of_charges(self): + return len(self.cells[0][0]) + def _is_valid_cell(self, cell: List[Molecule]) -> bool: if len(cell) == len(self.structure): for i, molecule in enumerate(cell): diff --git a/crystalpol/shared/system/molecule.py b/crystalpol/shared/system/molecule.py index bc499ce..d5eaa66 100644 --- a/crystalpol/shared/system/molecule.py +++ b/crystalpol/shared/system/molecule.py @@ -28,12 +28,7 @@ class Molecule: __slots__ = ( 'mol_name', 'atoms', - 'position', - 'energy', - 'gradient', - 'hessian', - 'total_mass', - 'com' + 'position' ) def __init__(self, mol_name: str) -> None: @@ -51,6 +46,9 @@ class Molecule: for atom in self.atoms: yield atom + def __len__(self): + return len(self.atoms) + def add_atom(self, a: Atom) -> None: """ Adds Atom instance to the molecule.