add support for point charges

cctk/__init__.py

@@ -14,3 +14,5 @@
 from .pdb_file import PDBFile
 
 from .si_file import SIFile
+
+from .point_charge import PointCharge

cctk/gaussian_file.py

@@ -5,6 +5,7 @@
 
 from cctk import File, Molecule, ConformationalEnsemble, OneIndexedArray
 from cctk.helper_functions import get_symbol, get_number, get_corrected_free_energy
+import cctk
 
 import cctk.parse_gaussian as parse
 
@@ -138,7 +139,7 @@ def __str__(self):
         return f"GaussianFile (title=\"{str(self.title)}\", {len(self.ensemble)} entries in Ensemble)"
 
     @classmethod
-    def write_molecule_to_file(cls, filename, molecule, route_card, link0={"mem": "32GB", "nprocshared": 16}, footer=None, title="title", append=False, print_symbol=False):
+    def write_molecule_to_file(cls, filename, molecule, route_card, link0={"mem": "32GB", "nprocshared": 16}, footer=None, title="title", append=False, print_symbol=False, point_charges=None):
         """
         Write a ``.gjf`` file using the given molecule.
 
@@ -161,6 +162,11 @@ def write_molecule_to_file(cls, filename, molecule, route_card, link0={"mem": "3
         if not re.match(r"^#p", route_card):
             warnings.warn(f"route card doesn't start with #p: {route_card}")
 
+        if point_charges is not None:
+            assert isinstance(point_charges, list), "point_charges must be list"
+            assert all([isinstance(pc, cctk.PointCharge) for pc in point_charges]), "point_charges must be list of point charges"
+            assert re.search(r"charge", route_card, flags=re.IGNORECASE), "charge must be in route_card if point_charges are present"
+
         #### generate the text
         text = ""
         if append:
@@ -185,6 +191,11 @@ def write_molecule_to_file(cls, filename, molecule, route_card, link0={"mem": "3
         if footer is not None:
             text += f"{footer.strip()}\n\n"
 
+        if point_charges is not None:
+            for point_charge in point_charges:
+                text += f"{point_charge.coordinates[0]:>13.8f} {point_charge.coordinates[1]:>13.8f} {point_charge.coordinates[2]:>13.8f} {point_charge.charge:.5f}\n"
+            text += "\n"
+
         #### write the file
         if append:
             super().append_to_file(filename, text)

cctk/point_charge.py

@@ -0,0 +1,18 @@
+import numpy as np
+
+class PointCharge():
+    """
+    Represents a point charge.
+
+    Attributes:
+        coordinates (np.ndarray): 3-element ndarray
+        charge (float): charge
+    """
+
+    def __init__(self, coordinates, charge):
+        assert isinstance(coordinates, (np.ndarray, list)), "coordinates must be list or ndarray!"
+        assert len(coordinates) == 3, "coordinates must have len 3!"
+        self.coordinates = np.array(coordinates)
+
+        assert isinstance(charge, (float, int)), "charge must be numeric"
+        self.charge = float(charge)

test/test_gaussian.py

@@ -193,5 +193,16 @@ def test_pathological(self):
 #        print(file)
 #        self.assertTrue(isinstance(file[0], cctk.GaussianFile))
 
+    def test_point_charge(self):
+        path = "test/static/Li.out"
+        mol = cctk.GaussianFile.read_file(path).get_molecule()
+
+        point_charge = cctk.PointCharge(coordinates=[0,1,0], charge=-1)
+        self.assertTrue(isinstance(point_charge, cctk.PointCharge))
+
+        new_path = "test/static/Li_pc.gjf"
+        cctk.GaussianFile.write_molecule_to_file(new_path, mol, route_card="#p opt b3lyp/6-31gd charge", point_charges=[point_charge])
+        os.remove(new_path)
+
 if __name__ == '__main__':
     unittest.main()