XDM intermolecular summation

.github/workflows/pip.yml

@@ -12,25 +12,48 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        platform: [windows-latest, macos-latest, ubuntu-latest]
-        python-version: ["3.8", "3.11"]
+        python-version: ["3.10", "3.12"]
 
-    runs-on: ${{ matrix.platform }}
+    runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         submodules: true
 
-    - uses: actions/setup-python@v4
+    - uses: conda-incubator/setup-miniconda@v3
       with:
+        activate-environment: ci
+        auto-activate-base: false
+        environment-file: environment.yml
         python-version: ${{ matrix.python-version }}
 
-    - name: Add requirements
+    - name: Add Python requirements
+      shell: bash -el {0}
       run: python -m pip install --upgrade wheel setuptools
 
-    - name: Build and install
+    - name: Build with CMake
+      shell: bash -el {0}
+      run: |
+        mkdir -p build
+        cd build
+        cmake -DCMAKE_BUILD_TYPE=Release ..
+        make -j$(nproc)
+
+    - name: Install Python package
+      shell: bash -el {0}
       run: pip install --verbose .[test]
 
+    - name: Add dftd4 binary to PATH
+      shell: bash -el {0}
+      run: echo "${{ github.workspace }}/build/dftd4/app" >> $GITHUB_PATH
+
+    - name: Verify dftd4 binary
+      shell: bash -el {0}
+      run: |
+        which dftd4
+        dftd4 --version || true
+
     - name: Test
+      shell: bash -el {0}
       run: python -m pytest

.gitmodules

@@ -4,4 +4,4 @@
     branch = master
 [submodule "dftd4"]
 	path = dftd4
-	url = git@github.com:Awallace3/dftd4.git
+	url = https://github.com/Awallace3/dftd4.git

environment.yml

@@ -0,0 +1,23 @@
+name: disp
+channels:
+  - conda-forge
+  - nodefaults
+dependencies:
+    # build
+  - c-compiler
+  - cmake
+  - cxx-compiler
+  - ninja
+  - eigen                     
+  - libblas=*=*mkl
+  - libboost-headers          
+  - llvm-openmp
+  - numpy
+  - pip                      
+  - python
+  - fortran-compiler
+  - zlib
+  - setuptools
+  - pytest
+  - pytest-xdist
+  - qcelemental

setup.py

@@ -262,7 +262,7 @@ def run(self):
     },
     zip_safe=False,
     extras_require={
-        "test": ["pytest>=6.0"],
+        "test": ["pytest>=6.0", "numpy", "qcelemental"],
     },
     python_requires=">=3.8",
     classifiers=[

src/binder.cpp

@@ -109,6 +109,15 @@ PYBIND11_MODULE(dispersion, m) {
           py::arg("pos"), py::arg("carts"), py::arg("C6s"), py::arg("C8s"),
           py::arg("C10s"), py::arg("RCs"), py::arg("params"));
 
+  m_d.def("disp_2B_XDM_inter", &disp::disp_2B_XDM_inter, R"pbdoc(
+        calculate intermolecular XDM dispersion energy using dimer-level C6s, C8s,
+          C10s, RCs matrices with monA/monB indices to restrict summation to
+          intermolecular atom pairs only. params = [a1, a2]
+        )pbdoc",
+          py::arg("pos"), py::arg("carts"), py::arg("C6s"), py::arg("C8s"),
+          py::arg("C10s"), py::arg("RCs"), py::arg("monAs"), py::arg("monBs"),
+          py::arg("params"));
+
   m_d.def("disp_2B_dimer", &disp::disp_2B_dimer, R"pbdoc(
         calculate 2-body -D4 dispersion energy from positions, cartesians, C6s, and params
         for a dimer broken into two monomers with C6s and C8s

src/disp.cpp

@@ -167,6 +167,45 @@ double disp_2B_XDM_scaled(Ref<VectorXi> pos, py::EigenDRef<MatrixXd> carts,
   return energy;
 };
 
+double disp_2B_XDM_inter(Ref<VectorXi> pos, py::EigenDRef<MatrixXd> carts,
+                         py::EigenDRef<MatrixXd> C6s, py::EigenDRef<MatrixXd> C8s,
+                         py::EigenDRef<MatrixXd> C10s, py::EigenDRef<MatrixXd> RCs,
+                         Ref<VectorXi> monAs, Ref<VectorXi> monBs,
+                         Ref<VectorXd> params) {
+  // Intermolecular XDM dispersion: only sum over pairs (i in monAs, j in monBs)
+  // Uses dimer-level C6s, C8s, C10s, RCs matrices indexed by dimer atom indices
+  // params[0] = a1, params[1] = a2
+  double energy = 0;
+  double dis;
+  int A, B, i, j;
+  double a1, a2, de;
+  double ang_to_bohr = 1.8897261245650624;
+  a1 = params[0];
+  a2 = params[1] * ang_to_bohr;
+#pragma omp parallel for shared(C6s, C8s, C10s, RCs, carts, a1, a2, pos, monAs, \
+                                monBs) private(A, B, i, j, dis, de)              \
+    reduction(+ : energy)
+  for (A = 0; A < monAs.size(); A++) {
+    i = monAs[A];
+    for (B = 0; B < monBs.size(); B++) {
+      j = monBs[B];
+      double dx = carts(i, 0) - carts(j, 0);
+      double dy = carts(i, 1) - carts(j, 1);
+      double dz = carts(i, 2) - carts(j, 2);
+      double d2 = dx * dx + dy * dy + dz * dz;
+      double d6 = d2 * d2 * d2;
+      double d8 = d6 * d2;
+      double d10 = d8 * d2;
+      dis = a1 * RCs(i, j) + a2;
+
+      de = -C6s(i, j) / (d6 + pow(dis, 6)) - C8s(i, j) / (d8 + pow(dis, 8)) -
+           C10s(i, j) / (d10 + pow(dis, 10));
+      energy += de;
+    }
+  };
+  return energy;
+};
+
 double disp_2B_dimer(Ref<VectorXi> pos, py::EigenDRef<MatrixXd> carts,
                      py::EigenDRef<MatrixXd> C6s, Ref<VectorXi> pA,
                      py::EigenDRef<MatrixXd> cA, py::EigenDRef<MatrixXd> C6s_A,
@@ -1312,7 +1351,6 @@ double disp_ATM_CHG_trimer_nambe(Ref<VectorXi> pos,
         r3 = r2 * r1;
         r5 = r3 * r2;
         fdmp = 1.0 / (1.0 + 6.0 * pow(r0 / r1, alph / 3.0));
-        printf("fdmp: %f\n", fdmp);
         ang = (0.375 * (dis_ij + dis_jk - dis_ik) * (dis_ij - dis_jk + dis_ik) *
                    (-dis_ij + dis_jk + dis_ik) / r5 +
                1.0 / r3);
@@ -1386,7 +1424,6 @@ double disp_ATM_TT_trimer_nambe(Ref<VectorXi> pos,
                    (-dis_ij + dis_jk + dis_ik) / r5 +
                1.0 / r3);
         fdmp = f6_ij * f6_ik * f6_jk;
-        printf("fdmp: %f\n", fdmp);
         energy -= 6 * (ang * fdmp * c9 * triple / 6.0);
       };
     };

src/disp.hpp

@@ -36,6 +36,12 @@ double disp_2B_XDM_scaled(Ref<VectorXi> pos, py::EigenDRef<MatrixXd> carts,
                           py::EigenDRef<MatrixXd> C10s, py::EigenDRef<MatrixXd> RCs,
                           Ref<VectorXd> params);
 
+double disp_2B_XDM_inter(Ref<VectorXi> pos, py::EigenDRef<MatrixXd> carts,
+                         py::EigenDRef<MatrixXd> C6s, py::EigenDRef<MatrixXd> C8s,
+                         py::EigenDRef<MatrixXd> C10s, py::EigenDRef<MatrixXd> RCs,
+                         Ref<VectorXi> monAs, Ref<VectorXi> monBs,
+                         Ref<VectorXd> params);
+
 double disp_2B_dimer(Ref<VectorXi> pos, py::EigenDRef<MatrixXd> carts,
                      py::EigenDRef<MatrixXd> C6s, Ref<VectorXi> pA,
                      py::EigenDRef<MatrixXd> cA, py::EigenDRef<MatrixXd> C6s_A,

tests/test_basic.py

@@ -34,7 +34,7 @@ def test_disp_ATM_CHG_trimer_nambe():
     }
     geom = np.hstack((data["symbols"].reshape(-1, 1), data["geometry"]))
     qcel_trimer = qcel.models.Molecule.from_data(geom, frags=data["frags"])
-    print(qcel_trimer)
+    # print(qcel_trimer)
     C6s, C8s, C6s_ATM = pydispersion.qcel_mol_acquire_c6s(qcel_trimer)
     # Same parameters as Yi in
     # https://pubs.aip.org/aip/jcp/article/158/9/094110/2881313/Assessment-of-three-body-dispersion-models-against
@@ -46,7 +46,7 @@ def test_disp_ATM_CHG_trimer_nambe():
             "s9": 1.0,
         },
     )
-    print(f"nambe_ATM: {nambe_ATM}")
+    # print(f"nambe_ATM: {nambe_ATM}")
     return
 
 
@@ -81,7 +81,7 @@ def test_disp_ATM_TT_trimer_nambe():
     }
     geom = np.hstack((data["symbols"].reshape(-1, 1), data["geometry"]))
     qcel_trimer = qcel.models.Molecule.from_data(geom, frags=data["frags"])
-    print(qcel_trimer)
+    # print(qcel_trimer)
     C6s, C8s, C6s_ATM = pydispersion.qcel_mol_acquire_c6s(qcel_trimer)
     # Same parameters as Yi in
     # https://pubs.aip.org/aip/jcp/article/158/9/094110/2881313/Assessment-of-three-body-dispersion-models-against
@@ -93,7 +93,7 @@ def test_disp_ATM_TT_trimer_nambe():
             "s9": 1.0,
         },
     )
-    print(f"nambe_ATM: {nambe_ATM}")
+    # print(f"nambe_ATM: {nambe_ATM}")
     return
 
 if __name__ == "__main__":