Merge pull request #220 from ahurta92/adrian-Dipole

NWChem: Run dipole moment calculation

.gitignore

@@ -74,3 +74,6 @@ timer.dat
 *.out
 iter.dat
 psi*clean
+
+# vscode folders
+.vscode

qcengine/programs/nwchem/germinate.py

@@ -105,12 +105,7 @@ def muster_modelchem(method: str, derint: int, use_tce: bool) -> Tuple[str, Dict
  opts = {}
 
  # Map the run type to
- runtyp = {
- 0: "energy",
- 1: "gradient",
- 2: "hessian",
- # 'properties': 'prop',
- }[derint]
+ runtyp = {"energy": "energy", "gradient": "gradient", "hessian": "hessian", "properties": "property"}[derint]
 
  # Write out the theory directive
  if method == "nwchem":

qcengine/programs/nwchem/harvester.py

@@ -614,11 +614,94 @@ def harvest_outfile_pass(outtext):
  if mobj:
  psivar["N ALPHA ELECTRONS"] = mobj.group(2)
  psivar["N BETA ELECTRONS"] = mobj.group(3)
+
+ if psivar["N ALPHA ELECTRONS"] == psivar["N BETA ELECTRONS"]:
+
+ # get HOMO and LUMO energy
+ mobj = re.search(
+ r"Vector"
+ + r"\s+"
+ + r"%d" % (psivar["N ALPHA ELECTRONS"])
+ + r"\s+"
+ + r"Occ="
+ + r".*"
+ + r"\s+"
+ + r"E="
+ + r"([+-]?\s?\d+[.]\d+)"
+ + r"[D]"
+ + r"([+-])"
+ + r"[0]"
+ + r"(\d+)",
+ outtext,
+ re.MULTILINE,
+ )
+ if mobj:
+ homo = float(mobj.group(1)) * (10 ** (-1 * float(mobj.group(3))))
+ psivar["HOMO"] = np.array([round(homo, 10)])
+ mobj = re.search(
+ r"Vector"
+ + r"\s+"
+ + r"%d" % (psivar["N ALPHA ELECTRONS"] + 1)
+ + r"\s+"
+ + r"Occ="
+ + r".*"
+ + r"\s+"
+ + r"E="
+ + r"([+-]?\s?\d+[.]\d+)"
+ + r"[D]"
+ + r"([+-])"
+ + r"[0]"
+ + r"(\d+)",
+ outtext,
+ re.MULTILINE,
+ )
+ if mobj:
+ lumo = float(mobj.group(1)) * (10 ** (-1 * float(mobj.group(3))))
+ psivar["LUMO"] = np.array([round(lumo, 10)])
+
  mobj = re.search(r"AO basis - number of functions:\s+(\d+)\s+number of shells:\s+(\d+)", outtext, re.MULTILINE)
  if mobj:
  psivar["N MO"] = mobj.group(2)
  psivar["N BASIS"] = mobj.group(1)
 
+ # Search for Center of charge
+ mobj = re.search(
+ r"Center of charge \(in au\) is the expansion point"
+ + r"\n"
+ + r"\s+"
+ + r"X\s+=\s+([+-]?\d+[.]\d+)"
+ + r"\s+"
+ + r"Y\s+=\s+([+-]?\d+[.]\d+)"
+ + r"\s+"
+ + r"Z\s+=\s+([+-]?\d+[.]\d+)",
+ outtext,
+ re.MULTILINE,
+ )
+ if mobj:
+ psivar["CENTER OF CHARGE"] = np.array([mobj.group(1), mobj.group(2), mobj.group(3)])
+
+ mobj = re.search(
+ r"Dipole moment"
+ + r".*?"
+ + r"A\.U\."
+ + r"\s+"
+ + r"DMX\s+([+-]?\d+[.]\d+)\s+"
+ + r"DMXEFC\s+[+-]?\d+[.]\d+\s+"
+ + r"DMY\s+([+-]?\d+[.]\d+)\s+"
+ + r"DMYEFC\s+[+-]?\d+[.]\d+\s+"
+ + r"DMZ\s+([+-]?\d+[.]\d+)\s+"
+ + r"DMZEFC\s+[+-]?\d+[.]\d+\s+"
+ + r"\-EFC\-"
+ + r".*?"
+ + r"A\.U\.\s+"
+ + r"Total dipole\s+([+-]?\d+[.]\d+\s+)",
+ outtext,
+ re.MULTILINE,
+ )
+ # + r"DMY\s+" + r"([+-]?\d+[.]\d+)", outtext, re.MULTILINE)
+ if mobj:
+ psivar["DIPOLE MOMENT"] = np.array([mobj.group(1), mobj.group(2), mobj.group(3)])
+ psivar["TOTAL DIPOLE MOMENT"] = mobj.group(4)
  # Process CURRENT energies (TODO: needs better way)
  if "HF TOTAL ENERGY" in psivar:
  psivar["SCF TOTAL ENERGY"] = psivar["HF TOTAL ENERGY"]

qcengine/programs/nwchem/runner.py

@@ -155,9 +155,7 @@ def build_input(
  opts.update(moldata["keywords"])
 
  # Handle calc type and quantum chemical method
- mdccmd, mdcopts = muster_modelchem(
- input_model.model.method, input_model.driver.derivative_int(), opts.pop("qc_module", False)
- )
+ mdccmd, mdcopts = muster_modelchem(input_model.model.method, input_model.driver, opts.pop("qc_module", False))
  opts.update(mdcopts)
 
  # Handle basis set
@@ -238,7 +236,11 @@ def parse_output(
  qcvars["CURRENT HESSIAN"] = nwhess
 
  # Normalize the output as a float or list of floats
- retres = qcvars[f"CURRENT {input_model.driver.upper()}"]
+ if input_model.driver.upper() == "PROPERTIES":
+ retres = qcvars[f"CURRENT ENERGY"]
+ else:
+ retres = qcvars[f"CURRENT {input_model.driver.upper()}"]
+
  if isinstance(retres, Decimal):
  retres = float(retres)
  elif isinstance(retres, np.ndarray):

qcengine/programs/tests/test_nwchem.py

@@ -97,3 +97,48 @@ def test_gradient(nh2):
  shif_det = np.linalg.det(shif_grads)
 
  assert np.allclose(orig_det, shif_det)
+
+
+@pytest.fixture
+def h20():
+ water = """
+-1 2
+O 0 0 0
+H 0 0 1
+H 0 1 0
+ """
+ return qcel.models.Molecule.from_data(water)
+
+
+@using("nwchem")
+def test_dipole(h20):
+ # Run NH2
+ resi = {
+ "molecule": h20,
+ "driver": "properties",
+ "model": {"method": "dft", "basis": "3-21g"},
+ "keywords": {"dft__xc": "b3lyp", "property__dipole": True},
+ }
+ res = qcng.compute(resi, "nwchem", raise_error=True, return_dict=True)
+
+ # Make sure the calculation completed successfully
+ assert compare_values(-75.764944, res["return_result"], atol=1e-3)
+ assert res["driver"] == "properties"
+ assert "provenance" in res
+ assert res["success"] is True
+
+ # Check the other status information
+ assert res["extras"]["qcvars"]["N ALPHA ELECTRONS"] == "6"
+ assert res["extras"]["qcvars"]["N ATOMS"] == "3"
+ assert res["extras"]["qcvars"]["N BASIS"] == "13"
+
+ # Make sure the properties parsed correctly
+ assert compare_values(-75.764944, res["properties"]["return_energy"], atol=1e-3)
+ assert res["properties"]["calcinfo_natom"] == 3
+ assert res["properties"]["calcinfo_nalpha"] == 6
+ assert res["properties"]["calcinfo_nbasis"] == 13
+ # Make sure Dipole Moment and center of charge parsed correctly
+ assert compare_values(0.272949872, float(res["extras"]["qcvars"]["TOTAL DIPOLE MOMENT"]), atol=1e-5)
+ assert compare_values(-0.00, float(res["extras"]["qcvars"]["DIPOLE MOMENT"][0]), atol=1e-3)
+ assert compare_values(-0.00, float(res["extras"]["qcvars"]["DIPOLE MOMENT"][1]), atol=1e-3)
+ assert compare_values(-0.272949872, float(res["extras"]["qcvars"]["DIPOLE MOMENT"][2]), atol=1e-5)