Feat: Single-fragment DMET-like DFT embedding framework

examples/embedding/48-dmet-embedding.py

@@ -0,0 +1,86 @@
+# Copyright 2025 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""
+Example 48: Standard Multi-Fragment Self-Consistent DMET Calculation.
+
+This script demonstrates how to partition a molecule into multiple fragments
+and optimize the global correlation potential (u_oao) to match high-level and 
+low-level local 1-RDM density matrices self-consistently.
+"""
+
+from pyscf import gto
+from gpu4pyscf.scf import hf as gpu_hf
+from gpu4pyscf.qmmm.embedding.embedding import DMET
+
+def run_dmet_example():
+    # 1. Define the system (Ethane molecule with 6-31G basis)
+    mol = gto.Mole()
+    mol.atom = '''
+        C      -0.76091    -0.00000     0.00000
+        C       0.76091    -0.00000     0.00000
+        H      -1.16001     1.02029     0.00000
+        H      -1.16001    -0.51014    -0.88357
+        H      -1.16001    -0.51014     0.88357
+        H       1.16001    -1.02029     0.00000
+        H       1.16001     0.51014     0.88357
+        H       1.16001     0.51014    -0.88357    
+    '''
+    mol.basis = '6-31g'
+    mol.verbose = 4  # Set verbose to see detailed DMET iteration logs
+    mol.build()
+
+    print("--- Step 1: Initialize Low-Level and High-Level Solver Templates ---")
+    # In this classic exact-back-to-exact test case, we nest RHF within RHF.
+    # DMET should converge the correlation potential to exactly zero.
+    mf_outer = gpu_hf.RHF(mol)
+    mf_outer.conv_tol = 1e-12
+
+    mf_inner_template = gpu_hf.RHF(mol)
+    mf_inner_template.conv_tol = 1e-12
+
+    print("\n--- Step 2: Define Molecular Fragments ---")
+    # Partition the Ethane molecule into two methyl fragments based on atom indices:
+    # Fragment 0: First Methyl group [C1, H1, H2, H3]
+    # Fragment 1: Second Methyl group [C2, H4, H5, H6]
+    fragments = [
+        [0, 2, 3, 4],
+        [1, 5, 6, 7]
+    ]
+    print(f"Fragment 0 atom indices: {fragments[0]}")
+    print(f"Fragment 1 atom indices: {fragments[1]}")
+
+    print("\n--- Step 3: Setup and Execute the Self-Consistent DMET Solver ---")
+    dmet_solver = DMET(
+        mf_outer=mf_outer,
+        mf_inner=mf_inner_template,
+        fragments=fragments,
+        threshold=1e-5,       # SVD eigenvalue threshold for bath selection
+        max_macro_iter=20,    # Max macro loops for correlation potential fitting
+        macro_tol=1e-4        # Convergence tolerance for the density matching cost
+    )
+
+    # Trigger the DMET macroscopic self-consistent optimization
+    e_dmet = dmet_solver.kernel()
+
+    print("\n--- Final Results Summary ---")
+    # Run the raw full system RHF as an exact reference
+    e_hf_ref = mf_outer.kernel()
+
+    print(f"Global Reference RHF Energy  : {e_hf_ref:.8f} Hartree") # -79.19706462
+    print(f"Macroscopic DMET Total Energy: {e_dmet:.8f} Hartree") # -79.19706462
+    print(f"Absolute Energy Deviation    : {abs(e_dmet - e_hf_ref):.2e} Hartree") # 9.15e-11 Hartree
+
+if __name__ == '__main__':
+    run_dmet_example()

examples/embedding/49-dft-dmet-embedding.py

@@ -0,0 +1,84 @@
+# Copyright 2025 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""
+Example 49: Single-Fragment Delta-Method DFT-in-DFT Embedding.
+
+This script demonstrates how to embed a high-level hybrid DFT functional (B3LYP)
+into a localized region of a low-level GGA DFT environment (PBE) using a 
+highly-optimized projection basis without macroscopic iterations.
+"""
+
+from pyscf import gto
+from gpu4pyscf.dft import rks
+from gpu4pyscf.qmmm.embedding.embedding_dft import SingleFragmentEmbedding
+
+def run_dft_embedding_example():
+    # 1. Define the system (Ethane molecule with 6-31G basis)
+    mol = gto.Mole()
+    mol.atom = '''
+        C      -0.76091    -0.00000     0.00000
+        C       0.76091    -0.00000     0.00000
+        H      -1.16001     1.02029     0.00000
+        H      -1.16001    -0.51014    -0.88357
+        H      -1.16001    -0.51014     0.88357
+        H       1.16001    -1.02029     0.00000
+        H       1.16001     0.51014     0.88357
+        H       1.16001     0.51014    -0.88357    
+    '''
+    mol.basis = '6-31g'
+    mol.verbose = 4  # Enable to monitor localized cluster basis dimensions and logs
+    mol.build()
+
+    print("--- Step 1: Prepare Environment (PBE) and Active Region (B3LYP) Solvers ---")
+    # Low-level full system solver (Environment description)
+    mf_outer = rks.RKS(mol, xc='PBE')
+    mf_outer.conv_tol = 1e-10
+
+    # High-level solver template (Active cluster description)
+    mf_inner_template = rks.RKS(mol, xc='B3LYP')
+    mf_inner_template.conv_tol = 1e-10
+
+    print("\n--- Step 2: Define Single Target Active Fragment ---")
+    # Select only one methyl group as the active QM region. 
+    # The other half will automatically serve as the embedding environment.
+    active_fragment = [0, 2, 3, 4]
+    print(f"Target QM Active Region atom indices: {active_fragment}")
+
+    print("\n--- Step 3: Initialize and Run Single Fragment Embedding ---")
+    # Construct the single-shot embedding object. Notice that mf_inner_template 
+    # will be cloned internally via .copy() to completely avoid cache poisoning.
+    emb_obj = SingleFragmentEmbedding(
+        mf_outer=mf_outer,
+        mf_inner=mf_inner_template,
+        fragment=active_fragment,
+        threshold=1e-5  # Filters out pure fragment states and numerical noise
+    )
+
+    # Compute the final multi-scale total energy via the delta method:
+    # E_tot = E_PBE(Full) + [E_B3LYP(Active) - E_PBE(Active)]
+    e_embedded_tot = emb_obj.kernel()
+
+    print("\n--- Step 4: Verification of Template Isolation ---")
+    # Verify that our protection armor works seamlessly: 
+    # Executing the template after embedding must converge successfully without any side effects.
+    print("Verifying inner template isolation status...")
+    mf_inner_template.kernel()
+    if mf_inner_template.converged:
+        print("Template isolation check passed successfully! No cache poisoning detected.")
+    else:
+        print("Warning: Template convergence failed, check cache isolation leaks.")
+
+if __name__ == '__main__':
+    run_dft_embedding_example()

examples/embedding/50-example_ml_density_embedding.py

@@ -0,0 +1,124 @@
+# Copyright 2021-2025 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""
+Example: ML-Driven DFT Embedding (ONIOM-like scheme)
+
+This example demonstrates how to use the `HarrisRKS` and `SingleFragmentEmbedding_ML` 
+classes to perform a multi-scale quantum chemistry calculation (QM/QM). 
+It uses a dummy ML density evaluator to simulate an ultra-fast global PBE calculation, 
+and then performs a rigorous B3LYP high-level calculation only on the active fragment.
+"""
+
+import numpy as np
+import cupy as cp
+from pyscf import gto
+from gpu4pyscf.dft import rks
+from gpu4pyscf.qmmm.embedding.embedding_dft_harris import HarrisRKS, SingleFragmentEmbedding_ML
+
+
+def dummy_eval_density_func(mol, xc, grids):
+    """
+    A pure DFT surrogate that mimics the behavior of an ML density predictor.
+    It performs a standard SCF to convergence and returns the exact potentials 
+    and energies, acting as the "Ground Truth" ML model.
+    """
+    print("\n[ML Surrogate] Generating density and effective potentials...")
+    mf = rks.RKS(mol)
+    mf.xc = xc
+    mf.grids = grids
+    mf.verbose = 0
+    mf.kernel()
+
+    dm = cp.asarray(mf.make_rdm1())
+    vj, vk = mf.get_jk(mol, dm)
+    e_j = 0.5 * float(cp.sum(dm * vj))
+
+    is_hybrid = mf._numint.libxc.is_hybrid_xc(xc)
+    if is_hybrid:
+        hyb = mf._numint.libxc.hybrid_coeff(xc, spin=mol.spin)
+        vk = vk * hyb
+        e_k = 0.5 * float(cp.sum(dm * vk))
+    else:
+        vk = None
+        e_k = 0.0
+
+    _, e_xc, vxc = mf._numint.nr_rks(mol, grids, xc, dm)
+    int_rho_vxc = float(cp.sum(dm * vxc))
+
+    print("[ML Surrogate] Potential generation completed.\n")
+    return vj, vk, vxc, e_j, e_k, float(e_xc), int_rho_vxc
+
+
+def main():
+    # 1. Build a target molecule (e.g., Hexane)
+    mol = gto.Mole()
+    mol.atom = '''
+        C   1.4522500000  -2.8230000000   0.0000000000
+        C   1.4522500000  -1.2830000000   0.0000000000
+        C   0.0002500000  -0.7700000000   0.0000000000
+        C   0.0002500000   0.7700000000   0.0000000000
+        C  -1.4517500000   1.2830000000   0.0000000000
+        C  -1.4517500000   2.8230000000   0.0000000000
+        H   2.4792500000  -3.1870000000   0.0000000000
+        H   0.9382500000  -3.1870000000   0.8900000000
+        H   0.9382500000  -3.1870000000  -0.8900000000
+        H   1.9652500000  -0.9200000000   0.8900000000
+        H   1.9652500000  -0.9200000000  -0.8900000000
+        H  -0.5137500000  -1.1330000000  -0.8900000000
+        H  -0.5137500000  -1.1330000000   0.8900000000
+        H   0.5132500000   1.1330000000   0.8900000000
+        H   0.5132500000   1.1330000000  -0.8900000000
+        H  -1.9657500000   0.9200000000  -0.8900000000
+        H  -1.9657500000   0.9200000000   0.8900000000
+        H  -2.4797500000   3.1870000000   0.0000000000
+        H  -0.9377500000   3.1870000000   0.8900000000
+        H  -0.9377500000   3.1870000000  -0.8900000000
+    '''
+    mol.basis = 'sto3g' # Use a small basis set for quick demonstration
+    mol.spin = 0
+    mol.verbose = 4
+    mol.build()
+
+    # 2. Define the active region (e.g., the terminal methyl group: C + 3xH)
+    methyl_fragment = [0, 6, 7, 8]
+
+    print("==================================================")
+    print("   Starting ML-Driven DFT Embedding Calculation   ")
+    print("==================================================")
+
+    # 3. Setup the Global Low-Level Solver (driven by ML)
+    # This evaluates the full system using the Harris functional approach in 1 step.
+    mf_outer = HarrisRKS(mol, dummy_eval_density_func, xc='PBE')
+
+    # 4. Setup the Local High-Level Solver (Standard rigorous DFT)
+    # This will only be executed within the embedded active space.
+    mf_inner = rks.RKS(mol, xc='B3LYP')
+
+    # 5. Initialize and execute the ML Embedding framework
+    emb_obj = SingleFragmentEmbedding_ML(mf_outer, mf_inner, methyl_fragment)
+    e_tot = emb_obj.kernel()
+
+    print("\n==================================================")
+    print("                 Summary of Results               ")
+    print("==================================================")
+    print(f"Global Low-Level E (ML-PBE) : {mf_outer.e_tot:.8f} Hartree")
+    print(f"High-Level Local E (B3LYP)  : {emb_obj.e_inner[0]:.8f} Hartree")
+    print(f"Low-Level Local E (PBE)     : {emb_obj.e_inner[0] - emb_obj.e_tot + mf_outer.e_tot:.8f} Hartree") # Reverse engineered for display
+    print(f"--------------------------------------------------")
+    print(f"FINAL ONIOM TOTAL ENERGY    : {e_tot:.8f} Hartree")
+    print("==================================================")
+
+if __name__ == '__main__':
+    main()

gpu4pyscf/qmmm/embedding/__init__.py

@@ -0,0 +1,17 @@
+# Copyright 2025 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+from .embedding import DMET
+from .embedding_dft import SingleFragmentEmbedding