ri_ccsdt_naive: code

Author: ajz34

src/lib.rs

@@ -6,6 +6,7 @@ pub mod diis;
 pub mod rhf_slow;
 pub mod ri_rccsd;
 pub mod ri_rccsdt;
+pub mod ri_rccsdt_naive;
 pub mod ri_rccsdt_slow;
 pub mod ri_rhf;
 pub mod ri_rhf_slow;
@@ -61,6 +62,41 @@ fn playground_ri_ccsd() {
 fn playground_ri_ccsdt() {
     use crate::prelude::*;
 
+    let cint_data = CInt::from_json("assets/h2o-tzvp.json");
+    let aux_cint_data = CInt::from_json("assets/h2o-def2_jk.json");
+    let rhf_results = ri_rhf::minimal_ri_rhf(&cint_data, &aux_cint_data);
+
+    let ccsd_info = RCCSDInfo {
+        cint_data,
+        aux_cint_data,
+        mo_coeff: rhf_results.mo_coeff.clone(),
+        mo_energy: rhf_results.mo_energy.clone(),
+    };
+    let ccsd_config = CCSDConfig::default();
+
+    let (ccsd_results, ccsd_intrm) = ri_rccsd::riccsd_iteration(&ccsd_info, &ccsd_config);
+    println!("CCSD Corr Energy: {}", ccsd_results.e_corr);
+
+    println!("======");
+
+    let ccsdt_results = ri_rccsdt_naive::get_riccsd_pt_energy(&ccsd_info, &ccsd_intrm, &ccsd_results);
+    println!("CCSD(T) Perturb Energy (naive algorithm): {}", ccsdt_results.e_corr_pt);
+
+    println!("======");
+
+    let ccsdt_results = ri_rccsdt_slow::get_riccsd_pt_energy(&ccsd_info, &ccsd_intrm, &ccsd_results);
+    println!("CCSD(T) Perturb Energy (slow  algorithm): {}", ccsdt_results.e_corr_pt);
+
+    println!("======");
+
+    let ccsdt_results = ri_rccsdt::get_riccsd_pt_energy(&ccsd_info, &ccsd_intrm, &ccsd_results);
+    println!("CCSD(T) Perturb Energy (fast  algorithm): {}", ccsdt_results.e_corr_pt);
+}
+
+#[test]
+fn playground_ri_ccsdt_efficiency() {
+    use crate::prelude::*;
+
     let cint_data = CInt::from_json("assets/h2o_5-pvdz.json");
     let aux_cint_data = CInt::from_json("assets/h2o_5-pvdz_ri.json");
     let rhf_results = ri_rhf::minimal_ri_rhf(&cint_data, &aux_cint_data);
@@ -71,21 +107,18 @@ fn playground_ri_ccsdt() {
         mo_coeff: rhf_results.mo_coeff.clone(),
         mo_energy: rhf_results.mo_energy.clone(),
     };
+    let ccsd_config = CCSDConfig::default();
 
-    // do not run actual CCSD, but use faked T1 amplitudes
-    let mut ccsd_intrm = RCCSDIntermediates::default();
-    ri_rccsd::get_riccsd_intermediates_cderi(&ccsd_info, &mut ccsd_intrm);
-    let mut ccsd_results = ri_rccsd::get_riccsd_initial_guess(&ccsd_info, &ccsd_intrm);
-    ccsd_results.t1 =
-        rt::arange((ccsd_results.t1.size() as f64, ccsd_results.t1.device())).into_shape(ccsd_results.t1.shape());
+    let (ccsd_results, ccsd_intrm) = ri_rccsd::riccsd_iteration(&ccsd_info, &ccsd_config);
+    println!("CCSD Corr Energy: {}", ccsd_results.e_corr);
 
     println!("======");
 
     let ccsdt_results = ri_rccsdt_slow::get_riccsd_pt_energy(&ccsd_info, &ccsd_intrm, &ccsd_results);
-    println!("CCSD(T) Perturb Energy: {}", ccsdt_results.e_corr_pt);
+    println!("CCSD(T) Perturb Energy (slow  algorithm): {}", ccsdt_results.e_corr_pt);
 
     println!("======");
 
     let ccsdt_results = ri_rccsdt::get_riccsd_pt_energy(&ccsd_info, &ccsd_intrm, &ccsd_results);
-    println!("CCSD(T) Perturb Energy: {}", ccsdt_results.e_corr_pt);
+    println!("CCSD(T) Perturb Energy (fast  algorithm): {}", ccsdt_results.e_corr_pt);
 }

src/ri_rccsdt_naive.rs

@@ -0,0 +1,66 @@
+use crate::prelude::*;
+
+pub fn get_riccsd_pt_energy(
+    mol_info: &RCCSDInfo,
+    ccsd_intermediates: &RCCSDIntermediates,
+    ccsd_results: &RCCSDResults,
+) -> RCCSDTResults {
+    let time_outer = std::time::Instant::now();
+
+    let nocc = mol_info.nocc();
+    let nvir = mol_info.nvir();
+    let nmo = nocc + nvir;
+    let mo_energy = &mol_info.mo_energy;
+
+    let time = std::time::Instant::now();
+    let intermediates = ri_rccsdt_slow::prepare_intermediates(mol_info, ccsd_intermediates, ccsd_results);
+    println!("Time elapsed (CCSD(T) preparation): {:?}", time.elapsed());
+
+    let t1_t = intermediates.t1_t.as_ref().unwrap();
+    let t2_t = intermediates.t2_t.as_ref().unwrap();
+    let eri_vvoo_t = intermediates.eri_vvoo_t.as_ref().unwrap();
+    let eri_vooo_t = intermediates.eri_vooo_t.as_ref().unwrap();
+    let eri_vvov_t = intermediates.eri_vvov_t.as_ref().unwrap();
+
+    let device = t1_t.device().clone();
+
+    let time = std::time::Instant::now();
+
+    let wp: Tsr = unsafe { rt::empty(([nvir, nvir, nvir, nocc, nocc, nocc], &device)) };
+    (0..nvir).into_par_iter().for_each(|a| {
+        (0..nvir).into_par_iter().for_each(|b| {
+            (0..nvir).into_par_iter().for_each(|c| unsafe {
+                let mut wp = wp.force_mut();
+                let wp_1 = eri_vvov_t.i((a, b)) % t2_t.i(c).reshape((nvir, -1));
+                let wp_2 = t2_t.i((a, b)).t() % eri_vooo_t.i(c).reshape((nocc, -1));
+                wp.i_mut([a, b, c]).assign(wp_1.reshape((nocc, nocc, nocc)) - wp_2.reshape((nocc, nocc, nocc)));
+            });
+        });
+    });
+    let w = wp.transpose((0, 1, 2, 3, 4, 5))
+        + wp.transpose((0, 2, 1, 3, 5, 4))
+        + wp.transpose((1, 0, 2, 4, 3, 5))
+        + wp.transpose((1, 2, 0, 4, 5, 3))
+        + wp.transpose((2, 0, 1, 5, 3, 4))
+        + wp.transpose((2, 1, 0, 5, 4, 3));
+    let v = &w
+        + t1_t.i((.., None, None, .., None, None)) * eri_vvoo_t.i((None, .., .., None, .., ..))
+        + t1_t.i((None, .., None, None, .., None)) * eri_vvoo_t.i((.., None, .., .., None, ..))
+        + t1_t.i((None, None, .., None, None, ..)) * eri_vvoo_t.i((.., .., None, .., .., None));
+    let (so, sv) = (slice!(0, nocc), slice!(nocc, nmo));
+    let d = -mo_energy.i((sv, None, None, None, None, None))
+        - mo_energy.i((None, sv, None, None, None, None))
+        - mo_energy.i((None, None, sv, None, None, None))
+        + mo_energy.i((None, None, None, so, None, None))
+        + mo_energy.i((None, None, None, None, so, None))
+        + mo_energy.i((None, None, None, None, None, so));
+    let wt = 4.0_f64 * &w + w.transpose((1, 2, 0, 3, 4, 5)) + w.transpose((2, 0, 1, 3, 4, 5));
+    let vt = &v - v.transpose((2, 1, 0, 3, 4, 5));
+    let e_corr_pt = (wt * vt / &d).sum() / 3.0;
+
+    println!("Time elapsed (CCSD(T) computation): {:?}", time.elapsed());
+    println!("Total time elapsed (CCSD(T) energy): {:?}", time_outer.elapsed());
+    println!("Time elapsed (CCSD(T) energy): {:?}", time.elapsed());
+
+    RCCSDTResults { e_corr_pt }
+}

src/ri_rccsdt_slow.rs

@@ -97,7 +97,7 @@ fn ccsd_t_energy_contribution(abc: [usize; 3], mol_info: &RCCSDInfo, intermediat
         + get_w([c, b, a], intermediates).transpose([2, 1, 0]);
     let v = &w
         + t1_t.i((a, .., None, None)) * eri_vvoo_t.i([b, c]).i((None, .., ..))
-        + t1_t.i((b, None, .., None)) * eri_vvoo_t.i([c, a]).t().i((.., None, ..))
+        + t1_t.i((b, None, .., None)) * eri_vvoo_t.i([a, c]).i((.., None, ..))
         + t1_t.i((c, None, None, ..)) * eri_vvoo_t.i([a, b]).i((.., .., None));
     let d = -(ev[[a]] + ev[[b]] + ev[[c]]) + d_ooo;
     let z = 4.0 * &w + w.transpose([1, 2, 0]) + w.transpose([2, 0, 1])