Added dispersion correction

cyllab · Jan 27, 2022 · d0fdd99 · d0fdd99
1 parent c9b33fd
commit d0fdd99
Show file tree

Hide file tree

Showing 10 changed files with 378 additions and 28 deletions.
diff --git a/ccinput/calculation.py b/ccinput/calculation.py
@@ -5,7 +5,7 @@
                                MissingParameter
 from ccinput.utilities import get_abs_software, get_method, get_abs_basis_set, \
                               get_abs_solvent, get_theory_level, standardize_memory, \
-                              get_npxyz, get_coord
+                              get_npxyz, get_coord, has_dispersion_parameters
 from ccinput.constants import ATOMIC_NUMBER, SYN_SOFTWARE
 from ccinput.logging import warn
 
@@ -96,7 +96,8 @@ class Parameters:
 
     def __init__(self, software, solvent="", solvation_model="",
             solvation_radii="", basis_set="", method="", specifications="",
-            density_fitting="", custom_basis_sets="", **kwargs):
+            density_fitting="", custom_basis_sets="", d3=False,
+            d3bj=False, **kwargs):
 
         self.solvent = get_abs_solvent(solvent)
         if self.solvent != "":
@@ -129,6 +130,20 @@ def __init__(self, software, solvent="", solvation_model="",
         else:
             self.basis_set = ""
 
+        if d3 and d3bj:
+            raise InvalidParameter("Cannot use both D3(0) and D3BJ dispersion corrections")
+
+        self.d3 = d3
+        self.d3bj = d3bj
+
+        if d3 and not has_dispersion_parameters(self.method, version='d3'):
+            warn(f"Your calculation requests a method ({self.method}) that may not have D3 "
+                 "parameters. Be aware that the calculation might result in an error")
+
+        if d3bj and not has_dispersion_parameters(self.method, version='d3bj'):
+            warn(f"Your calculation requests a method ({self.method}) that may not have D3BJ "
+                 "parameters. Be aware that the calculation might result in an error")
+
         self.specifications = specifications
         self.density_fitting = density_fitting
         self.custom_basis_sets = custom_basis_sets

diff --git a/ccinput/constants.py b/ccinput/constants.py
@@ -1955,4 +1955,103 @@ class CalcType(Enum):
                 },
 }
 
+# The two lists below use the "absolute" method identifications.
+# See SYN_METHODS for the full list of names and their synonyms.
+
+# These lists contain all the functionals with D3 parameters developed
+# by the Grimme group, which are not software specific. As such, the
+# parameters might not be available in every QM package (although it could).
+
+# Some functionals were not found in the list of known functionals or synonyms.
+# These are indicated with the question mark.
+FUNCTIONALS_WITH_DISPERSION_PARAMETERS = {
+        'd3': [
+                'b1b95',
+                'b2gp-plyp',
+                'b3lyp',
+                'b97-d',
+                'bhlyp', # BHandHLYP?
+                'blyp',
+                'bp86', # ?
+                'bpbe', # HISSbPBE?
+                'mpwlyp',
+                'pbe',
+                'pbe0',
+                'pw6b95',
+                'pwb6k', # ?
+                'revpbe',
+                'tpss',
+                'tpss0',
+                'tpssh',
+                'bop', # ?
+                'mpw1b95', # ?
+                'mpwb1k', # ?
+                'olyp',
+                'opbe',
+                'otpss', # ?
+                'pbe38', # ?
+                'pbesol', # ?
+                'revssb', # ?
+                'ssb', # ?
+                'b3pw91', # ?
+                'bmk',
+                'cam-b3lyp',
+                'lcwpbe', # ?
+                'm052x',
+                'm05',
+                'm062x',
+                'm06hf',
+                'm06l',
+                'm06',
+                'hcth120', # ?
+                'b2plyp',
+                'dsdblyp', # ?
+                'ptpss', # ?
+                'pwpb95',
+                'revpbe0', # ?
+                'revpbe38', # ?
+                'rpw86pbe',
+            ],
+        'd3bj': [
+                '1b95',
+                'b2gp-plyp',
+                'b3pw91', # ?
+                'bhlyp', # ?
+                'bmk',
+                'bop', # ?
+                'bpbe', # ?
+                'cam-b3lyp',
+                'lcwpbe', # ?
+                'mpw1b95', # ?
+                'mpwb1k', # ?
+                'mpwlyp',
+                'olyp',
+                'opbe',
+                'otpss', # ?
+                'pbe38', # ?
+                'pbesol', # ?
+                'ptpss', # ?
+                'pwb6k', # ?
+                'revssb', # ?
+                'ssb', # ?
+                'tpssh',
+                'hcth120', # ?
+                'b2plyp',
+                'b3lyp',
+                'b97-d',
+                'blyp',
+                'bp86', # ?
+                'dsdblyp', # ?
+                'pbe0',
+                'pbe',
+                'pw6b95',
+                'pwpb95',
+                'revpbe0', # ?
+                'revpbe38', # ?
+                'revpbe',
+                'rpw86pbe',
+                'tpss0',
+                'tpss',
+            ],
+}
 
diff --git a/ccinput/packages/gaussian.py b/ccinput/packages/gaussian.py
@@ -124,6 +124,11 @@ def add_spec(key, option):
             spec_formatted = f'{spec}({specs_str}) '
             self.additional_commands.append(spec_formatted)
 
+        if self.calc.parameters.d3:
+            self.additional_commands.append("EmpiricalDispersion=GD3 ")
+        elif self.calc.parameters.d3bj:
+            self.additional_commands.append("EmpiricalDispersion=GD3BJ ")
+
     def handle_command(self):
         cmd = ""
         base_specs = []

diff --git a/ccinput/packages/orca.py b/ccinput/packages/orca.py
@@ -62,31 +62,36 @@ def clean(self, s):
 
     def handle_specifications(self):
         _specifications = self.clean(self.calc.parameters.specifications).lower().strip()
-        if _specifications == '':
-            return
 
         specifications_list = []
-        sspecs = _specifications.split()
-        ind = 0
-        while ind < len(sspecs):
-            spec = sspecs[ind]
-            if spec == "--phirshfeld":
-                HIRSHFELD_BLOCK = """%output
-                Print[ P_Hirshfeld] 1
-                end"""
-                self.blocks.append(HIRSHFELD_BLOCK)
-            elif spec == "--nimages":
-                nimages = sspecs[ind+1]
-                try:
-                    nimages = int(nimages)
-                except ValueError:
-                    raise InvalidParameter("Invalid specifications")
-                self.specifications['nimages'] = nimages
+
+        if _specifications != '':
+            sspecs = _specifications.split()
+            ind = 0
+            while ind < len(sspecs):
+                spec = sspecs[ind]
+                if spec == "--phirshfeld":
+                    HIRSHFELD_BLOCK = """%output
+                    Print[ P_Hirshfeld] 1
+                    end"""
+                    self.blocks.append(HIRSHFELD_BLOCK)
+                elif spec == "--nimages":
+                    nimages = sspecs[ind+1]
+                    try:
+                        nimages = int(nimages)
+                    except ValueError:
+                        raise InvalidParameter("Invalid specifications")
+                    self.specifications['nimages'] = nimages
+                    ind += 1
+                elif spec not in specifications_list:
+                    specifications_list.append(spec)
+
                 ind += 1
-            elif spec not in specifications_list:
-                specifications_list.append(spec)
 
-            ind += 1
+        if self.calc.parameters.d3:
+            specifications_list.append("d3zero")
+        elif self.calc.parameters.d3bj:
+            specifications_list.append("d3bj")
 
         if len(specifications_list) > 0:
             self.additional_commands = " ".join(specifications_list)

diff --git a/ccinput/tests/test_cli.py b/ccinput/tests/test_cli.py
@@ -289,3 +289,63 @@ def test_opt_freq2(self):
         line = "gaussian opt-freq HF -bs Def2SVP --xyz 'Cl 0 0 0' -n 1 --mem 1G -c -1"
         self.assertTrue(self.are_equivalent(args, line))
 
+    def test_d3_gaussian(self):
+        args = {
+            'software': "gaussian",
+            'type': "opt",
+            'method': "M06",
+            'basis_set': "Def2SVP",
+            'xyz': "Cl 0 0 0\n",
+            'nproc': 1,
+            'mem': "1G",
+            'charge': -1,
+            'd3': True,
+        }
+        line = "gaussian opt M06 -bs Def2SVP --xyz 'Cl 0 0 0' -n 1 --mem 1G -c -1 --d3"
+        self.assertTrue(self.are_equivalent(args, line))
+
+    def test_d3bj_gaussian(self):
+        args = {
+            'software': "gaussian",
+            'type': "opt",
+            'method': "PBE0",
+            'basis_set': "Def2SVP",
+            'xyz': "Cl 0 0 0\n",
+            'nproc': 1,
+            'mem': "1G",
+            'charge': -1,
+            'd3bj': True,
+        }
+        line = "gaussian opt PBE0 -bs Def2SVP --xyz 'Cl 0 0 0' -n 1 --mem 1G -c -1 --d3bj"
+        self.assertTrue(self.are_equivalent(args, line))
+
+    def test_d3_orca(self):
+        args = {
+            'software': "ORCA",
+            'type': "opt",
+            'method': "M06",
+            'basis_set': "Def2SVP",
+            'xyz': "Cl 0 0 0\n",
+            'nproc': 1,
+            'mem': "1G",
+            'charge': -1,
+            'd3': True,
+        }
+        line = "orca opt M06 -bs Def2SVP --xyz 'Cl 0 0 0' -n 1 --mem 1G -c -1 --d3"
+        self.assertTrue(self.are_equivalent(args, line))
+
+    def test_d3bj_orca(self):
+        args = {
+            'software': "ORCA",
+            'type': "opt",
+            'method': "PBE0",
+            'basis_set': "Def2SVP",
+            'xyz': "Cl 0 0 0\n",
+            'nproc': 1,
+            'mem': "1G",
+            'charge': -1,
+            'd3bj': True,
+        }
+        line = "orca opt PBE0 -bs Def2SVP --xyz 'Cl 0 0 0' -n 1 --mem 1G -c -1 --d3bj"
+        self.assertTrue(self.are_equivalent(args, line))
+
diff --git a/ccinput/tests/test_gaussian.py b/ccinput/tests/test_gaussian.py
@@ -1,4 +1,5 @@
 from ccinput.tests.testing_utilities import InputTests
+from ccinput.exceptions import InvalidParameter
 
 class GaussianTests(InputTests):
 
@@ -2587,3 +2588,80 @@ def test_opt_freq2(self):
 
         self.assertTrue(self.is_equivalent(REF, inp.input_file))
 
+    def test_d3(self):
+        params = {
+                'nproc': 8,
+                'mem': '10GB',
+                'type': 'Single-Point Energy',
+                'in_file': 'Cl.xyz',
+                'software': 'Gaussian',
+                'method': 'M062X',
+                'basis_set': '3-21G',
+                'charge': '-1',
+                'd3': True,
+                }
+
+        inp = self.generate_calculation(**params)
+
+        REF = """
+        %chk=calc.chk
+        %nproc=8
+        %mem=10000MB
+        #p sp M062X/3-21G EmpiricalDispersion=GD3
+
+        File created by ccinput
+
+        -1 1
+        Cl 0.0 0.0 0.0
+
+        """
+
+        self.assertTrue(self.is_equivalent(REF, inp.input_file))
+
+    def test_d3bj(self):
+        params = {
+                'nproc': 8,
+                'mem': '10GB',
+                'type': 'Single-Point Energy',
+                'in_file': 'Cl.xyz',
+                'software': 'Gaussian',
+                'method': 'PBE0',
+                'basis_set': '3-21G',
+                'charge': '-1',
+                'd3bj': True,
+                }
+
+        inp = self.generate_calculation(**params)
+
+        REF = """
+        %chk=calc.chk
+        %nproc=8
+        %mem=10000MB
+        #p sp PBE1PBE/3-21G EmpiricalDispersion=GD3BJ
+
+        File created by ccinput
+
+        -1 1
+        Cl 0.0 0.0 0.0
+
+        """
+
+        self.assertTrue(self.is_equivalent(REF, inp.input_file))
+
+    def test_d3_d3bj_crash(self):
+        params = {
+                'nproc': 8,
+                'mem': '10GB',
+                'type': 'Single-Point Energy',
+                'in_file': 'Cl.xyz',
+                'software': 'Gaussian',
+                'method': 'PBE0',
+                'basis_set': '3-21G',
+                'charge': '-1',
+                'd3': True,
+                'd3bj': True,
+                }
+
+        with self.assertRaises(InvalidParameter):
+            self.generate_calculation(**params)
+