Added truly custom solvation radii

cyllab · Feb 1, 2022 · 9c0ee3d · 9c0ee3d
1 parent 345d8b0
commit 9c0ee3d
Show file tree

Hide file tree

Showing 7 changed files with 536 additions and 61 deletions.
diff --git a/ccinput/calculation.py b/ccinput/calculation.py
@@ -94,24 +94,25 @@ class Parameters:
         (e.g. regarding the charge and multiplicity) and can be reused.
     """
 
-    def __init__(self, software, solvent="", solvation_model="",
-            solvation_radii="", basis_set="", method="", specifications="",
-            density_fitting="", custom_basis_sets="", d3=False,
-            d3bj=False, **kwargs):
+    def __init__(self, software, solvent="", solvation_model="", solvation_radii="",
+            custom_solvation_radii="", basis_set="", method="", specifications="",
+            density_fitting="", custom_basis_sets="", d3=False, d3bj=False, **kwargs):
 
         self.solvent = get_abs_solvent(solvent)
         if self.solvent != "":
             self.solvation_model = solvation_model.lower()
             self.solvation_radii = solvation_radii.lower()
+            self.custom_solvation_radii = custom_solvation_radii.lower()
 
             if self.solvation_model.strip() == "":
-                raise InvalidParameter("No solvation model specified," +
+                raise InvalidParameter("No solvation model specified, " +
                             "although solvation is requested")
             if self.solvation_radii.strip() == "":
                 warn("No solvation radii specified; using default radii")
         else:
             self.solvation_model = ""
             self.solvation_radii = ""
+            self.custom_solvation_radii = ""
 
         self.software = get_abs_software(software)
 

diff --git a/ccinput/exceptions.py b/ccinput/exceptions.py
@@ -17,3 +17,6 @@ class MissingParameter(CCInputException):
 
 class InternalError(CCInputException):
     pass
+
+class UnimplementedError(CCInputException):
+    pass
diff --git a/ccinput/packages/gaussian.py b/ccinput/packages/gaussian.py
@@ -3,7 +3,7 @@
 
 from ccinput.utilities import get_method, get_solvent, get_basis_set, clean_xyz, \
                               get_distance, get_angle, get_dihedral, get_npxyz
-from ccinput.constants import CalcType, ATOMIC_NUMBER
+from ccinput.constants import CalcType, ATOMIC_NUMBER, LOWERCASE_ATOMIC_SYMBOLS
 from ccinput.exceptions import InvalidParameter
 
 class GaussianCalculation:
@@ -292,30 +292,72 @@ def handle_xyz(self):
         lines = [i + '\n' for i in clean_xyz(self.calc.xyz).split('\n') if i != '' ]
         self.xyz_structure = ''.join(lines)
 
+    def get_custom_solvation_radii_appendix(self):
+        radii_appendix = ""
+        for radius in self.calc.parameters.custom_solvation_radii.split(';'):
+            if radius.strip() == "":
+                continue
+            sradius = radius.split('=')
+            if len(sradius) != 2:
+                raise InvalidParameter(f"Invalid custom solvation radius specification: '{radius}': must follow the pattern '<atom1>=<radius1>;...'")
+
+            element, rad = sradius
+            if element not in LOWERCASE_ATOMIC_SYMBOLS:
+                raise InvalidParameter(f"Invalid element in custom solvation radius specifications: '{element}'")
+
+            _element = LOWERCASE_ATOMIC_SYMBOLS[element] # Add the proper case back
+
+            try:
+                _rad = float(rad)
+            except ValueError:
+                raise InvalidParameter(f"Invalid custom solvation radius for element {element}: '{rad}'")
+            radii_appendix += f"{_element} {_rad:.2f}\n"
+
+        return radii_appendix
+
     def handle_solvation(self):
         if self.calc.parameters.solvent.lower() not in ["", "vacuum"]:
             solvent_keyword = get_solvent(self.calc.parameters.solvent, self.calc.parameters.software, solvation_model=self.calc.parameters.solvation_model)
-            if self.calc.parameters.solvation_model == "smd":
-                if self.calc.parameters.solvation_radii in ["",  "default"]:
-                    self.command_line += f"SCRF(SMD, Solvent={solvent_keyword}) "
-                else:
-                    self.command_line += f"SCRF(SMD, Solvent={solvent_keyword}, Read) "
-                    self.appendix.append(self.SMD18_APPENDIX)
-            elif self.calc.parameters.solvation_model == "pcm":
-                if self.calc.parameters.solvation_radii in ["uff", ""]:
-                    self.command_line += f"SCRF(PCM, Solvent={solvent_keyword}) "
-                else:
-                    self.command_line += f"SCRF(PCM, Solvent={solvent_keyword}, Read) "
-                    self.appendix.append(f"Radii={self.calc.parameters.solvation_radii}\n")
-            elif self.calc.parameters.solvation_model == "cpcm":
-                if self.calc.parameters.solvation_radii in ["uff", ""]:
-                    self.command_line += f"SCRF(CPCM, Solvent={solvent_keyword}) "
+            model = self.calc.parameters.solvation_model
+            radii_set = self.calc.parameters.solvation_radii
+            custom_radii = self.calc.parameters.custom_solvation_radii
+
+            DEFAULT_RADII_SETS = {
+                        'smd': ["", "default"],
+                        'pcm': ["", "uff"],
+                        'cpcm': ["", "uff"],
+            }
+
+            radii_appendix = self.get_custom_solvation_radii_appendix()
+            solv_appendix = ""
+
+            if radii_set not in DEFAULT_RADII_SETS[model] or radii_appendix != "":
+                self.command_line += f"SCRF({model.upper()}, Solvent={solvent_keyword}, Read) "
+            else:
+                self.command_line += f"SCRF({model.upper()}, Solvent={solvent_keyword}) "
+
+            if model == 'smd':
+                if radii_set == "smd18":
+                    solv_appendix += self.SMD18_APPENDIX
+                    solv_appendix += radii_appendix
                 else:
-                    self.command_line += f"SCRF(CPCM, Solvent={solvent_keyword}, Read) "
-                    self.appendix.append(f"Radii={self.calc.parameters.solvation_radii}\n")
+                    if radii_appendix != "":
+                        solv_appendix += "modifysph\n\n"
+                        solv_appendix += radii_appendix
+
+            elif model in ['pcm', 'cpcm']:
+                if radii_set not in DEFAULT_RADII_SETS[model]:
+                    solv_appendix += f"Radii={radii_set}\n"
+
+                if radii_appendix != "":
+                    solv_appendix += "modifysph\n\n"
+                    solv_appendix += radii_appendix
             else:
                 raise InvalidParameter(f"Invalid solvation method for Gaussian: '{self.calc.parameters.solvation_model}'")
 
+            self.appendix.append(solv_appendix)
+
+
     def create_input_file(self):
         additional_commands = " ".join([i.strip() for i in self.additional_commands]).strip()
         self.confirmed_specifications += additional_commands

diff --git a/ccinput/packages/orca.py b/ccinput/packages/orca.py
@@ -1,6 +1,6 @@
 import basis_set_exchange
 
-from ccinput.constants import CalcType, ATOMIC_NUMBER
+from ccinput.constants import CalcType, ATOMIC_NUMBER, LOWERCASE_ATOMIC_SYMBOLS
 from ccinput.utilities import get_method, get_basis_set, get_solvent, clean_xyz
 from ccinput.exceptions import InvalidParameter
 
@@ -265,38 +265,72 @@ def handle_pal(self):
 
         self.blocks.append(pal_block)
 
+    def get_custom_solvation_radii_specs(self):
+        radii_specs = ""
+        for radius in self.calc.parameters.custom_solvation_radii.split(';'):
+            if radius.strip() == "":
+                continue
+            sradius = radius.split('=')
+            if len(sradius) != 2:
+                raise InvalidParameter(f"Invalid custom solvation radius specification: '{radius}': must follow the pattern '<atom1>=<radius1>;...'")
+
+            element, rad = sradius
+            if element not in LOWERCASE_ATOMIC_SYMBOLS:
+                raise InvalidParameter(f"Invalid element in custom solvation radius specifications: '{element}'")
+
+            _element = ATOMIC_NUMBER[LOWERCASE_ATOMIC_SYMBOLS[element]]
+
+            try:
+                _rad = float(rad)
+            except ValueError:
+                raise InvalidParameter(f"Invalid custom solvation radius for element {element}: '{rad}'")
+            radii_specs += f"radius[{_element}] {_rad:.2f}\n"
+
+        return radii_specs
+
     def handle_solvation(self):
         if self.calc.parameters.solvent.lower() not in ["vacuum", ""]:
             solvent_keyword = get_solvent(self.calc.parameters.solvent,
                                           self.calc.parameters.software,
                                           solvation_model=self.calc.parameters.solvation_model)
+            model = self.calc.parameters.solvation_model
+            radii_set = self.calc.parameters.solvation_radii
+            custom_radii = self.calc.parameters.custom_solvation_radii
+
+            solv_block = ""
 
             if self.calc.parameters.method[:3] == 'gfn':
                 self.command_line += f" ALPB({solvent_keyword})"
-            elif self.calc.parameters.solvation_model == "smd":
-                if self.calc.parameters.solvation_radii in ["default", ""]:
-                    smd_block = f'''%cpcm
-                    smd true
-                    SMDsolvent "{solvent_keyword}"
-                    end'''
-                    self.blocks.append(smd_block)
+            elif model == "smd":
+                radii_specs = ""
 
                 # Refined solvation radii
                 # E. Engelage, N. Schulz, F. Heinen, S. M. Huber, D. G. Truhlar,
                 # C. J. Cramer, Chem. Eur. J. 2018, 24, 15983–15987.
-                elif self.calc.parameters.solvation_radii == "smd18":
-                    smd_block = f'''%cpcm
-                    smd true
-                    SMDsolvent "{solvent_keyword}"
-                    radius[53] 2.74
-                    radius[35] 2.60
-                    end'''
-                    self.blocks.append(smd_block)
-            elif self.calc.parameters.solvation_model == "cpcm":
+                if model == "smd" and radii_set == "smd18":
+                    radii_specs += "radius[53] 2.74\nradius[35] 2.60\n"
+
+                if custom_radii != "":
+                    radii_specs += self.get_custom_solvation_radii_specs()
+
+                solv_block = f'''%cpcm
+                smd true
+                SMDsolvent "{solvent_keyword}"
+                {radii_specs}end'''
+                self.blocks.append(solv_block)
+            elif model == "cpcm":
                 self.command_line += f"CPCM({solvent_keyword}) "
-                ###CPCM radii
+                radii_specs = self.get_custom_solvation_radii_specs()
+                if radii_specs != "":
+                    solv_block = f'''%cpcm
+                    {radii_specs}end'''
+                    self.blocks.append(solv_block)
+                if radii_set not in ["", "default", "bondi"]:
+                    raise UnimplementedError("As of now, ccinput does not know how to request "
+                            "other solvation radii than the default ones in ORCA with CPCM")
             else:
-                raise InvalidParameter(f"Invalid solvation model for ORCA: '{self.calc.parameters.solvation_model}'")
+                raise InvalidParameter(f"Invalid solvation model for ORCA: "
+                        "'{self.calc.parameters.solvation_model}'")
 
     def create_input_file(self):
         self.block_lines = '\n'.join(self.blocks)