MBIS valence charges (psi4#3226)

* MBIS valence charges Addresses psi4#3140 Adds the output to MBIS tests * Updating the glossary * Update codedeps.yaml * Update codedeps.yaml * Update azure-pipelines-windows.yml --------- Co-authored-by: Elizabeth <decolven@deshawresearch.com> Co-authored-by: Jonathon Misiewicz <jpmisiewicz@vt.edu> Co-authored-by: Lori A. Burns <lori.burns@gmail.com>
EricaCMitchell · Nov 7, 2024 · d94a68a · d94a68a
1 parent 36ddfd8
commit d94a68a
Show file tree

Hide file tree

Showing 10 changed files with 56 additions and 2 deletions.
diff --git a/.azure-pipelines/azure-pipelines-windows.yml b/.azure-pipelines/azure-pipelines-windows.yml
@@ -142,7 +142,7 @@ jobs:
                         pytest=7 ^
                         pytest-xdist ^
                         python=%PYTHON_VERSION% ^
-                        dftd3-python ^
+                        dftd3-python=0.6.0 ^
                         gcp-correction ^
                         gau2grid ^
                         libxc-c ^

diff --git a/codedeps.yaml b/codedeps.yaml
@@ -403,7 +403,7 @@ data:
         linux-64: "=*=*mkl"
         osx-64: "=*=*mkl"
         osx-arm64: "=*=*accelerate"
-        win-64: "=*=*mkl"
+        win-64: "=*=24*mkl"
       aux_build_names:
         - blas-devel
       aux_build_names_note:

diff --git a/doc/sphinxman/source/glossary_psivariables.rst b/doc/sphinxman/source/glossary_psivariables.rst
@@ -1066,6 +1066,10 @@ PSI Variables by Alpha
    Per-atom expectation value of r^3 [a0^3], equivalent to the volume
    of the MBIS-partitioned density.
 
+.. psivar:: MBIS VALENCE CHARGES
+
+   Per-atom valence charges [e] computed from an MBIS partitioned density.
+
 .. psivar:: MBIS VALENCE WIDTHS
 
    Per-atom density width [a0] of the associated valence charge computed

diff --git a/psi4/src/psi4/libmints/oeprop.cc b/psi4/src/psi4/libmints/oeprop.cc
@@ -2146,10 +2146,13 @@ std::tuple<SharedMatrix, SharedMatrix, SharedMatrix, SharedMatrix> PopulationAna
     }
 
     auto valence_widths = std::make_shared<Matrix>("MBIS Valence Widths", num_atoms, 1);
+    auto valence_charges = std::make_shared<Matrix>("MBIS Valence Charges", num_atoms, 1);
     for (int atom = 0; atom < num_atoms; atom++) {
         valence_widths->set(atom, 0, Sai[atom][mA[atom] - 1]);
+        valence_charges->set(atom, 0, -Nai[atom][mA[atom] - 1]);
     }
     wfn_->set_array_variable("MBIS VALENCE WIDTHS", valence_widths);
+    wfn_->set_array_variable("MBIS VALENCE CHARGES", valence_charges);
 
     // Compute the volume widths, only for molecules
     bool free_atom = (num_atoms == 1);
@@ -2186,6 +2189,15 @@ std::tuple<SharedMatrix, SharedMatrix, SharedMatrix, SharedMatrix> PopulationAna
             outfile->Printf("  %5d      %2s %4d   %9.6f\n", a + 1, mol->label(a).c_str(), static_cast<int>(mol->Z(a)),
                             valence_widths->get(a, 0));
         }
+
+        outfile->Printf("\n\n  MBIS Valence Charges: (a.u.)\n");
+        outfile->Printf("   Center  Symbol  Z     Charge\n");
+
+        for (int a = 0; a < num_atoms; a++) {
+            outfile->Printf("  %5d      %2s %4d   %9.6f\n", a + 1, mol->label(a).c_str(), static_cast<int>(mol->Z(a)),
+                            valence_charges->get(a, 0));
+        }
+
         if (free_atom_volumes) {
             if (free_atom == false) {
                 outfile->Printf("\n\n  MBIS Volume Ratios: \n");

diff --git a/tests/mbis-1/input.dat b/tests/mbis-1/input.dat
@@ -56,6 +56,11 @@ vwidths_ref = psi4.Matrix.from_list([ # TEST
  [0.35779695],   #TEST
  [0.35779695]])  #TEST
 
+vcharges_ref = psi4.Matrix.from_list([ # TEST
+ [-7.233675],    #TEST
+ [-0.577139],    #TEST
+ [-0.577139]])   #TEST
+
 molecule mol {
   0 1
   O
@@ -85,13 +90,15 @@ quadrupoles = wfn.array_variable('MBIS QUADRUPOLES')
 octupoles = wfn.array_variable('MBIS OCTUPOLES')
 avols = wfn.array_variable('MBIS RADIAL MOMENTS <R^3>')
 vwidths = wfn.array_variable('MBIS VALENCE WIDTHS')
+vcharges = wfn.array_variable('MBIS VALENCE CHARGES')
 
 compare_matrices(charges_ref, charges, 5, "MBIS Charges")             #TEST
 compare_matrices(dipoles_ref, dipoles, 5, "MBIS Dipoles")             #TEST
 compare_matrices(quadrupoles_ref, quadrupoles, 5, "MBIS Quadrupoles") #TEST
 compare_matrices(octupoles_ref, octupoles, 5, "MBIS Octupoles")       #TEST
 compare_matrices(avols_ref, avols, 5, "MBIS Radial Moments <r^3>")    #TEST
 compare_matrices(vwidths_ref, vwidths, 5, "MBIS Valence Widths")      #TEST
+compare_matrices(vcharges_ref, vcharges, 5, "MBIS Valence Charges")   #TEST
 
 oeprop(wfn, "MBIS_VOLUME_RATIOS", title='H20 SCF')
 

diff --git a/tests/mbis-2/input.dat b/tests/mbis-2/input.dat
@@ -60,6 +60,10 @@ vwidths_ref = np.array([  #TEST
  [0.4199382],             #TEST
  [0.4140744]])            #TEST
 
+vcharges_ref = np.array([ #TEST
+ [-7.558524],
+ [-0.808223]])
+
 vratios_ref = np.array([
 [1.568969],
 [0.457094]])
@@ -92,6 +96,7 @@ quadrupoles = wfn.variable("MBIS QUADRUPOLES")
 octupoles = wfn.variable("MBIS OCTUPOLES")
 avols = wfn.variable("MBIS RADIAL MOMENTS <R^3>")
 vwidths = wfn.variable("MBIS VALENCE WIDTHS")
+vcharges = wfn.variable("MBIS VALENCE CHARGES")
 vratios = wfn.variable("MBIS VOLUME RATIOS")
 
 compare_values(charges_ref, charges, "MBIS Charges", atol=1.e-5)             #TEST
@@ -100,4 +105,5 @@ compare_values(quadrupoles_ref, quadrupoles, "MBIS Quadrupoles", atol=1.e-5) #TE
 compare_values(octupoles_ref, octupoles, "MBIS Octupoles", atol=1.e-5)       #TEST
 compare_values(avols_ref, avols, "MBIS Radial Moments <r^3>", atol=1.e-5)    #TEST
 compare_values(vwidths_ref, vwidths, "MBIS Valence Widths", atol=1.e-5)      #TEST
+compare_values(vcharges_ref, vcharges, "MBIS Valence Charges", atol=1.e-5)    #TEST
 compare_values(vratios_ref, vratios, "MBIS Volume Ratios", atol=1.e-5)      #TEST
diff --git a/tests/mbis-3/input.dat b/tests/mbis-3/input.dat
@@ -38,6 +38,10 @@ vwidths_ref = psi4.Matrix.from_list([  #TEST
  [0.39070772],    #TEST
  [0.35579084]])   #TEST
 
+vcharges_ref = psi4.Matrix.from_list([  #TEST
+ [-6.781942],    #TEST
+ [-0.610676]])   #TEST
+
 vratios_ref = psi4.Matrix.from_list([  #TEST
 [1.139889],
 [0.211763]])
@@ -71,12 +75,14 @@ quadrupoles = wfn.array_variable('MBIS QUADRUPOLES')
 octupoles = wfn.array_variable('MBIS OCTUPOLES')
 avols = wfn.array_variable('MBIS RADIAL MOMENTS <R^3>')
 vwidths = wfn.array_variable('MBIS VALENCE WIDTHS')
+vcharges = wfn.array_variable('MBIS VALENCE CHARGES')
 vratios = wfn.array_variable('MBIS VOLUME RATIOS')
 
 compare_matrices(charges_ref, charges, 5, "MBIS Charges")             #TEST
 compare_matrices(dipoles_ref, dipoles, 5, "MBIS Dipoles")             #TEST
 compare_matrices(avols_ref, avols, 5, "MBIS Radial Moments <r^3>")    #TEST
 compare_matrices(vwidths_ref, vwidths, 5, "MBIS Valence Widths")      #TEST
+compare_matrices(vcharges_ref, vcharges, 5, "MBIS Valence Charges")      #TEST
 compare_matrices(vratios, vratios, 5, "MBIS Volume Ratios")      #TEST
 
 # lands on different degenerate solutions due to high-symmetry wfn, so test stable values. (happens for molecular quadrupole, too, not just mbis.)

diff --git a/tests/mbis-4/input.dat b/tests/mbis-4/input.dat
@@ -39,6 +39,10 @@ vwidths_ref = psi4.Matrix.from_list([  #TEST
  [1.13464706],     #TEST
  [0.56829083]])    #TEST
 
+vcharges_ref = psi4.Matrix.from_list([  #TEST
+ [-0.032768],      #TEST
+ [-9.231223]])     #TEST
+
 vratios_ref = psi4.Matrix.from_list([
 [0.088273],
 [1.490111]])
@@ -72,6 +76,7 @@ quadrupoles = wfn.array_variable('MBIS QUADRUPOLES')
 octupoles = wfn.array_variable('MBIS OCTUPOLES')
 avols = wfn.array_variable('MBIS RADIAL MOMENTS <R^3>')
 vwidths = wfn.array_variable('MBIS VALENCE WIDTHS')
+vcharges = wfn.array_variable('MBIS VALENCE CHARGES')
 vratios = wfn.array_variable('MBIS VOLUME RATIOS')
 
 
@@ -81,4 +86,5 @@ compare_matrices(quadrupoles_ref, quadrupoles, 5, "MBIS Quadrupoles") #TEST
 compare_matrices(octupoles_ref, octupoles, 5, "MBIS Octupoles")       #TEST
 compare_matrices(avols_ref, avols, 5, "MBIS Radial Moments <r^3>")    #TEST
 compare_matrices(vwidths_ref, vwidths, 5, "MBIS Valence Widths")      #TEST
+compare_matrices(vcharges_ref, vcharges, 5, "MBIS Valence Charges")   #TEST
 compare_matrices(vratios_ref, vratios, 5, "MBIS Volume Ratios")      #TEST
diff --git a/tests/mbis-5/input.dat b/tests/mbis-5/input.dat
@@ -38,6 +38,10 @@ vwidths_ref = psi4.Matrix.from_list([  #TEST
  [0.70558769],     #TEST
  [0.40965805]])    #TEST
 
+vcharges_ref = psi4.Matrix.from_list([  #TEST
+ [-2.455785],      #TEST
+ [-6.935246]])     #TEST
+
 vratios_ref = psi4.Matrix.from_list([  #TEST
 [0.652824],
 [1.379211]])
@@ -70,6 +74,7 @@ quadrupoles = wfn.array_variable('MBIS QUADRUPOLES')
 octupoles = wfn.array_variable('MBIS OCTUPOLES')
 avols = wfn.array_variable('MBIS RADIAL MOMENTS <R^3>')
 vwidths = wfn.array_variable('MBIS VALENCE WIDTHS')
+vcharges = wfn.array_variable('MBIS VALENCE CHARGES')
 vratios = wfn.array_variable('MBIS VOLUME RATIOS')
 
 compare_matrices(charges_ref, charges, 5, "MBIS Charges")             #TEST
@@ -78,4 +83,5 @@ compare_matrices(quadrupoles_ref, quadrupoles, 5, "MBIS Quadrupoles") #TEST
 compare_matrices(octupoles_ref, octupoles, 5, "MBIS Octupoles")       #TEST
 compare_matrices(avols_ref, avols, 5, "MBIS Radial Moments <r^3>")    #TEST
 compare_matrices(vwidths_ref, vwidths, 5, "MBIS Valence Widths")      #TEST
+compare_matrices(vcharges_ref, vcharges, 5, "MBIS Valence Charges")   #TEST
 compare_matrices(vratios_ref, vratios, 5, "MBIS Volume Ratios")      #TEST
diff --git a/tests/mbis-6/input.dat b/tests/mbis-6/input.dat
@@ -42,6 +42,11 @@ vwidths_ref = psi4.Matrix.from_list([ # TEST
 [0.36811050869 ],
 [0.36811050869 ]])
 
+vcharges_ref = psi4.Matrix.from_list([ # TEST
+[-7.163559],
+[-0.614615],
+[-0.614615]])
+
 molecule mol {
   0 1
   O
@@ -72,6 +77,7 @@ quadrupoles = wfn.array_variable('MBIS QUADRUPOLES')
 octupoles = wfn.array_variable('MBIS OCTUPOLES')
 avols = wfn.array_variable('MBIS RADIAL MOMENTS <R^3>')
 vwidths = wfn.array_variable('MBIS VALENCE WIDTHS')
+vcharges = wfn.array_variable('MBIS VALENCE CHARGES')
 vratios = wfn.array_variable('MBIS VOLUME RATIOS')
 
 compare_matrices(charges_ref, charges, 5, "MBIS Charges")             #TEST
@@ -80,4 +86,5 @@ compare_matrices(quadrupoles_ref, quadrupoles, 5, "MBIS Quadrupoles") #TEST
 compare_matrices(octupoles_ref, octupoles, 5, "MBIS Octupoles")       #TEST
 compare_matrices(avols_ref, avols, 5, "MBIS Radial Moments <r^3>")    #TEST
 compare_matrices(vwidths_ref, vwidths, 5, "MBIS Valence Widths")      #TEST
+compare_matrices(vcharges_ref, vcharges, 5, "MBIS Valence Charges")   #TEST
 compare_matrices(vratios_ref, vratios, 5, "MBIS Volume Ratios")      #TEST