diff --git a/adcc/OperatorIntegrals.py b/adcc/OperatorIntegrals.py index 270529eb..0bfb0b9a 100644 --- a/adcc/OperatorIntegrals.py +++ b/adcc/OperatorIntegrals.py @@ -102,7 +102,8 @@ def available(self): "electric_dipole", "magnetic_dipole", "nabla", - "pe_induction_elec" + "pe_induction_elec", + "pcm_potential_elec" ) return [integral for integral in integrals if hasattr(self.provider_ao, integral)] @@ -187,6 +188,19 @@ def pe_induction_elec(self): callback = self.provider_ao.pe_induction_elec return self.__import_density_dependent_operator(callback) + @property + def pcm_potential_elec(self): + """ + Returns a function to obtain the (density-dependent) + electronic PCM potential operator in the molecular orbital basis + """ + if "pcm_potential_elec" not in self.available: + raise NotImplementedError("electronic PCM potential operator " + "not implemented " + f"in {self.provider_ao.backend} backend.") + callback = self.provider_ao.pcm_potential_elec + return self.__import_density_dependent_operator(callback) + @property def timer(self): ret = Timer() diff --git a/adcc/adc_pp/environment.py b/adcc/adc_pp/environment.py index 1ee02d42..aad254e0 100644 --- a/adcc/adc_pp/environment.py +++ b/adcc/adc_pp/environment.py @@ -41,3 +41,22 @@ def apply(ampl): vpe = op.pe_induction_elec(tdm) return AmplitudeVector(ph=vpe.ov) return AdcBlock(apply, 0) + + +def block_ph_ph_0_pcm(hf, mp, intermediates): + """ + Constructs an :py:class:`AdcBlock` that describes the + linear response coupling to the polarizable environment + from PCM via a CIS-like transition density as described + in 10.1021/ct300763v, eq 63. Since the contribution + depends on the input amplitude itself, + a diagonal term cannot be formulated. + """ + op = hf.operators + + def apply(ampl): + tdm = OneParticleOperator(mp, is_symmetric=False) + tdm.vo = ampl.ph.transpose() + vpcm = op.pcm_potential_elec(tdm) + return AmplitudeVector(ph=vpcm.ov) + return AdcBlock(apply, 0) diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py index 94b3aa24..4d72844f 100644 --- a/adcc/backends/psi4.py +++ b/adcc/backends/psi4.py @@ -64,6 +64,16 @@ def pe_induction_elec_ao(dm): )[1] return pe_induction_elec_ao + @property + def pcm_potential_elec(self): + if self.wfn.PCM_enabled(): + def pcm_potential_elec_ao(dm): + return psi4.core.PCM.compute_V( + self.wfn.get_PCM(), + psi4.core.Matrix.from_array(dm.to_ndarray()) + ) + return pcm_potential_elec_ao + class Psi4EriBuilder(EriBuilder): def __init__(self, wfn, n_orbs, n_orbs_alpha, n_alpha, n_beta, restricted): @@ -106,6 +116,13 @@ def pe_energy(self, dm, elec_only=True): elec_only=elec_only) return e_pe + def pcm_energy(self, dm): + psi_dm = psi4.core.Matrix.from_array(dm.to_ndarray()) + # computes the Fock matrix contribution + # By contraction with the tdm, the electronic contribution is obtained + V_pcm = psi4.core.PCM.compute_V(self.wfn.get_PCM(), psi_dm).to_array() + return np.einsum("uv,uv->", dm.to_ndarray(), V_pcm) + @property def excitation_energy_corrections(self): ret = [] @@ -121,6 +138,12 @@ def excitation_energy_corrections(self): elec_only=True) ) ret.extend([ptlr, ptss]) + if self.environment == "pcm": + ptlr = EnergyCorrection( + "pcm_ptlr_correction", + lambda view: self.pcm_energy(view.transition_dm_ao) + ) + ret.extend([ptlr]) return {ec.name: ec for ec in ret} @property @@ -128,6 +151,8 @@ def environment(self): ret = None if hasattr(self.wfn, "pe_state"): ret = "pe" + elif self.wfn.PCM_enabled(): + ret = "pcm" return ret def get_backend(self): @@ -245,6 +270,8 @@ def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11, "ccpvdz": "cc-pvdz", } + # needed for PE and PCM tests + psi4.core.clean_options() mol = psi4.geometry(f""" {charge} {multiplicity} {xyz} diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py index 90128750..f3d947e4 100644 --- a/adcc/backends/pyscf.py +++ b/adcc/backends/pyscf.py @@ -29,7 +29,8 @@ from ..exceptions import InvalidReference from ..ExcitedStates import EnergyCorrection -from pyscf import ao2mo, gto, scf, solvent +from pyscf import ao2mo, gto, scf +from pyscf.solvent import pol_embed, ddcosmo class PyScfOperatorIntegralProvider: @@ -59,13 +60,28 @@ def nabla(self): @property def pe_induction_elec(self): if hasattr(self.scfres, "with_solvent"): - if isinstance(self.scfres.with_solvent, solvent.pol_embed.PolEmbed): + if isinstance(self.scfres.with_solvent, pol_embed.PolEmbed): def pe_induction_elec_ao(dm): return self.scfres.with_solvent._exec_cppe( dm.to_ndarray(), elec_only=True )[1] return pe_induction_elec_ao + @property + def pcm_potential_elec(self): + if hasattr(self.scfres, "with_solvent"): + if isinstance(self.scfres.with_solvent, ddcosmo.DDCOSMO): + def pcm_potential_elec_ao(dm): + # Since eps (dielectric constant) is the only solvent parameter + # in pyscf and there is no solvent data available in the + # program, the user needs to adjust scfres.with_solvent.eps + # manually to the optical dielectric constant (if non + # equilibrium solvation is desired). + return self.scfres.with_solvent._B_dot_x( + dm.to_ndarray() + ) + return pcm_potential_elec_ao + # TODO: refactor ERI builder to be more general # IntegralBuilder would be good @@ -125,6 +141,15 @@ def pe_energy(self, dm, elec_only=True): e_pe, _ = pe_state.kernel(dm.to_ndarray(), elec_only=elec_only) return e_pe + def pcm_energy(self, dm): + # Since eps (dielectric constant) is the only solvent parameter + # in pyscf and there is no solvent data available in the + # program, the user needs to adjust scfres.with_solvent.eps + # manually to the optical dielectric constant (if non + # equilibrium solvation is desired). + V_pcm = self.scfres.with_solvent._B_dot_x(dm.to_ndarray()) + return np.einsum("uv,uv->", dm.to_ndarray(), V_pcm) + @property def excitation_energy_corrections(self): ret = [] @@ -140,14 +165,22 @@ def excitation_energy_corrections(self): elec_only=True) ) ret.extend([ptlr, ptss]) + elif self.environment == "pcm": + ptlr = EnergyCorrection( + "pcm_ptlr_correction", + lambda view: self.pcm_energy(view.transition_dm_ao) + ) + ret.extend([ptlr]) return {ec.name: ec for ec in ret} @property def environment(self): ret = None if hasattr(self.scfres, "with_solvent"): - if isinstance(self.scfres.with_solvent, solvent.pol_embed.PolEmbed): + if isinstance(self.scfres.with_solvent, pol_embed.PolEmbed): ret = "pe" + elif isinstance(self.scfres.with_solvent, ddcosmo.DDCOSMO): + ret = "pcm" return ret def get_backend(self): diff --git a/adcc/backends/test_backends_pcm.py b/adcc/backends/test_backends_pcm.py new file mode 100644 index 00000000..f98090b2 --- /dev/null +++ b/adcc/backends/test_backends_pcm.py @@ -0,0 +1,206 @@ +import pytest +import unittest +import itertools +import adcc +import adcc.backends +import os + +from adcc.misc import expand_test_templates +from adcc.testdata import static_data +from adcc.testdata.cache import psi4_data, pyscf_data +from numpy.testing import assert_allclose +from adcc.exceptions import InputError + +from adcc.adc_pp.environment import block_ph_ph_0_pcm +from adcc.AdcMatrix import AdcExtraTerm + +backends = [b for b in adcc.backends.available() + if b in ["psi4", "pyscf"]] +basissets = ["sto-3g", "cc-pvdz"] +methods = ["adc1", "adc2", "adc3"] + + +@pytest.mark.skipif(len(backends) == 0, reason="Psi4 and Pyscf backends not found.") +@expand_test_templates(list(itertools.product(basissets, methods, backends))) +class TestPCM(unittest.TestCase): + def template_pcm_ptlr_formaldehyde(self, basis, method, backend): + if method != "adc1": + pytest.skip("Data only available for adc1.") + + if backend == "psi4": + data = psi4_data + run_hf = psi4_run_pcm_hf + pcm_options = {"weight": 0.3, "pcm_method": "IEFPCM", "neq": True, + "solvent": "Water"} + elif backend == "pyscf": + data = pyscf_data + run_hf = pyscf_run_pcm_hf + pcm_options = {"eps": 78.3553, "eps_opt": 1.78} + basename = f"formaldehyde_{basis}_pcm_{method}" + result = data[basename] + + scfres = run_hf(static_data.xyz["formaldehyde"], basis, charge=0, + multiplicity=1, conv_tol=1e-12, conv_tol_grad=1e-11, + max_iter=150, pcm_options=pcm_options) + + assert_allclose(scfres.energy_scf, result["energy_scf"], atol=1e-8) + + state = adcc.run_adc(scfres, method=method, n_singlets=5, + conv_tol=1e-7, environment="ptlr") + + # compare ptLR result to LR data + assert_allclose(state.excitation_energy, + result["lr_excitation_energy"], atol=5 * 1e-3) + + # Consistency check with values obtained with ADCc + assert_allclose(state.excitation_energy, + result["ptlr_adcc_excitation_energy"], atol=1e-6) + + # remove cavity files from PSI4 PCM calculations + if backend == "psi4": + for cavityfile in os.listdir(os.getcwd()): + if cavityfile.startswith(("cavity.off_", "PEDRA.OUT_")): + os.remove(cavityfile) + + def template_pcm_linear_response_formaldehyde(self, basis, method, backend): + if method != "adc1": + pytest.skip("Reference only exists for adc1.") + + if backend == "psi4": + data = psi4_data + run_hf = psi4_run_pcm_hf + pcm_options = {"weight": 0.3, "pcm_method": "IEFPCM", "neq": True, + "solvent": "Water"} + elif backend == "pyscf": + data = pyscf_data + run_hf = pyscf_run_pcm_hf + pcm_options = {"eps": 78.3553, "eps_opt": 1.78} + basename = f"formaldehyde_{basis}_pcm_{method}" + result = data[basename] + + scfres = run_hf(static_data.xyz["formaldehyde"], basis, charge=0, + multiplicity=1, conv_tol=1e-12, conv_tol_grad=1e-11, + max_iter=150, pcm_options=pcm_options) + + assert_allclose(scfres.energy_scf, result["energy_scf"], atol=1e-8) + + matrix = adcc.AdcMatrix(method, scfres) + solvent = AdcExtraTerm(matrix, {'ph_ph': block_ph_ph_0_pcm}) + + matrix += solvent + assert len(matrix.extra_terms) + + state = adcc.run_adc(matrix, n_singlets=5, conv_tol=1e-7, + environment=False) + assert_allclose( + state.excitation_energy_uncorrected, + result["lr_excitation_energy"], + atol=1e-5 + ) + + state_cis = adcc.ExcitedStates(state, property_method="adc0") + assert_allclose( + state_cis.oscillator_strength, + result["lr_osc_strength"], atol=1e-3 + ) + + # invalid combination + with pytest.raises(InputError): + adcc.run_adc(scfres, method=method, n_singlets=5, + environment={"linear_response": True, "ptlr": True}) + + # no environment specified + with pytest.raises(InputError): + adcc.run_adc(scfres, method=method, n_singlets=5) + + # automatically add coupling term + state = adcc.run_adc(scfres, method=method, n_singlets=5, + conv_tol=1e-7, environment="linear_response") + assert_allclose( + state.excitation_energy_uncorrected, + result["lr_excitation_energy"], + atol=1e-5 + ) + + if backend == "psi4": + # remove cavity files from PSI4 PCM calculations + for cavityfile in os.listdir(os.getcwd()): + if cavityfile.startswith(("cavity.off_", "PEDRA.OUT_")): + os.remove(cavityfile) + + +def psi4_run_pcm_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-12, + conv_tol_grad=1e-11, max_iter=150, pcm_options=None): + import psi4 + + # needed for PE and PCM tests + psi4.core.clean_options() + mol = psi4.geometry(f""" + {charge} {multiplicity} + {xyz} + symmetry c1 + units au + no_reorient + no_com + """) + + psi4.core.be_quiet() + psi4.set_options({ + 'basis': basis, + 'scf_type': 'pk', + 'e_convergence': conv_tol, + 'd_convergence': conv_tol_grad, + 'maxiter': max_iter, + 'reference': "RHF", + "pcm": True, + "pcm_scf_type": "total", + }) + psi4.pcm_helper(f""" + Units = AU + Cavity {{Type = Gepol + Area = {pcm_options.get("weight", 0.3)}}} + Medium {{Solvertype = {pcm_options.get("pcm_method", "IEFPCM")} + Nonequilibrium = {pcm_options.get("neq", True)} + Solvent = {pcm_options.get("solvent", "Water")}}}""") + + if multiplicity != 1: + psi4.set_options({ + 'reference': "UHF", + 'maxiter': max_iter + 500, + 'soscf': 'true' + }) + + _, wfn = psi4.energy('SCF', return_wfn=True, molecule=mol) + psi4.core.clean() + return adcc.backends.import_scf_results(wfn) + + +def pyscf_run_pcm_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11, + conv_tol_grad=1e-10, max_iter=150, pcm_options=None): + from pyscf import scf, gto + from pyscf.solvent import ddCOSMO + + mol = gto.M( + atom=xyz, + unit="Bohr", + basis=basis, + # spin in the pyscf world is 2S + spin=multiplicity - 1, + charge=charge, + # Disable commandline argument parsing in pyscf + parse_arg=False, + dump_input=False, + verbose=0, + ) + + mf = ddCOSMO(scf.RHF(mol)) + # default eps + mf.with_solvent.eps = pcm_options.get("eps") + mf.conv_tol = conv_tol + mf.conv_tol_grad = conv_tol_grad + mf.max_cycle = max_iter + + mf.kernel() + # replace eps with eps_opt for the ADC calculation + mf.with_solvent.eps = pcm_options.get("eps_opt") + return adcc.backends.import_scf_results(mf) diff --git a/adcc/testdata/cache.py b/adcc/testdata/cache.py index 4dbfab89..95887f32 100644 --- a/adcc/testdata/cache.py +++ b/adcc/testdata/cache.py @@ -258,3 +258,5 @@ def read_yaml_data(fname): qchem_data = read_yaml_data("qchem_dump.yml") tmole_data = read_yaml_data("tmole_dump.yml") +psi4_data = read_yaml_data("psi4_dump.yml") +pyscf_data = read_yaml_data("pyscf_dump.yml") diff --git a/adcc/testdata/generate_psi4_data.py b/adcc/testdata/generate_psi4_data.py new file mode 100644 index 00000000..42de4394 --- /dev/null +++ b/adcc/testdata/generate_psi4_data.py @@ -0,0 +1,109 @@ +import psi4 +import adcc +from psi4.driver.procrouting.response.scf_response import tdscf_excitations +from adcc.testdata import static_data +import yaml +import os + + +def run_psi4_tdscf(xyz, basis, charge=0, multiplicity=1, + conv_tol=1e-12, conv_tol_grad=1e-11, max_iter=150, + pcm_options=None): + mol = psi4.geometry(f""" + {charge} {multiplicity} + {xyz} + units au + symmetry c1 + """) + psi4.set_options({ + 'basis': basis, + 'scf_type': "pK", + 'e_convergence': conv_tol, + 'd_convergence': conv_tol_grad, + 'maxiter': max_iter, + 'reference': "RHF", + 'save_jk': True, + }) + if pcm_options: + # think its better to use a dict for the pcm options, because there + # are already enough function arguments. Of course one could just + # fix the options... would be sufficient for now too. + psi4.set_options({'pcm': True, 'pcm_scf_type': "total"}) + psi4.pcm_helper(f""" + Units = AU + Cavity {{ + Type = Gepol + Area = {pcm_options.get("weight", 0.3)} + }} + Medium {{ + Solvertype = {pcm_options.get("pcm_method", "IEFPCM")} + Solvent = {pcm_options.get("solvent", "Water")} + Nonequilibrium = {pcm_options.get("neq", True)} + }} + """) + psi4.core.be_quiet() + _, wfn = psi4.energy('scf', return_wfn=True, molecule=mol) + + res = tdscf_excitations(wfn, states=5, triplets="NONE", tda=True, + r_convergence=1e-7) + + # remove cavity files from PSI4 PCM calculations + if pcm_options: + for cavityfile in os.listdir(os.getcwd()): + if cavityfile.startswith(("cavity.off_", "PEDRA.OUT_")): + os.remove(cavityfile) + return wfn, res + + +def run_adcc_ptlr(wfn): + return adcc.run_adc(wfn, method="adc1", n_singlets=5, + conv_tol=1e-7, environment="ptlr") + + +def dump_results(molecule, basis, **kwargs): + pcm_options = kwargs.get("pcm_options", None) + if pcm_options: + name = f"{molecule}_{basis}_pcm_adc1" + else: + name = f"{molecule}_{basis}_adc1" + + geom = static_data.xyz[molecule] + + ret = {"basis": basis, + "method": "adc1", + "molecule": molecule} + + wfn, res = run_psi4_tdscf(geom, basis, pcm_options=pcm_options) + + ret["energy_scf"] = wfn.energy() + # yaml safe_dump doesn't like np.floats + ret["lr_excitation_energy"] = [round(float(r["EXCITATION ENERGY"]), 9) + for r in res] + ret["lr_osc_strength"] = [round(float(r["OSCILLATOR STRENGTH (LEN)"]), 5) + for r in res] + + if pcm_options: + state = run_adcc_ptlr(wfn) + ret["ptlr_adcc_excitation_energy"] = [ + round(float(e), 9) for e in state.excitation_energy + ] + return name, ret + + +def main(): + basis_set = ["sto-3g", "cc-pvdz"] + pcm_options = {"weight": 0.3, "pcm_method": "IEFPCM", "neq": True, + "solvent": "Water"} + psi4_results = {} + for basis in basis_set: + key, ret = dump_results("formaldehyde", basis, pcm_options=pcm_options) + psi4_results[key] = ret + print(f"Dumped {key}.") + + with open("psi4_dump.yml", 'w') as yamlfile: + yaml.safe_dump(psi4_results, yamlfile, + sort_keys=False, default_flow_style=None) + + +if __name__ == "__main__": + main() diff --git a/adcc/testdata/generate_pyscf_data.py b/adcc/testdata/generate_pyscf_data.py new file mode 100644 index 00000000..83cafb6c --- /dev/null +++ b/adcc/testdata/generate_pyscf_data.py @@ -0,0 +1,94 @@ +import adcc +from pyscf import gto, scf, tdscf +from pyscf.solvent import ddCOSMO +from adcc.testdata import static_data +import yaml + + +def run_pyscf_tdscf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-12, + conv_tol_grad=1e-11, max_iter=150, pcm_options=None): + mol = gto.M( + atom=xyz, + basis=basis, + unit="Bohr", + charge=charge, + # spin in the pyscf world is 2S + spin=multiplicity - 1, + verbose=0, + # Disable commandline argument parsing in pyscf + parse_arg=False, + dump_input=False, + ) + + if pcm_options: + mf = ddCOSMO(scf.RHF(mol)) + mf.with_solvent.eps = pcm_options.get("eps") + else: + mf = scf.RHF(mol) + mf.conv_tol = conv_tol + mf.conv_tol_grad = conv_tol_grad + mf.max_cycle = max_iter + + mf.kernel() + + if pcm_options: + mf.with_solvent.eps = pcm_options.get("eps_opt") + # for n_eq solvation only PTE implemented + mf.with_solvent.equilibrium_solvation = True + cis = ddCOSMO(tdscf.TDA(mf)) + else: + cis = tdscf.TDA(mf) + cis.nstates = 5 + cis.conv_tol = 1e-7 + cis.kernel() + return mf, cis + + +def run_adcc_ptlr(wfn): + return adcc.run_adc(wfn, method="adc1", n_singlets=5, + conv_tol=1e-7, environment="ptlr") + + +def dump_results(molecule, basis, **kwargs): + pcm_options = kwargs.get("pcm_options", None) + if pcm_options: + name = f"{molecule}_{basis}_pcm_adc1" + else: + name = f"{molecule}_{basis}_adc1" + + geom = static_data.xyz[molecule] + + ret = {"basis": basis, + "method": "adc1", + "molecule": molecule} + + hf, cis = run_pyscf_tdscf(geom, basis, pcm_options=pcm_options) + + ret["energy_scf"] = float(hf.e_tot) + ret["lr_excitation_energy"] = [float(round(s, 9)) for s in cis.e] + ret["lr_osc_strength"] = [float(round(s, 5)) for s in cis.oscillator_strength()] + + if pcm_options: + state = run_adcc_ptlr(hf) + ret["ptlr_adcc_excitation_energy"] = [ + round(float(e), 9) for e in state.excitation_energy + ] + return name, ret + + +def main(): + basis_set = ["sto-3g", "cc-pvdz"] + pcm_options = {"eps": 78.3553, "eps_opt": 1.78} + pyscf_results = {} + for basis in basis_set: + key, ret = dump_results("formaldehyde", basis, pcm_options=pcm_options) + pyscf_results[key] = ret + print(f"Dumped {key}.") + + with open("pyscf_dump.yml", 'w') as yamlfile: + yaml.safe_dump(pyscf_results, yamlfile, + sort_keys=False, default_flow_style=None) + + +if __name__ == "__main__": + main() diff --git a/adcc/testdata/psi4_dump.yml b/adcc/testdata/psi4_dump.yml new file mode 100644 index 00000000..fb506bfc --- /dev/null +++ b/adcc/testdata/psi4_dump.yml @@ -0,0 +1,18 @@ +formaldehyde_sto-3g_pcm_adc1: + basis: sto-3g + method: adc1 + molecule: formaldehyde + energy_scf: -112.35458302984613 + lr_excitation_energy: [0.164670987, 0.360903879, 0.465492571, 0.50874055, 0.693857844] + lr_osc_strength: [0.0, 0.01075, 0.34344, 0.0, 0.35858] + ptlr_adcc_excitation_energy: [0.164676366, 0.360930626, 0.467575499, 0.508742077, + 0.694474733] +formaldehyde_cc-pvdz_pcm_adc1: + basis: cc-pvdz + method: adc1 + molecule: formaldehyde + energy_scf: -113.88427067069416 + lr_excitation_energy: [0.178995043, 0.38176984, 0.388019451, 0.392159876, 0.435439339] + lr_osc_strength: [0.0, 7.0e-05, 0.24871, 0.23489, 0.0] + ptlr_adcc_excitation_energy: [0.179023803, 0.381851509, 0.389582948, 0.392594025, + 0.435463213] diff --git a/adcc/testdata/pyscf_dump.yml b/adcc/testdata/pyscf_dump.yml new file mode 100644 index 00000000..5788c519 --- /dev/null +++ b/adcc/testdata/pyscf_dump.yml @@ -0,0 +1,18 @@ +formaldehyde_sto-3g_pcm_adc1: + basis: sto-3g + method: adc1 + molecule: formaldehyde + energy_scf: -112.35408506197172 + lr_excitation_energy: [0.163950251, 0.360005329, 0.465135225, 0.509146753, 0.693115282] + lr_osc_strength: [0.0, 0.01098, 0.34455, 0.0, 0.35987] + ptlr_adcc_excitation_energy: [0.163957725, 0.360034834, 0.467270151, 0.509148683, + 0.693769922] +formaldehyde_cc-pvdz_pcm_adc1: + basis: cc-pvdz + method: adc1 + molecule: formaldehyde + energy_scf: -113.88287611072602 + lr_excitation_energy: [0.177403519, 0.379894734, 0.386588779, 0.390625944, 0.435452412] + lr_osc_strength: [0.0, 0.00012, 0.24985, 0.2374, 0.0] + ptlr_adcc_excitation_energy: [0.177443175, 0.379990015, 0.388208087, 0.391106802, + 0.435484664]