ExternalPotential in Bohr (#2515)

* extern in Bohr

* docs and fix mdi

* review comments

* Update doc/sphinxman/source/scf.rst

Co-authored-by: Andy Simmonett <[email protected]>

Co-authored-by: Andy Simmonett <[email protected]>

doc/sphinxman/source/scf.rst

@@ -947,19 +947,43 @@ computations, |PSIfour| can perform more rudimentary QM/MM procedures via the
 |scf__extern| keyword. The following snippet, extracted from the
 :srcsample:`extern1` test case, demonstrates its use for a TIP3P external potential::
 
+ import numpy as np
+ external_potentials = [
+ [-0.834, np.array([1.649232019048,0.0,-2.356023604706]) / psi_bohr2angstroms],
+ [ 0.417, np.array([0.544757019107,0.0,-3.799961446760]) / psi_bohr2angstroms],
+ [ 0.417, np.array([0.544757019107,0.0,-0.912085762652]) / psi_bohr2angstroms]]
+
+ gradient('scf', external_potentials=external_potentials)
+
+The ``external_potentials`` array has three rows for three separate
+particles, and it is passed to the SCF code on the last line. The
+rows are composed of the atomic charge, x coordinate, y coordinate,
+and z coordinate in that order. The atomic charge and coordinates are
+specified in atomic units, [e] and [a0]. Add as many particle rows as
+needed to describe the full MM region.
+
+.. caution:: In |PSIfour| previous to Spring 2022 and v1.6, setting an
+ external potential like the above looked like ::
+
  Chrgfield = QMMM()
  Chrgfield.extern.addCharge(-0.834, 1.649232019048, 0.0, -2.356023604706)
  Chrgfield.extern.addCharge( 0.417, 0.544757019107, 0.0, -3.799961446760)
  Chrgfield.extern.addCharge( 0.417, 0.544757019107, 0.0, -0.912085762652)
  psi4.set_global_option_python('EXTERN', Chrgfield.extern)
 
-First a QMMM object is created, then three separate particles are added to this
-object before the SCF code is told about its existence on the last line. The
-calls to ``addCharge`` take the atomic charge, x coordinate, y coordinate, and
-z coordinate in that order. The atomic charge is specified in atomic units,
-and the coordinates always use the same units as the geometry specification in
-the regular QM region. Additional MM molecules may be specified by adding
-extra calls to ``addCharge`` to describe the full MM region.
+ gradient('scf')
+
+ The main differences are that (1) the specification of
+ charge locations in the old way used the units of the active
+ molecule, whereas the new way always uses Bohr and (2) the
+ specification of the charge and locations in the old way used the
+ :py:class:`psi4.driver.QMMM` class directly and added one charge
+ per command, whereas the new way consolidates all into an array and
+ passes it by keyword argument to the calculation.
+
+ The successor to the :py:class:`psi4.driver.QMMM` class,
+ :py:class:`psi4.driver.QMMMbohr`, is operable, but it is discouraged
+ from being used directly.
 
 To run a computation in a constant dipole field, the |scf__perturb_h|,
 |scf__perturb_with| and |scf__perturb_dipole| keywords can be used. As an

psi4/driver/__init__.py

@@ -41,7 +41,7 @@
 from psi4.driver.p4util.fcidump import *
 from psi4.driver.p4util.fchk import *
 from psi4.driver.p4util.text import *
-from psi4.driver.qmmm import QMMM
+from psi4.driver.qmmm import QMMM, QMMMbohr
 from psi4.driver.pluginutil import *
 
 from psi4.driver import gaussian_n

psi4/driver/driver.py

@@ -47,6 +47,7 @@
 from psi4.driver import driver_findif
 from psi4.driver import p4util
 from psi4.driver import qcdb
+from psi4.driver import qmmm
 from psi4.driver.procrouting import *
 from psi4.driver.p4util.exceptions import *
 from psi4.driver.mdi_engine import mdi_run
@@ -554,6 +555,10 @@ def energy(name, **kwargs):
  #for precallback in hooks['energy']['pre']:
  # precallback(lowername, **kwargs)
 
+ ep = kwargs.get('external_potentials', None)
+ if ep is not None and not isinstance(ep, dict):
+ set_external_potential(kwargs['external_potentials'])
+
  optstash = driver_util._set_convergence_criterion('energy', lowername, 6, 8, 6, 8, 6)
  optstash2 = p4util.OptionsState(['SCF', 'GUESS'])
 
@@ -733,6 +738,10 @@ def gradient(name, **kwargs):
  molecule = kwargs.pop('molecule', core.get_active_molecule())
  molecule.update_geometry()
 
+ ep = kwargs.get('external_potentials', None)
+ if ep is not None and not isinstance(ep, dict):
+ set_external_potential(kwargs['external_potentials'])
+
  # Does dertype indicate an analytic procedure both exists and is wanted?
  if dertype == 1:
  core.print_out("""gradient() will perform analytic gradient computation.\n""")
@@ -896,6 +905,10 @@ def properties(*args, **kwargs):
  if "/" in lowername:
  return driver_cbs._cbs_gufunc(properties, lowername, ptype='properties', **kwargs)
 
+ ep = kwargs.get('external_potentials', None)
+ if ep is not None and not isinstance(ep, dict):
+ set_external_potential(kwargs['external_potentials'])
+
  return_wfn = kwargs.pop('return_wfn', False)
  props = kwargs.get('properties', ['dipole', 'quadrupole'])
 
@@ -1541,6 +1554,10 @@ def hessian(name, **kwargs):
  """hessian() switching to finite difference by gradients for partial Hessian calculation.\n""")
  dertype = 1
 
+ ep = kwargs.get('external_potentials', None)
+ if ep is not None and not isinstance(ep, dict):
+ set_external_potential(kwargs['external_potentials'])
+
  # At stationary point?
  if 'ref_gradient' in kwargs:
  core.print_out("""hessian() using ref_gradient to assess stationary point.\n""")
@@ -2218,6 +2235,28 @@ def molden(wfn, filename=None, density_a=None, density_b=None, dovirtual=None):
 def tdscf(wfn, **kwargs):
  return proc.run_tdscf_excitations(wfn,**kwargs)
 
+
+def set_external_potential(external_potential):
+ """Initialize :py:class:`psi4.core.ExternalPotential` object from charges and locations.
+
+ Parameters
+ ----------
+ external_potential
+ List-like structure where each row corresponds to a charge. Lines can be composed of ``q, [x, y, z]`` or
+ ``q, x, y, z``. Locations are in [a0].
+
+ """
+ Chrgfield = qmmm.QMMMbohr()
+ for qxyz in external_potential:
+ if len(qxyz) == 2:
+ Chrgfield.extern.addCharge(qxyz[0], qxyz[1][0], qxyz[1][1], qxyz[1][2])
+ elif len(qxyz) == 4:
+ Chrgfield.extern.addCharge(qxyz[0], qxyz[1], qxyz[2], qxyz[3])
+ else:
+ raise ValidationError(f"Point charge '{qxyz}' not mapping into 'chg, [x, y, z]' or 'chg, x, y, z'")
+ core.set_global_option_python('EXTERN', Chrgfield.extern)
+
+
 # Aliases
 opt = optimize
 freq = frequency

psi4/driver/driver_nbody_helper.py

@@ -190,13 +190,11 @@ def electrostatic_embedding(metadata, pair):
  if not metadata['return_total_data']:
  raise Exception('Cannot return interaction data when using embedding scheme.')
  # Add embedding point charges
- Chrgfield = qmmm.QMMM()
+ Chrgfield = qmmm.QMMMbohr()
  for p in metadata['embedding_charges']:
  if p in pair[1]: continue
  mol = metadata['molecule'].extract_subsets([p])
  for i in range(mol.natom()):
  geom = np.array([mol.x(i), mol.y(i), mol.z(i)])
- if mol.units() == 'Angstrom':
- geom *= constants.bohr2angstroms
  Chrgfield.extern.addCharge(metadata['embedding_charges'][p][i], geom[0], geom[1], geom[2])
  core.set_global_option_python('EXTERN', Chrgfield.extern)

psi4/driver/mdi_engine.py

@@ -76,7 +76,7 @@ def __init__(self, scf_method, **kwargs):
  self.nlattice = 0 # number of lattice point charges
  self.clattice = [] # list of lattice coordinates
  self.lattice = [] # list of lattice charges
- self.lattice_field = psi4.QMMM() # Psi4 chargefield
+ self.lattice_field = psi4.QMMMbohr() # Psi4 chargefield
 
  # MPI variables
  self.mpi_world = None
@@ -281,12 +281,11 @@ def recv_masses(self, masses=None):
  def set_lattice_field(self):
  """ Set a field of lattice point charges using information received through MDI
  """
- self.lattice_field = psi4.QMMM()
- unit_conv = self.length_conversion()
+ self.lattice_field = psi4.QMMMbohr()
  for ilat in range(self.nlattice):
- latx = self.clattice[3 * ilat + 0] * unit_conv
- laty = self.clattice[3 * ilat + 1] * unit_conv
- latz = self.clattice[3 * ilat + 2] * unit_conv
+ latx = self.clattice[3 * ilat + 0]
+ laty = self.clattice[3 * ilat + 1]
+ latz = self.clattice[3 * ilat + 2]
  self.lattice_field.extern.addCharge(self.lattice[ilat], latx, laty, latz)
  psi4.core.set_global_option_python('EXTERN', self.lattice_field.extern)
  self.set_lattice = True

psi4/driver/procrouting/proc.py

@@ -1249,12 +1249,23 @@ def scf_wavefunction_factory(name, ref_wfn, reference, **kwargs):
  # For the dimer SAPT calculation, we need to account for the external potential
  # in all of the subsystems A, B, C. So we add them all in total_external_potential
  # and set the external potential to the dimer wave function
+ from psi4.driver.qmmm import QMMMbohr
+
  total_external_potential = core.ExternalPotential()
 
  for frag in kwargs['external_potentials']:
  if frag.upper() in "ABC":
- wfn.set_potential_variable(frag.upper(), kwargs['external_potentials'][frag].extern)
- total_external_potential.appendCharges(kwargs['external_potentials'][frag].extern.getCharges())
+ chrgfield = QMMMbohr()
+ for qxyz in kwargs['external_potentials'][frag]:
+ if len(qxyz) == 2:
+ chrgfield.extern.addCharge(qxyz[0], qxyz[1][0], qxyz[1][1], qxyz[1][2])
+ elif len(qxyz) == 4:
+ chrgfield.extern.addCharge(qxyz[0], qxyz[1], qxyz[2], qxyz[3])
+ else:
+ raise ValidationError(f"Point charge '{qxyz}' not mapping into 'chg, [x, y, z]' or 'chg, x, y, z'")
+
+ wfn.set_potential_variable(frag.upper(), chrgfield.extern)
+ total_external_potential.appendCharges(chrgfield.extern.getCharges())
 
  else:
  core.print_out("\n Warning! Unknown key for the external_potentials argument: %s" % frag)

psi4/driver/procrouting/sapt/fisapt_proc.py

@@ -141,14 +141,18 @@ def fisapt_fdrop(self, external_potentials=None):
  fh.write(xyz)
 
  # write external potential geometries
- if external_potentials is not None:
+ if external_potentials is not None and isinstance(external_potentials, dict):
  for frag in "ABC":
  potential = external_potentials.get(frag, None)
  if potential is not None:
- charges = potential.extern.getCharges()
- xyz = str(len(charges)) + "\n\n"
- for charge in charges:
- xyz += "Ch %f %f %f\n" % (charge[1], charge[2], charge[3])
+ xyz = str(len(potential)) + "\n\n"
+ for qxyz in potential:
+ if len(qxyz) == 2:
+ xyz += "Ch %f %f %f\n" % (qxyz[1][0], qxyz[1][1], qxyz[1][2])
+ elif len(qxyz) == 4:
+ xyz += "Ch %f %f %f\n" % (qxyz[1], qxyz[2], qxyz[3])
+ else:
+ raise ValidationError(f"Point charge '{qxyz}' not mapping into 'chg, [x, y, z]' or 'chg, x, y, z'")
 
  with open(filepath + os.sep + "Extern_%s.xyz" % frag, "w") as fh:
  fh.write(xyz)

psi4/driver/qmmm.py

@@ -120,7 +120,16 @@ def populateExtern(self, extern):
  extern.addCharge(self.molecule.Z(A), self.molecule.x(A), self.molecule.y(A), self.molecule.z(A))
 
 
-class QMMM(object):
+class QMMM():
+ """Hold charges and :py:class:`psi4.core.ExternalPotential`. Use :py:class:`psi4.driver.QMMMbohr` instead."""
+
+ def __init__(self):
+ raise UpgradeHelper(self.__class__.__name__, "QMMMbohr", 1.6, ' Replace object with a list of charges and locations in Bohr passed as keyword argument, e.g., `energy(..., embedding_charges=[[0.5, [0, 0, 1]], [-0.5, [0, 0, -1]]])`.')
+
+
+class QMMMbohr():
+ """Hold charges and :py:class:`psi4.core.ExternalPotential`. To add external charges to a calculation, prefer
+ passing the array of charges with kwarg ``external_potentials``, as in extern2 example."""
 
  def __init__(self):
  self.charges = []

psi4/src/psi4/fisapt/fisapt.cc

@@ -721,12 +721,7 @@ void FISAPT::nuclear() {
  for (int p1 = 0; p1 < pot_list.size(); p1++) {
  for (int p2 = p1+1; p2 < pot_list.size(); p2++) {
 
- bool in_angstrom = false;
- if (mol->units() == Molecule::Angstrom) {
- in_angstrom = true;
- }
-
- double IE = pot_list[p1]->computeExternExternInteraction(pot_list[p2], in_angstrom);
+ double IE = pot_list[p1]->computeExternExternInteraction(pot_list[p2]);
  // store half the interaction so that Eij + Eji = Etotal
  extern_extern_IE_matp[pot_ids[p1]][pot_ids[p2]] = IE * 0.5;
  extern_extern_IE_matp[pot_ids[p2]][pot_ids[p1]] = IE * 0.5;

psi4/src/psi4/libmints/extern.cc

@@ -67,7 +67,7 @@ void ExternalPotential::print(std::string out) const {
 
  // Charges
  if (charges_.size()) {
- printer->Printf(" > Charges [e] [molecule length units] < \n\n");
+ printer->Printf(" > Charges [e] [a0] < \n\n");
  printer->Printf(" %10s %10s %10s %10s\n", "Z", "x", "y", "z");
  for (size_t i = 0; i < charges_.size(); i++) {
  printer->Printf(" %10.5f %10.5f %10.5f %10.5f\n", std::get<0>(charges_[i]), std::get<1>(charges_[i]),
@@ -103,15 +103,12 @@ SharedMatrix ExternalPotential::computePotentialMatrix(std::shared_ptr<BasisSet>
  nthreads = Process::environment.get_n_threads();
 #endif
 
- double convfac = 1.0;
- if (basis->molecule()->units() == Molecule::Angstrom) convfac /= pc_bohr2angstroms;
-
  // Monopoles
  std::vector<std::pair<double, std::array<double, 3>>> Zxyz;
  for (size_t i=0; i< charges_.size(); ++i) {
- Zxyz.push_back({std::get<0>(charges_[i]),{{convfac * std::get<1>(charges_[i]),
- convfac * std::get<2>(charges_[i]),
- convfac * std::get<3>(charges_[i])}}});
+ Zxyz.push_back({std::get<0>(charges_[i]),{{std::get<1>(charges_[i]),
+ std::get<2>(charges_[i]),
+ std::get<3>(charges_[i])}}});
  }
 
  std::vector<SharedMatrix> V_charge;
@@ -210,14 +207,11 @@ SharedMatrix ExternalPotential::computePotentialGradients(std::shared_ptr<BasisS
  auto grad = std::make_shared<Matrix>("External Potential Gradient", natom, 3);
  double **Gp = grad->pointer();
 
- double convfac = 1.0;
- if (mol->units() == Molecule::Angstrom) convfac /= pc_bohr2angstroms;
-
  std::vector<std::pair<double, std::array<double, 3>>> Zxyz;
  for (size_t i=0; i< charges_.size(); ++i) {
- Zxyz.push_back({std::get<0>(charges_[i]),{{convfac * std::get<1>(charges_[i]),
- convfac * std::get<2>(charges_[i]),
- convfac * std::get<3>(charges_[i])}}});
+ Zxyz.push_back({std::get<0>(charges_[i]),{{std::get<1>(charges_[i]),
+ std::get<2>(charges_[i]),
+ std::get<3>(charges_[i])}}});
  }
 
  // Start with the nuclear contribution
@@ -313,9 +307,6 @@ SharedMatrix ExternalPotential::computePotentialGradients(std::shared_ptr<BasisS
 
 double ExternalPotential::computeNuclearEnergy(std::shared_ptr<Molecule> mol) {
  double E = 0.0;
- double convfac = 1.0;
-
- if (mol->units() == Molecule::Angstrom) convfac /= pc_bohr2angstroms;
 
  // Nucleus-charge interaction
  for (int A = 0; A < mol->natom(); A++) {
@@ -327,9 +318,9 @@ double ExternalPotential::computeNuclearEnergy(std::shared_ptr<Molecule> mol) {
  if (ZA > 0) { // skip Ghost interaction
  for (size_t B = 0; B < charges_.size(); B++) {
  double ZB = std::get<0>(charges_[B]);
- double xB = convfac * std::get<1>(charges_[B]);
- double yB = convfac * std::get<2>(charges_[B]);
- double zB = convfac * std::get<3>(charges_[B]);
+ double xB = std::get<1>(charges_[B]);
+ double yB = std::get<2>(charges_[B]);
+ double zB = std::get<3>(charges_[B]);
 
  double dx = xA - xB;
  double dy = yA - yB;
@@ -369,25 +360,21 @@ double ExternalPotential::computeNuclearEnergy(std::shared_ptr<Molecule> mol) {
  return E;
 }
 
-double ExternalPotential::computeExternExternInteraction(std::shared_ptr<ExternalPotential> other_extern, bool in_angstrom) {
+double ExternalPotential::computeExternExternInteraction(std::shared_ptr<ExternalPotential> other_extern) {
  double E = 0.0;
- double convfac = 1.0; // assume the geometry of the point charges is in Bohr.
- if (in_angstrom) {
- convfac /= pc_bohr2angstroms; // Here we convert to Angstrom
- }
 
  // charge-charge interaction
  for (auto self_charge: charges_) {
  double ZA = std::get<0>(self_charge);
- double xA = convfac * std::get<1>(self_charge);
- double yA = convfac * std::get<2>(self_charge);
- double zA = convfac * std::get<3>(self_charge);
+ double xA = std::get<1>(self_charge);
+ double yA = std::get<2>(self_charge);
+ double zA = std::get<3>(self_charge);
 
  for (auto other_charge: other_extern->charges_) {
  double ZB = std::get<0>(other_charge);
- double xB = convfac * std::get<1>(other_charge);
- double yB = convfac * std::get<2>(other_charge);
- double zB = convfac * std::get<3>(other_charge);
+ double xB = std::get<1>(other_charge);
+ double yB = std::get<2>(other_charge);
+ double zB = std::get<3>(other_charge);
 
  double dx = xA - xB;
  double dy = yA - yB;

psi4/src/psi4/libmints/extern.h

@@ -72,6 +72,7 @@ class PSI_API ExternalPotential {
  void setName(const std::string& name) { name_ = name; }
 
  /// Add a charge Z at (x,y,z)
+ /// Apr 2022: now always [a0] rather than units of Molecule
  void addCharge(double Z, double x, double y, double z);
 
  /// get the vector of charges
@@ -98,7 +99,7 @@ class PSI_API ExternalPotential {
  double computeNuclearEnergy(std::shared_ptr<Molecule> mol);
 
  // Compute the interaction of this potential with an external potential
- double computeExternExternInteraction(std::shared_ptr<ExternalPotential> other_extern, bool in_angstrom=false);
+ double computeExternExternInteraction(std::shared_ptr<ExternalPotential> other_extern);
 
  /// Print a trace of the external potential
  void print(std::string out_fname = "outfile") const;

pytest.ini

@@ -32,6 +32,7 @@ markers =
  df: "tests that employ density-fitting"
  dft
  dipole
+ extern: "tests that use the ExternalPotential object"
  fcidump
  freq
  fsapt

tests/extern1/CMakeLists.txt

@@ -1,3 +1,4 @@
 include(TestingMacros)
 
-add_regression_test(extern1 "psi;quicktests;scf")
+add_regression_test(extern1 "psi;quicktests;scf;extern")
+