Initial Polarization Process

config.example.yml

@@ -1,5 +1,5 @@
 crystal_pol:
-    mem: 28
+    mem: 24
     n_procs: 20
     level: "b3lyp/aug-cc-pVDZ"
     pop: "chelpg"

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

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]

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 = \

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

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.