diff --git a/src/apps/mp2/mp2.cc b/src/apps/mp2/mp2.cc index 8f7753bf0a8..b15de2ba34c 100644 --- a/src/apps/mp2/mp2.cc +++ b/src/apps/mp2/mp2.cc @@ -80,11 +80,19 @@ int main(int argc, char** argv) { if (world.rank() == 0) printf("\nstarting at time %.1fs\n", wall_time()); - const double hf_energy = mp2.get_hf().value(); - const double mp2_energy = mp2.value(); - if (world.rank() == 0) { - printf("final hf/mp2/total energy %12.8f %12.8f %12.8f\n", - hf_energy, mp2_energy, hf_energy + mp2_energy); + const double hf_energy=mp2.get_hf().value(); + const double mp2_energy=mp2.value(); +// const double mp2_energy=0.0; + if(world.rank() == 0) { + printf("final hf/mp2/total energy %12.8f %12.8f %12.8f\n", + hf_energy,mp2_energy,hf_energy+mp2_energy); + } + double mp3_correction=mp2.mp3(); + print("mp3 correction",mp3_correction); + double mp3_energy=mp3_correction+mp2_energy; + if(world.rank() == 0) { + printf("final hf/mp3/total energy %12.8f %12.8f %12.8f\n", + hf_energy,mp3_energy,hf_energy+mp3_energy); } } catch (std::exception& e) { diff --git a/src/madness/chem/CCPotentials.cc b/src/madness/chem/CCPotentials.cc index e6f76d3b160..b5de6c1b837 100644 --- a/src/madness/chem/CCPotentials.cc +++ b/src/madness/chem/CCPotentials.cc @@ -24,10 +24,12 @@ CCPotentials::CCPotentials(World& world_, std::shared_ptr nemo, const CCP //mo_ket_(make_mo_ket(nemo)), //mo_bra_(make_mo_bra(nemo)), //orbital_energies_(init_orbital_energies(nemo)) - g12(world, OT_G12, param), f12(world, OT_F12, param), +// g12(std::shared_ptrget_calc()->molecule), get_potentials(world, param), output(world) { + g12=std::shared_ptr(new CCConvolutionOperator(world,OpType::OT_G12,param)); + f12=std::shared_ptr(new CCConvolutionOperator(world,OpType::OT_F12,param)); output.debug = parameters.debug(); // reset_nemo(nemo); // g12.update_elements(mo_bra_, mo_ket_); @@ -87,7 +89,7 @@ CCPotentials::make_pair_gs(const real_function_6d& u, const CC_vecfunction& tau, CCPairFunction u_part(u); functions.push_back(u_part); if (parameters.decompose_Q()) { - CCPairFunction f_part(&f12, t(i), t(j)); + CCPairFunction f_part(f12, t(i), t(j)); functions.push_back(f_part); CCPairFunction Ot1 = apply_Ot(f_part, pt, 1); CCPairFunction Ot2 = apply_Ot(f_part, pt, 2); @@ -232,9 +234,9 @@ CCPotentials::make_pair_ex(const real_function_6d& u, const CC_vecfunction& tau, CCPairFunction u_part(u); functions.push_back(u_part); if (parameters.decompose_Q()) { - CCPairFunction f_xt(&f12, x(i), t(j)); + CCPairFunction f_xt(f12, x(i), t(j)); functions.push_back(f_xt); - CCPairFunction f_tx(&f12, t(i), x(j)); + CCPairFunction f_tx(f12, t(i), x(j)); functions.push_back(f_tx); { CCPairFunction Ot1_xt = apply_Ot(f_xt, pt, 1); // O1t(f|xt>) @@ -257,7 +259,7 @@ CCPotentials::make_pair_ex(const real_function_6d& u, const CC_vecfunction& tau, functions.push_back(QtOt_tx.invert_sign()); // - " } if (parameters.QtAnsatz()) { - CCPairFunction ftt(&f12, t(i), t(j)); // f|tt> + CCPairFunction ftt(f12, t(i), t(j)); // f|tt> CCPairFunction O1x_tt = apply_Ot(ftt, x, 1); // O1x(f|tt>) CCPairFunction OxQt_tt = apply_Qt(O1x_tt, pt, 2); // O1xQt(f|tt>) functions.push_back(OxQt_tt.invert_sign()); // - " @@ -322,11 +324,11 @@ CCPotentials::compute_pair_correlation_energy(const CCPair& u, const CC_vecfunct for (size_t mm = 0; mm < u.functions.size(); mm++) { double tmp = 0.0; - const double part1 = make_xy_op_u(mobi, mobj, g12, u.functions[mm]); + const double part1 = make_xy_op_u(mobi, mobj, *g12, u.functions[mm]); if (symmetric) tmp = part1; else //if(world.rank()==0) std::cout << std::fixed << std::setprecision(10) << part1 << "\n"; { - const double part2 = make_xy_op_u(mobj, mobi, g12, u.functions[mm]); + const double part2 = make_xy_op_u(mobj, mobi, *g12, u.functions[mm]); tmp = 2.0 * (2.0 * part1 - part2); // non symmetric pairs -> offdiagonal -> count twice } result += tmp; @@ -336,8 +338,8 @@ CCPotentials::compute_pair_correlation_energy(const CCPair& u, const CC_vecfunct } if (u.ctype == CT_CC2 && !singles.get_vecfunction().empty()) { MADNESS_ASSERT(singles.type == PARTICLE); - const double omega_s = 2.0 * mobi.inner(g12(mobj, singles(u.j)) * singles(u.i).function) - - mobi.inner(g12(mobj, singles(u.i)) * singles(u.j).function); + const double omega_s = 2.0 * mobi.inner((*g12)(mobj, singles(u.j)) * singles(u.i).function) - + mobi.inner((*g12)(mobj, singles(u.i)) * singles(u.j).function); if (world.rank() == 0) std::cout << std::setw(15) << "from singles=" << std::setfill(' ') << std::fixed << std::setprecision(10) << omega_s << "\n\n"; @@ -408,11 +410,11 @@ CCPotentials::compute_excited_pair_energy(const CCPair& d, const CC_vecfunction& const CCFunction& xbi = xbra(d.i); const CCFunction& mobj = mo_bra_(d.j); double result = 0.0; - double s2b = 2.0 * make_xy_op_u(xbi, mobj, g12, d.functions) - make_xy_op_u(mobj, xbi, g12, d.functions); + double s2b = 2.0 * make_xy_op_u(xbi, mobj, *g12, d.functions) - make_xy_op_u(mobj, xbi, *g12, d.functions); double s2c = 0.0; for (const auto& ktmp : x.functions) { const size_t k = ktmp.first; - const real_function_3d j_igk = g12(mo_bra_(d.i), mo_ket_(k)) * mo_bra_(d.j).function; + const real_function_3d j_igk = (*g12)(mo_bra_(d.i), mo_ket_(k)) * mo_bra_(d.j).function; s2c -= 2.0 * make_xy_u(xbra(k), j_igk, d.functions) - make_xy_u(j_igk, xbra(k), d.functions); } result = s2b + s2c; @@ -438,23 +440,23 @@ CCPotentials::compute_cispd_energy(const CC_vecfunction& x, const Pairs for (const auto& ktmp : x.functions) { // s2b part: const size_t k = ktmp.first; - s2b += 2.0 * make_xy_op_u(xbra(i), mo_bra_(k), g12, get_pair_function(cispd, i, k)) - - make_xy_op_u(mo_bra_(k), xbra(i), g12, get_pair_function(cispd, i, k)); - const real_function_3d kgi = g12(mo_bra_(k), mo_ket_(i)); - const real_function_3d kgxi = g12(mo_bra_(k), x(i)); - const real_function_3d kgxk = g12(mo_bra_(k), x(k)); + s2b += 2.0 * make_xy_op_u(xbra(i), mo_bra_(k), *g12, get_pair_function(cispd, i, k)) - + make_xy_op_u(mo_bra_(k), xbra(i), *g12, get_pair_function(cispd, i, k)); + const real_function_3d kgi = (*g12)(mo_bra_(k), mo_ket_(i)); + const real_function_3d kgxi = (*g12)(mo_bra_(k), x(i)); + const real_function_3d kgxk = (*g12)(mo_bra_(k), x(k)); for (const auto& ltmp : x.functions) { // s2c part: const size_t l = ltmp.first; - const real_function_3d k_lgxk = mo_bra_(k).function * g12(mo_bra_(l), x(k)); + const real_function_3d k_lgxk = mo_bra_(k).function * (*g12)(mo_bra_(l), x(k)); const real_function_3d l_kgxk = mo_bra_(l).function * kgxk; const real_function_3d l_kgi = mo_bra_(l).function * kgi; const real_function_3d l_kgxi = mo_bra_(l).function * kgxi; s2c += 2.0 * make_xy_u(xbra(i), l_kgi, get_pair_function(cispd, k, l)) - make_xy_u(l_kgi, xbra(i), get_pair_function(cispd, k, l)); const double xil = xbra(i).function.inner(x(l).function); - s4a += xil * (2.0 * make_xy_op_u(mo_bra_(l), mo_bra_(k), g12, get_pair_function(mp2, i, k)) - - make_xy_op_u(mo_bra_(k), mo_bra_(l), g12, get_pair_function(mp2, i, k))); + s4a += xil * (2.0 * make_xy_op_u(mo_bra_(l), mo_bra_(k), *g12, get_pair_function(mp2, i, k)) - + make_xy_op_u(mo_bra_(k), mo_bra_(l), *g12, get_pair_function(mp2, i, k))); s4b += 2.0 * make_xy_u(xbra(i), l_kgxi, get_pair_function(mp2, k, l)) - make_xy_u(l_kgxi, xbra(i), get_pair_function(mp2, k, l)); s4c += 4.0 * make_xy_u(xbra(i), l_kgxk, get_pair_function(mp2, i, l)) - @@ -514,7 +516,7 @@ CCPotentials::fock_residue_6d(const CCPair& u) const { // make the coulomb and local Un part with the composite factory real_function_3d hartree_potential = real_factory_3d(world); for (const auto& tmp : mo_ket_.functions) - hartree_potential += g12(mo_bra_(tmp.first), mo_ket_(tmp.first)); + hartree_potential += (*g12)(mo_bra_(tmp.first), mo_ket_(tmp.first)); real_function_3d local_part = (2.0 * hartree_potential + nemo_->ncf->U2()); if (parameters.debug()) local_part.print_size("vlocal"); @@ -853,7 +855,7 @@ CCPotentials::make_constant_part_cc2_gs(const CCPair& u, const CC_vecfunction& t CCTimer time_Vcc(world, "G(Coulomb Coupling Potential)"); real_function_6d GVcc = real_factory_6d(world); // make the g12|titj> function as op_decomposed function (not constructed in 6D) - CCPairFunction gtt(&g12, ti, tj); + CCPairFunction gtt(g12, ti, tj); // make Otau(1)(g12|titj>) CCPairFunction Otau1_gtt = apply_Ot(gtt, tau, 1); // make Otau1Q2 part and the Otau1Otau2. Otau1Otau2 part IS NOT used in the symmetry exploit @@ -920,7 +922,7 @@ CCPotentials::make_constant_part_cc2_Qt_gs(const CCPair& u, const CC_vecfunction CCTimer time_comm(world, "commutator"); const vector_real_function_3d Vtmp = get_potentials(tau, POT_singles_); const CC_vecfunction V(Vtmp, UNDEFINED, parameters.freeze()); - const CCPairFunction ftt(&f12, ti, tj); + const CCPairFunction ftt(f12, ti, tj); const CCPairFunction O1ftt = apply_Ot(ftt, V, 1); const CCPairFunction O1Q2ftt = apply_Qt(O1ftt, t, 2); const real_function_6d part1 = -2.0 * apply_G(O1Q2ftt, G); @@ -999,7 +1001,7 @@ CCPotentials::make_constant_part_cispd(const CCPair& u, const CC_vecfunction& x, real_function_6d GVcc; { time_cr.start(); - const CCPairFunction gij(&g12, moi, moj); + const CCPairFunction gij(g12, moi, moj); const CCPairFunction O1x_gij = apply_Ot(gij, x, 1); const CCPairFunction OQ_part = apply_Qt(O1x_gij, mo_ket_, 2); const real_function_6d GOQ = -2.0 * apply_G(OQ_part, G); @@ -1122,7 +1124,7 @@ CCPotentials::make_constant_part_cispd_Qt(const CCPair& u, const CC_vecfunction& const vector_real_function_3d Vxtmp = sub(world, get_potentials(x, POT_singles_), x.omega * x.get_vecfunction()); const CC_vecfunction Vx(Vxtmp, UNDEFINED, parameters.freeze()); - CCPairFunction ftt(&f12, moi, moj); + CCPairFunction ftt(f12, moi, moj); real_function_6d tmp1; real_function_6d tmp2; { @@ -1212,7 +1214,7 @@ CCPotentials::make_constant_part_cc2_ex(const CCPair& u, const CC_vecfunction& t real_function_6d tmp2; // make the xt parts of the functional and the QtOx part of the projector response { - CCPairFunction gxt(&g12, xi, tj); + CCPairFunction gxt(g12, xi, tj); // make QOtau*g*|xt> CCPairFunction O2tmp = apply_Ot(gxt, tau, 2); CCPairFunction QO = apply_Qt(O2tmp, mo_ket_, 1); @@ -1225,7 +1227,7 @@ CCPotentials::make_constant_part_cc2_ex(const CCPair& u, const CC_vecfunction& t CCPairFunction OO = apply_Ot(O1tmp, tau, 2); const real_function_6d part3 = -2.0 * apply_G(OO, G); // QtOx*g|titj> - CCPairFunction gtt(&g12, ti, tj); + CCPairFunction gtt(g12, ti, tj); CCPairFunction O2x = apply_Ot(gtt, x, 2); CCPairFunction QtOx = apply_Qt(O2x, t, 1); const real_function_6d part4 = -2.0 * apply_G(QtOx, G); @@ -1233,7 +1235,7 @@ CCPotentials::make_constant_part_cc2_ex(const CCPair& u, const CC_vecfunction& t } if (symmetric) tmp2 = swap_particles(tmp1); else { - CCPairFunction gtx(&g12, ti, xj); + CCPairFunction gtx(g12, ti, xj); // make QOtau*g*|tx> CCPairFunction O2tmp = apply_Ot(gtx, tau, 2); CCPairFunction QO = apply_Qt(O2tmp, mo_ket_, 1); @@ -1246,7 +1248,7 @@ CCPotentials::make_constant_part_cc2_ex(const CCPair& u, const CC_vecfunction& t CCPairFunction OO = apply_Ot(O1tmp, tau, 2); const real_function_6d part3 = -2.0 * apply_G(OO, G); // OxQt*g|titj> - CCPairFunction gtt(&g12, ti, tj); + CCPairFunction gtt(g12, ti, tj); CCPairFunction O1x = apply_Ot(gtt, x, 1); CCPairFunction OxQt = apply_Qt(O1x, t, 2); const real_function_6d part4 = -2.0 * apply_G(OxQt, G); @@ -1352,7 +1354,7 @@ CCPotentials::make_constant_part_cc2_Qt_ex(const CCPair& u, const CC_vecfunction const vector_real_function_3d Vtmp = get_potentials(tau, POT_singles_); const CC_vecfunction V(Vtmp, UNDEFINED, parameters.freeze()); { - const CCPairFunction fxt(&f12, xi, tj); + const CCPairFunction fxt(f12, xi, tj); const CCPairFunction O1V = apply_Ot(fxt, V, 1); const CCPairFunction OQ = apply_Qt(O1V, t, 2); const CCPairFunction O2V = apply_Ot(fxt, V, 2); @@ -1364,7 +1366,7 @@ CCPotentials::make_constant_part_cc2_Qt_ex(const CCPair& u, const CC_vecfunction real_function_6d part2; // the tx parts if (symmetric) part2 = swap_particles(part1); else { - const CCPairFunction ftx(&f12, ti, xj); + const CCPairFunction ftx(f12, ti, xj); const CCPairFunction O1V = apply_Ot(ftx, V, 1); const CCPairFunction OQ = apply_Qt(O1V, t, 2); const CCPairFunction O2V = apply_Ot(ftx, V, 2); @@ -1392,7 +1394,7 @@ CCPotentials::make_constant_part_cc2_Qt_ex(const CCPair& u, const CC_vecfunction const vector_real_function_3d Vttmp = get_potentials(tau, POT_singles_); const CC_vecfunction Vx(Vxtmp, UNDEFINED, parameters.freeze()); const CC_vecfunction Vt(Vttmp, UNDEFINED, parameters.freeze()); - CCPairFunction ftt(&f12, ti, tj); + CCPairFunction ftt(f12, ti, tj); real_function_6d tmp1; real_function_6d tmp2; { @@ -1600,7 +1602,7 @@ CCPotentials::apply_transformed_Ue(const CCFunction& x, const CCFunction& y, con const double a = inner(Uxy, tmp); const real_function_3d xx = (x.function * x.function * nemo_->ncf->square()); const real_function_3d yy = (y.function * y.function * nemo_->ncf->square()); - const real_function_3d gxx = g12(xx); + const real_function_3d gxx = (*g12)(xx); const double aa = inner(yy, gxx); const double error = std::fabs(a - aa); const double diff = a - aa; @@ -1862,7 +1864,7 @@ CCPotentials::apply_gf(const real_function_3d& f) const { BSHOperatorPtr3D(world, parameters.gamma(), parameters.lo(), parameters.thresh_poisson())); double bsh_prefactor = 4.0 * constants::pi; double prefactor = 1.0 / (2.0 * parameters.gamma()); - return prefactor * (g12(f) - bsh_prefactor * (*fBSH)(f)).truncate(); + return prefactor * ((*g12)(f) - bsh_prefactor * (*fBSH)(f)).truncate(); } double @@ -1884,11 +1886,11 @@ CCPotentials::make_xy_op_u(const CCFunction& x, const CCFunction& y, const CCCon result = inner(u.get_function(), xy_op); } else if (u.component->is_decomposed()) { if (u.component->has_operator()) { - if (op.type() == OT_G12 and u.decomposed().get_operator_ptr()->type() == OT_F12) + if (op.type() == OpType::OT_G12 and u.decomposed().get_operator_ptr()->type() == OpType::OT_F12) result = make_xy_gf_ab(x, y, u.decomposed().get_a()[0], u.decomposed().get_b()[0]); - else if (op.type() == OT_F12 and u.decomposed().get_operator_ptr()->type() == OT_G12) + else if (op.type() == OpType::OT_F12 and u.decomposed().get_operator_ptr()->type() == OpType::OT_G12) result = make_xy_gf_ab(x, y, u.decomposed().get_a()[0], u.decomposed().get_b()[0]); - else if (op.type() == OT_F12 and u.decomposed().get_operator_ptr()->type() == OT_F12) + else if (op.type() == OpType::OT_F12 and u.decomposed().get_operator_ptr()->type() == OpType::OT_F12) result = make_xy_ff_ab(x, y, u.decomposed().get_a()[0], u.decomposed().get_b()[0]); else MADNESS_EXCEPTION(("xy_" + op.name() + u.name() + " not implemented").c_str(), 1); } else { @@ -1939,9 +1941,9 @@ CCPotentials::apply_s2b_operation(const CCFunction& bra, const CCPairFunction& u real_function_3d result; MADNESS_ASSERT(particle == 1 || particle == 2); if (u.is_pure()) { - result = u.dirac_convolution(bra, g12, particle); + result = u.dirac_convolution(bra, *g12, particle); } else if (u.is_decomposed_no_op()) { - result = u.dirac_convolution(bra, g12, particle); + result = u.dirac_convolution(bra, *g12, particle); } else if (u.is_op_decomposed()) { // retunrns _particle CCFunction a; @@ -2404,7 +2406,7 @@ CCPotentials::fock_residue_closed_shell(const CC_vecfunction& singles) const { const CCFunction& taui = tmpi.second; real_function_3d hartree_potential = real_function_3d(world); for (const auto& tmpk : mo_ket_.functions) - hartree_potential += g12(mo_bra_(tmpk.first), tmpk.second); + hartree_potential += (*g12)(mo_bra_(tmpk.first), tmpk.second); const real_function_3d Ji = hartree_potential * taui.function; J.push_back(Ji); } @@ -2461,7 +2463,7 @@ madness::real_function_3d CCPotentials::K(const CCFunction& f) const { real_function_3d result = real_factory_3d(world); for (const auto k_iterator : mo_ket_.functions) { - result += g12(mo_bra_(k_iterator.first), f) * mo_ket_(k_iterator.first).function; + result += (*g12)(mo_bra_(k_iterator.first), f) * mo_ket_(k_iterator.first).function; } return result; } @@ -2487,7 +2489,7 @@ CCPotentials::apply_K(const real_function_6d& u, const size_t& particle) const { real_function_6d copyu = copy(u); real_function_6d X = (multiply(copyu, copy(mo_bra_(k).function), particle)).truncate(); // real_function_6d Y=(*poisson)(X); - real_function_6d Y = g12(X, particle); // overwrite X to save space + real_function_6d Y = (*g12)(X, particle); // overwrite X to save space result += (multiply(copy(Y), copy(mo_ket_(k).function), particle)).truncate(); // this will destroy X, but I d not intend to use it again so I choose here to save this copy } @@ -2559,7 +2561,7 @@ CCPotentials::make_f_xy(const CCFunction& x, const CCFunction& y, const real_con if (x.type != UNDEFINED && y.type != UNDEFINED) { CCTimer timer_db(world, "f|xy> sanity check"); const double test1 = (mo_bra_(y.i).function).inner(fxy.project_out(mo_bra_(x.i).function, 0)); - const double test2 = (mo_bra_(y.i).function).inner(f12(mo_bra_(x.i), x) * y.function); + const double test2 = (mo_bra_(y.i).function).inner((*f12)(mo_bra_(x.i), x) * y.function); const double sanity = test1 - test2; if (fabs(sanity) > FunctionDefaults<6>::get_thresh()) { if (world.rank() == 0) @@ -2620,11 +2622,11 @@ CCPotentials::ccs_unprojected(const CC_vecfunction& ti, const CC_vecfunction& tk for (const auto& itmp : ti.functions) { real_function_3d kgtk = real_factory_3d(world); for (const auto& ktmp : tk.functions) - kgtk += g12(mo_bra_(ktmp.first), ktmp.second); + kgtk += (*g12)(mo_bra_(ktmp.first), ktmp.second); const real_function_3d kgtk_ti = kgtk * ti(itmp.first).function; real_function_3d kgti_tk = real_factory_3d(world); for (const auto& ktmp : tk.functions) - kgti_tk += g12(mo_bra_(ktmp.first), ti(itmp.first)) * tk(ktmp.first).function; + kgti_tk += (*g12)(mo_bra_(ktmp.first), ti(itmp.first)) * tk(ktmp.first).function; const real_function_3d resulti = 2.0 * kgtk_ti - kgti_tk; result.push_back(resulti); } @@ -2654,8 +2656,8 @@ CCPotentials::x_s3a(const CC_vecfunction& x, const CC_vecfunction& t) const { for (const auto ktmp : mo_ket_.functions) { // unfrozen summation !!!!!! important !!!! const size_t k = ktmp.first; - const double gpart = make_xy_op_ab(x(i), mo_bra_(k), g12, t(i), mo_ket_(k)); - const double xpart = make_xy_op_ab(x(i), mo_bra_(k), g12, mo_ket_(k), t(i)); + const double gpart = make_xy_op_ab(x(i), mo_bra_(k), *g12, t(i), mo_ket_(k)); + const double xpart = make_xy_op_ab(x(i), mo_bra_(k), *g12, mo_ket_(k), t(i)); pot += (2.0 * gpart - xpart); } } @@ -2683,8 +2685,8 @@ CCPotentials::x_s3c(const CC_vecfunction& x, const CC_vecfunction& t) const { const size_t i = itmp.first; for (const auto ktmp : t.functions) { const size_t k = ktmp.first; - result += (2.0 * make_xy_op_ab(x(i), mo_bra_(k), g12, mo_ket_(i), t(k)) - - make_xy_op_ab(x(i), mo_bra_(k), g12, t(k), mo_ket_(i))); + result += (2.0 * make_xy_op_ab(x(i), mo_bra_(k), *g12, mo_ket_(i), t(k)) - + make_xy_op_ab(x(i), mo_bra_(k), *g12, t(k), mo_ket_(i))); } } return result; @@ -2699,8 +2701,8 @@ CCPotentials::x_s5b(const CC_vecfunction& x, const CC_vecfunction& t1, const CC_ const size_t i = itmp.first; for (const auto ktmp : t1.functions) { const size_t k = ktmp.first; - result += (2.0 * make_xy_op_ab(x(i), mo_bra_(k), g12, t1(i), t2(k)) - - make_xy_op_ab(x(i), mo_bra_(k), g12, t2(k), t1(i))); + result += (2.0 * make_xy_op_ab(x(i), mo_bra_(k), *g12, t1(i), t2(k)) - + make_xy_op_ab(x(i), mo_bra_(k), *g12, t2(k), t1(i))); } } return result; @@ -2717,8 +2719,8 @@ CCPotentials::x_s5c(const CC_vecfunction& x, const CC_vecfunction& t1, const CC_ const size_t k = ktmp.first; for (const auto ltmp : t2.functions) { const size_t l = ltmp.first; - result += (2.0 * make_xy_op_ab(mo_bra_(l), mo_bra_(k), g12, mo_ket_(i), t1(k)) - - make_xy_op_ab(mo_bra_(l), mo_bra_(k), g12, t1(k), mo_ket_(i))) * + result += (2.0 * make_xy_op_ab(mo_bra_(l), mo_bra_(k), *g12, mo_ket_(i), t1(k)) - + make_xy_op_ab(mo_bra_(l), mo_bra_(k), *g12, t1(k), mo_ket_(i))) * x(i).function.inner(t2(l).function); } } @@ -2738,8 +2740,8 @@ CCPotentials::x_s6(const CC_vecfunction& x, const CC_vecfunction& t1, const CC_v const size_t k = ktmp.first; for (const auto ltmp : t2.functions) { const size_t l = ltmp.first; - result += (2.0 * make_xy_op_ab(mo_bra_(l), mo_bra_(k), g12, t3(i), t1(k)) - - make_xy_op_ab(mo_bra_(l), mo_bra_(k), g12, t1(k), t3(i))) * + result += (2.0 * make_xy_op_ab(mo_bra_(l), mo_bra_(k), *g12, t3(i), t1(k)) - + make_xy_op_ab(mo_bra_(l), mo_bra_(k), *g12, t1(k), t3(i))) * x(i).function.inner(t2(l).function); } } @@ -2754,8 +2756,8 @@ CCPotentials::x_s2b(const CC_vecfunction& x, const Pairs& u) const { const size_t i = itmp.first; for (const auto ktmp : x.functions) { const size_t k = ktmp.first; - result += (2.0 * make_xy_op_u(x(i), mo_bra_(k), g12, get_pair_function(u, i, k)) - - make_xy_op_u(mo_bra_(k), x(i), g12, get_pair_function(u, i, k))); + result += (2.0 * make_xy_op_u(x(i), mo_bra_(k), *g12, get_pair_function(u, i, k)) - + make_xy_op_u(mo_bra_(k), x(i), *g12, get_pair_function(u, i, k))); } } return result; @@ -2768,7 +2770,7 @@ CCPotentials::x_s2c(const CC_vecfunction& x, const Pairs& u) const { const size_t i = itmp.first; for (const auto ktmp : x.functions) { const size_t k = ktmp.first; - const real_function_3d kgi = g12(mo_bra_(k), mo_ket_(i)); + const real_function_3d kgi = (*g12)(mo_bra_(k), mo_ket_(i)); for (const auto ltmp : x.functions) { const size_t l = ltmp.first; real_function_3d l_kgi = (mo_bra_(l).function * kgi).truncate(); @@ -2790,8 +2792,8 @@ CCPotentials::x_s4a(const CC_vecfunction& x, const CC_vecfunction& t, const Pair const size_t k = ktmp.first; for (const auto ltmp : x.functions) { const size_t l = ltmp.first; - result += (2.0 * make_xy_op_u(mo_bra_(l), mo_bra_(k), g12, get_pair_function(u, i, k)) - - make_xy_op_u(mo_bra_(k), mo_bra_(l), g12, get_pair_function(u, i, k))) * + result += (2.0 * make_xy_op_u(mo_bra_(l), mo_bra_(k), *g12, get_pair_function(u, i, k)) - + make_xy_op_u(mo_bra_(k), mo_bra_(l), *g12, get_pair_function(u, i, k))) * x(i).function.inner(t(l).function); } } @@ -2807,7 +2809,7 @@ CCPotentials::x_s4b(const CC_vecfunction& x, const CC_vecfunction& t, const Pair const size_t i = itmp.first; for (const auto ktmp : x.functions) { const size_t k = ktmp.first; - const real_function_3d kgti = g12(mo_bra_(k), t(i)); + const real_function_3d kgti = (*g12)(mo_bra_(k), t(i)); for (const auto ltmp : x.functions) { const size_t l = ltmp.first; real_function_3d l_kgti = (mo_bra_(l).function * kgti).truncate(); @@ -2827,10 +2829,10 @@ CCPotentials::x_s4c(const CC_vecfunction& x, const CC_vecfunction& t, const Pair const size_t i = itmp.first; for (const auto ktmp : x.functions) { const size_t k = ktmp.first; - const real_function_3d kgtk = g12(mo_bra_(k), t(k)); + const real_function_3d kgtk = (*g12)(mo_bra_(k), t(k)); for (const auto ltmp : x.functions) { const size_t l = ltmp.first; - const real_function_3d lgtk = g12(mo_bra_(l), t(k)); + const real_function_3d lgtk = (*g12)(mo_bra_(l), t(k)); const real_function_3d k_lgtk = (mo_bra_(k).function * lgtk).truncate(); const real_function_3d l_kgtk = (mo_bra_(l).function * kgtk).truncate(); result += (4.0 * make_xy_u(x(i), l_kgtk, get_pair_function(u, i, l)) - @@ -2897,7 +2899,7 @@ CCPotentials::s2c(const CC_vecfunction& singles, const Pairs& doubles) c real_function_3d resulti_r = real_factory_3d(world); for (const auto& ktmp : singles.functions) { const size_t k = ktmp.first; - const real_function_3d kgi = g12(mo_bra_(k), mo_ket_(i)); + const real_function_3d kgi = (*g12)(mo_bra_(k), mo_ket_(i)); for (const auto& ltmp : singles.functions) { const size_t l = ltmp.first; const real_function_3d l_kgi = mo_bra_(l).function * kgi; @@ -2961,7 +2963,7 @@ CCPotentials::s4b(const CC_vecfunction& singles, const Pairs& doubles) c real_function_3d resulti = real_factory_3d(world); for (const auto& ktmp : singles.functions) { const size_t k = ktmp.first; - const real_function_3d kgi = g12(mo_bra_(k), singles(i)); + const real_function_3d kgi = (*g12)(mo_bra_(k), singles(i)); vector_real_function_3d l_kgi = mul_sparse(world, kgi, active_mo_bra, parameters.thresh_3D()); truncate(world, l_kgi); for (const auto& ltmp : singles.functions) { @@ -2992,7 +2994,7 @@ CCPotentials::s4c(const CC_vecfunction& singles, const Pairs& doubles) c real_function_3d kgtauk = real_factory_3d(world); for (const auto& ktmp : singles.functions) { const size_t k = ktmp.first; - kgtauk += g12(mo_bra_(k), singles(k)); + kgtauk += (*g12)(mo_bra_(k), singles(k)); } vector_real_function_3d l_kgtauk = mul(world, kgtauk, active_mo_bra); truncate(world, l_kgtauk); @@ -3005,7 +3007,7 @@ CCPotentials::s4c(const CC_vecfunction& singles, const Pairs& doubles) c } for (const auto& ktmp : singles.functions) { const size_t k = ktmp.first; - const real_function_3d k_lgtauk = (mo_bra_(k).function * g12(mo_bra_(l), singles(k))).truncate(); + const real_function_3d k_lgtauk = (mo_bra_(k).function * (*g12)(mo_bra_(l), singles(k))).truncate(); for (size_t mm = 0; mm < uil.size(); mm++) { part3 += uil[mm].project_out(k_lgtauk, 2); part4 += uil[mm].project_out(k_lgtauk, 1); @@ -3104,7 +3106,7 @@ void CCPotentials::test_pair_consistency(const CCPairFunction& u, const size_t i v1.push_back(u); std::vector v2; v2.push_back(u); - CCPairFunction ftt(&f12, mo_ket_(i), mo_ket_(j)); + CCPairFunction ftt(f12, mo_ket_(i), mo_ket_(j)); CCPairFunction O1xftt = apply_Ot(ftt, x, 1); CCPairFunction OxQftt = apply_Qt(O1xftt, mo_ket_, 2); CCPairFunction OxQ = OxQftt.invert_sign(); @@ -3129,7 +3131,7 @@ void CCPotentials::test_pair_consistency(const CCPairFunction& u, const size_t i v1.push_back(u); std::vector v2; v2.push_back(u); - CCPairFunction ftt(&f12, mo_ket_(i), mo_ket_(j)); + CCPairFunction ftt(f12, mo_ket_(i), mo_ket_(j)); CCPairFunction O1xftt = apply_Ot(ftt, x, 1); CCPairFunction OxQftt = apply_Qt(O1xftt, mo_ket_, 2); v2.push_back(OxQftt); @@ -3165,8 +3167,8 @@ bool CCPotentials::test_compare_pairs(const CCPair& pair1, const CCPair& pair2) } else output("Test Passed, diff=" + std::to_string(diff)); // test energy integration - double energy_1 = make_xy_op_u(mo_bra_(pair1.i), mo_bra_(pair1.j), g12, pair1.functions); - double energy_2 = make_xy_op_u(mo_bra_(pair2.i), mo_bra_(pair2.j), g12, pair2.functions); + double energy_1 = make_xy_op_u(mo_bra_(pair1.i), mo_bra_(pair1.j), *g12, pair1.functions); + double energy_2 = make_xy_op_u(mo_bra_(pair2.i), mo_bra_(pair2.j), *g12, pair2.functions); double diff_energy = energy_1 - energy_2; if (world.rank() == 0) std::cout << std::fixed << std::setprecision(10) @@ -3207,7 +3209,7 @@ real_function_6d CCPotentials::make_6D_pair(const CCPair& pair) const { result += ab; } } else if (f.is_op_decomposed()) { - MADNESS_ASSERT(f.get_operator().type() == OT_F12); + MADNESS_ASSERT(f.get_operator().type() == OpType::OT_F12); real_function_6d fxy = make_f_xy(f.get_a()[0], f.get_b()[0]); result += fxy; } else MADNESS_EXCEPTION("Unknown type of CCPairFunction", 1); @@ -3412,13 +3414,11 @@ void CCPotentials::test_singles_potential() const { void CCPotentials::test() { output.section("Testing enums"); - PairFormat test1 = PT_DECOMPOSED; CalcType test2 = CT_MP2; - OpType test3 = OT_G12; + OpType test3 = OpType::OT_G12; FuncType test4 = HOLE; CCState test5 = GROUND_STATE; PotentialType test6 = POT_F3D_; - assign_name(test1); assign_name(test2); assign_name(test3); assign_name(test4); @@ -3490,7 +3490,7 @@ void CCPotentials::test() { const double lo = parameters.thresh_6D(); const double hi = parameters.thresh_3D(); { - CCPairFunction fab(&f12, a, b); + CCPairFunction fab(f12, a, b); const double test1 = overlap(fab, fab); const double prefactor = 1.0 / (4 * y * y); const double test2 = prefactor * (aR.inner(a) * bR.inner(b) - 2.0 * bb.inner(af2a) + bb.inner(affa)); @@ -3563,7 +3563,7 @@ void CCPotentials::test() { if (fabs(diff) > hi) passed_hi = false; } { - CCPairFunction fab(&f12, a, b); + CCPairFunction fab(f12, a, b); CCPairFunction ab2(ab_6d); const double test1 = overlap(fab, ab2); const double test2 = bb.inner(afa); @@ -3575,7 +3575,7 @@ void CCPotentials::test() { if (fabs(diff) > hi) passed_hi = false; } { - CCPairFunction fab(&f12, a, b); + CCPairFunction fab(f12, a, b); CCPairFunction ab2(mo_ket_.get_vecfunction(), mo_ket_.get_vecfunction()); const double test1 = overlap(fab, ab2); const double test2 = bb.inner(afa); diff --git a/src/madness/chem/CCPotentials.h b/src/madness/chem/CCPotentials.h index 6258aaf535f..18d4c227cd3 100644 --- a/src/madness/chem/CCPotentials.h +++ b/src/madness/chem/CCPotentials.h @@ -794,9 +794,9 @@ class CCPotentials { // update the intermediates void update_intermediates(const CC_vecfunction& t) { - g12.update_elements(mo_bra_, t); + g12->update_elements(mo_bra_, t); // g12.sanity(); - f12.update_elements(mo_bra_, t); + f12->update_elements(mo_bra_, t); // f12.sanity(); } @@ -829,9 +829,9 @@ class CCPotentials { std::vector orbital_energies_; /// the coulomb operator with all intermediates public: - CCConvolutionOperator g12; + std::shared_ptr g12; /// the f12 operator with all intermediates - CCConvolutionOperator f12; + std::shared_ptr f12; /// the correlation factor, holds necessary regularized potentials CorrelationFactor corrfac; /// Manager for stored intermediate potentials which are s2c, s2b and the whole singles potentials without fock-residue for GS and EX state diff --git a/src/madness/chem/CCStructures.cc b/src/madness/chem/CCStructures.cc index 0ba3e24d52b..998ad8cd05a 100644 --- a/src/madness/chem/CCStructures.cc +++ b/src/madness/chem/CCStructures.cc @@ -322,7 +322,7 @@ CCConvolutionOperator::operator()(const CCFunction& bra, const CCFunction& ket, real_function_6d CCConvolutionOperator::operator()(const real_function_6d& u, const size_t particle) const { MADNESS_ASSERT(particle == 1 or particle == 2); - MADNESS_ASSERT(operator_type == OT_G12); + MADNESS_ASSERT(operator_type == OpType::OT_G12); MADNESS_ASSERT(op); op->particle() = particle; return (*op)(u); @@ -331,7 +331,7 @@ real_function_6d CCConvolutionOperator::operator()(const real_function_6d& u, co real_function_3d CCConvolutionOperator::operator()(const CCFunction& bra, const real_function_6d& u, const size_t particle) const { MADNESS_ASSERT(particle == 1 or particle == 2); - MADNESS_ASSERT(operator_type == OT_G12); + MADNESS_ASSERT(operator_type == OpType::OT_G12); MADNESS_ASSERT(op); const real_function_6d tmp = multiply(copy(u), copy(bra.function), particle); op->particle() = particle; @@ -406,33 +406,47 @@ size_t CCConvolutionOperator::info() const { SeparatedConvolution * CCConvolutionOperator::init_op(const OpType& type, const Parameters& parameters) const { + bool debug=false; + bool printme=(world.rank()==0) and debug; switch (type) { - case OT_G12 : { - if (world.rank() == 0) + case OpType::OT_G12 : { + if (printme) std::cout << "Creating " << assign_name(type) << " Operator with thresh=" << parameters.thresh_op << " and lo=" << parameters.lo << std::endl; return CoulombOperatorPtr(world, parameters.lo, parameters.thresh_op); } - case OT_F12 : { - if (world.rank() == 0) + case OpType::OT_F12 : { + if (printme) std::cout << "Creating " << assign_name(type) << " Operator with thresh=" << parameters.thresh_op << " and lo=" << parameters.lo << " and Gamma=" << parameters.gamma << std::endl; return SlaterF12OperatorPtr(world, parameters.gamma, parameters.lo, parameters.thresh_op); } - case OT_SLATER : { - if (world.rank() == 0) + case OpType::OT_SLATER : { + if (printme) std::cout << "Creating " << assign_name(type) << " Operator with thresh=" << parameters.thresh_op << " and lo=" << parameters.lo << " and Gamma=" << parameters.gamma << std::endl; return SlaterOperatorPtr(world, parameters.gamma, parameters.lo, parameters.thresh_op); } - case OT_BSH : { - if (world.rank() == 0) + case OpType::OT_BSH : { + if (printme) std::cout << "Creating " << assign_name(type) << " Operator with thresh=" << parameters.thresh_op << " and lo=" << parameters.lo << " and Gamma=" << parameters.gamma << std::endl; return BSHOperatorPtr3D(world, parameters.gamma, parameters.lo, parameters.thresh_op); } - case OT_ONE : { - if (world.rank() == 0) + case OpType::OT_FG12: { + if (printme) + std::cout << "Creating " << assign_name(type) << " Operator with thresh=" << parameters.thresh_op + << " and lo=" << parameters.lo << " and Gamma=" << parameters.gamma << std::endl; + return FGOperatorPtr(world, parameters.gamma, parameters.lo, parameters.thresh_op); + } + case OpType::OT_F2G12: { + if (printme) + std::cout << "Creating " << assign_name(type) << " Operator with thresh=" << parameters.thresh_op + << " and lo=" << parameters.lo << " and Gamma=" << parameters.gamma << std::endl; + return F2GOperatorPtr(world, parameters.gamma, parameters.lo, parameters.thresh_op); + } + case OpType::OT_ONE : { + if (printme) std::cout << "Creating " << assign_name(type) << " Operator " << std::endl; return nullptr; } @@ -444,25 +458,6 @@ CCConvolutionOperator::init_op(const OpType& type, const Parameters& parameters) } -/// Assigns strings to enums for formated output -std::string -assign_name(const PairFormat& input) { - switch (input) { - case PT_FULL: - return "full"; - case PT_DECOMPOSED: - return "decomposed"; - case PT_OP_DECOMPOSED: - return "operator-decomposed"; - default: { - MADNESS_EXCEPTION("Unvalid enum assignement!", 1); - return "undefined"; - } - } - MADNESS_EXCEPTION("assign_name:pairtype, should not end up here", 1); - return "unknown pairtype"; -} - /// Assigns strings to enums for formated output std::string assign_name(const CCState& input) { @@ -484,17 +479,17 @@ assign_name(const CCState& input) { std::string assign_name(const OpType& input) { switch (input) { - case OT_G12: + case OpType::OT_G12: return "g12"; - case OT_F12: + case OpType::OT_F12: return "f12"; - case OT_SLATER: + case OpType::OT_SLATER: return "slater"; - case OT_FG12: + case OpType::OT_FG12: return "fg12"; - case OT_BSH: + case OpType::OT_BSH: return "bsh"; - case OT_ONE: + case OpType::OT_ONE: return "identity"; default: { MADNESS_EXCEPTION("Unvalid enum assignement!", 1); diff --git a/src/madness/chem/CCStructures.h b/src/madness/chem/CCStructures.h index 99470bf0e43..3189b9c4c98 100644 --- a/src/madness/chem/CCStructures.h +++ b/src/madness/chem/CCStructures.h @@ -21,18 +21,6 @@ namespace madness { -/// Operatortypes used by the CCConvolutionOperator Class -enum OpType { - OT_UNDEFINED, - OT_ONE, /// indicates the identity - OT_G12, /// 1/r - OT_SLATER, /// exp(r) - OT_F12, /// 1-exp(r) - OT_FG12, /// (1-exp(r))/r - OT_F212, /// (1-exp(r))^2 - OT_BSH /// exp(r)/r -}; - /// Calculation Types used by CC2 enum CalcType { CT_UNDEFINED, CT_MP2, CT_CC2, CT_LRCCS, CT_LRCC2, CT_CISPD, CT_ADC2, CT_TDHF, CT_TEST @@ -62,10 +50,6 @@ enum PotentialType { POT_singles_ }; -/// Assigns strings to enums for formated output -std::string -assign_name(const PairFormat& input); - /// Assigns strings to enums for formated output std::string assign_name(const CCState& input); @@ -744,66 +728,25 @@ struct CCConvolutionOperator { CCConvolutionOperator(const CCConvolutionOperator& other) = default; - friend bool can_combine(const CCConvolutionOperator& left, const CCConvolutionOperator& right) { - return (combine_OT(left,right).first!=OT_UNDEFINED); +protected: + + friend CCConvolutionOperator combine(const CCConvolutionOperator& a, const CCConvolutionOperator& b) { + auto info= SeparatedConvolution::combine_OT((*a.get_op()),(*b.get_op())); + Parameters param; + param.gamma=info.mu; + param.thresh_op=info.thresh; + param.lo=info.lo; + param.freeze=a.parameters.freeze; + return CCConvolutionOperator(a.world, info.type, param); } - friend std::pair combine_OT(const CCConvolutionOperator& left, const CCConvolutionOperator& right) { - OpType type=OT_UNDEFINED; - Parameters param=left.parameters; - if ((left.type()==OT_F12) and (right.type()==OT_G12)) { - type=OT_FG12; - } - if ((left.type()==OT_G12) and (right.type()==OT_F12)) { - type=OT_FG12; - param.gamma=right.parameters.gamma; - } - if ((left.type()==OT_F12) and (right.type()==OT_F12)) { - type=OT_F212; - // keep the original gamma - // (f12)^2 = (1- slater12)^2 = 1/(4 gamma) (1 - 2 exp(-gamma) + exp(-2 gamma)) - MADNESS_CHECK(right.parameters.gamma == left.parameters.gamma); - } - return std::make_pair(type,param); - } - - - /// combine 2 convolution operators to one - - /// @return a vector of pairs: factor and convolution operator - friend std::vector> combine(const CCConvolutionOperator& left, const CCConvolutionOperator& right) { - MADNESS_CHECK(can_combine(left,right)); - MADNESS_CHECK(left.world.id()==right.world.id()); - auto [type,param]=combine_OT(left,right); - std::vector> result; - if (type==OT_FG12) { - // fg = (1 - exp(-gamma r12)) / r12 = 1/r12 - exp(-gamma r12)/r12 = coulomb - bsh - - // coulombfit return 1/r - // we need 1/(2 gamma) 1/r - result.push_back(std::make_pair(1.0/(2.0*param.gamma),CCConvolutionOperator(left.world, OT_G12, param))); - - // bshfit returns 1/(4 pi) exp(-gamma r)/r - // we need 1/(2 gamma) exp(-gamma r)/r - const double factor = 4.0 * constants::pi /(2.0*param.gamma); - result.push_back(std::make_pair(-factor,CCConvolutionOperator(left.world, OT_BSH, param))); - } else if (type==OT_F212) { -// we use the slater operator which is S = e^(-y*r12), y=gamma -// the f12 operator is: 1/2y*(1-e^(-y*r12)) = 1/2y*(1-S) -// so the squared f12 operator is: f*f = 1/(4*y*y)(1-2S+S*S), S*S = S(2y) = e(-2y*r12) -// we have then: = 1/(4*y*y)*( - 2* + ) -// we have then: =( - 1/(4*y*y)*2* - MADNESS_CHECK(left.parameters.gamma==right.parameters.gamma); - const double prefactor = 1.0 / (4.0 * param.gamma); // Slater has no 1/(2 gamma) per se. - Parameters param2=param; - param2.gamma*=2.0; - result.push_back(std::make_pair(1.0*prefactor,CCConvolutionOperator(left.world, OT_ONE, param))); - result.push_back(std::make_pair(-2.0*prefactor,CCConvolutionOperator(left.world, OT_SLATER, left.parameters))); - result.push_back(std::make_pair(1.0*prefactor,CCConvolutionOperator(left.world, OT_SLATER, param2))); - } - return result; + friend std::shared_ptr combine(const std::shared_ptr& a, + const std::shared_ptr& b) { + auto bla=combine(*a,*b); + return std::shared_ptr(new CCConvolutionOperator(combine(*a,*b))); } +public: /// @param[in] f: a 3D function /// @param[out] the convolution op(f), no intermediates are used real_function_3d operator()(const real_function_3d& f) const { @@ -834,8 +777,8 @@ struct CCConvolutionOperator { // @param[in] f: a vector of 3D functions // @param[out] the convolution of op with each function, no intermeditates are used vector_real_function_3d operator()(const vector_real_function_3d& f) const { - if (op) return apply(world, (*op), f); - return f; + MADNESS_CHECK(op); + return apply(world, (*op), f); } // @param[in] bra: a 3D CC_function, if nuclear-correlation factors are used they have to be applied before @@ -891,9 +834,10 @@ struct CCConvolutionOperator { /// create a TwoElectronFactory with the operatorkernel TwoElectronFactory get_kernel() const { - if (type() == OT_G12) return TwoElectronFactory(world).dcut(1.e-7); - else if (type() == OT_F12) return TwoElectronFactory(world).dcut(1.e-7).f12().gamma(parameters.gamma); - else if (type() == OT_FG12) return TwoElectronFactory(world).dcut(1.e-7).BSH().gamma(parameters.gamma); + if (type() == OpType::OT_G12) return TwoElectronFactory(world).dcut(1.e-7); + else if (type() == OpType::OT_F12) return TwoElectronFactory(world).dcut(1.e-7).f12().gamma(parameters.gamma); + else if (type() == OpType::OT_FG12) return TwoElectronFactory(world).dcut(1.e-7).BSH().gamma(parameters.gamma); + else if (type() == OpType::OT_SLATER) return TwoElectronFactory(world).dcut(1.e-7).slater().gamma(parameters.gamma); else error("no kernel of type " + name() + " implemented"); return TwoElectronFactory(world); } @@ -908,7 +852,7 @@ struct CCConvolutionOperator { /// the world World& world; /// the operatortype, currently this can be g12_ or f12_ - const OpType operator_type = OT_UNDEFINED; + const OpType operator_type = OpType::OT_UNDEFINED; /// @param[in] optype: can be f12_ or g12_ depending on which operator shall be intitialzied /// @param[in] parameters: parameters (thresholds etc) @@ -928,8 +872,16 @@ struct CCConvolutionOperator { << "!!!!!\n\n" << std::endl; MADNESS_EXCEPTION(msg.c_str(), 1); } +public: }; +template +std::shared_ptr CCConvolutionOperatorPtr(World& world, const OpType type, + CCConvolutionOperator::Parameters param) { + return std::shared_ptr(new CCConvolutionOperator(world,type,param)); +} + + class CCPair : public archive::ParallelSerializableObject { public: CCPair(){}; diff --git a/src/madness/chem/CMakeLists.txt b/src/madness/chem/CMakeLists.txt index 345216e7cc3..9cb9cc0298c 100644 --- a/src/madness/chem/CMakeLists.txt +++ b/src/madness/chem/CMakeLists.txt @@ -24,6 +24,7 @@ set(MADCHEM_HEADERS gth_pseudopotential.h GuessFactory.h localizer.h + lowrankfunction.h masks_and_boxes.h molecularbasis.h molecular_functors.h @@ -139,8 +140,8 @@ if(BUILD_TESTING) SET(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) # The list of unit test source files set(CHEM_TEST_SOURCES_SHORT test_pointgroupsymmetry.cc test_masks_and_boxes.cc - test_qc.cc test_MolecularOrbitals.cc test_BSHApply.cc) - set(CHEM_TEST_SOURCES_LONG test_localizer.cc test_ccpairfunction.cc) + test_qc.cc test_MolecularOrbitals.cc test_BSHApply.cc) + set(CHEM_TEST_SOURCES_LONG test_localizer.cc test_ccpairfunction.cc test_low_rank_function.cc) if (LIBXC_FOUND) list(APPEND CHEM_TEST_SOURCES_SHORT test_dft.cc ) list(APPEND CHEM_TEST_SOURCES_LONG test_SCFOperators.cc) diff --git a/src/madness/chem/PNOF12Potentials.cpp b/src/madness/chem/PNOF12Potentials.cpp index ce3b5659386..6b5e5092d65 100644 --- a/src/madness/chem/PNOF12Potentials.cpp +++ b/src/madness/chem/PNOF12Potentials.cpp @@ -81,8 +81,6 @@ F12Potentials::F12Potentials(World& world,const Nemo& nemo, const BasisFunctions fop = std::shared_ptr < real_convolution_3d > (SlaterF12OperatorPtr(world, param.gamma(), lo, eps)); slaterop = std::shared_ptr < real_convolution_3d > (SlaterF12OperatorPtr(world, param.gamma(), lo, eps)); slaterop_sq = std::shared_ptr < real_convolution_3d > (SlaterF12OperatorPtr(world, param.gamma() * 2.0, lo, eps)); - slaterop->is_slaterf12 = false; - slaterop_sq->is_slaterf12 = false; fop = std::shared_ptr < real_convolution_3d > (SlaterF12OperatorPtr(world, param.gamma(), lo, eps)); //test diff --git a/src/madness/chem/TDHF.cc b/src/madness/chem/TDHF.cc index 764896ece0b..a21e282ebf4 100644 --- a/src/madness/chem/TDHF.cc +++ b/src/madness/chem/TDHF.cc @@ -69,7 +69,7 @@ void TDHF::initialize() { msg.section("Initialize TDHF Class"); msg.debug = parameters.debug(); - g12=std::make_shared(world, OT_G12, parameters.get_ccc_parameters(get_calcparam().lo())); + g12=std::make_shared(world, OpType::OT_G12, parameters.get_ccc_parameters(get_calcparam().lo())); const double old_thresh = FunctionDefaults<3>::get_thresh(); if (old_thresh > parameters.thresh() * 0.1 and old_thresh > 1.e-5) { diff --git a/src/madness/chem/ccpairfunction.cc b/src/madness/chem/ccpairfunction.cc index 46bcd016907..2c645534bff 100644 --- a/src/madness/chem/ccpairfunction.cc +++ b/src/madness/chem/ccpairfunction.cc @@ -5,6 +5,7 @@ #include #include #include +#include using namespace madness; @@ -17,6 +18,99 @@ CCPairFunction::invert_sign() { return *this; } +bool CCPairFunction::is_convertible_to_pure_no_op() const { + if (has_operator()) { + const auto type=get_operator().type(); + if (not (type==OpType::OT_SLATER or type==OpType::OT_F12)) return false; + } + if (is_decomposed() and (get_a().size()>1)) return false; + return true; +}; + + +void CCPairFunction::convert_to_pure_no_op_inplace() { + pureT tmp; + pureT result=real_factory_6d(world()); + if (is_pure_no_op()) { + return; + } else if (is_pure()) { + tmp= CompositeFactory(world()) + .g12(get_operator().get_kernel()) + .ket(get_function()); + tmp.fill_tree(); + result=tmp; + } else if (is_decomposed()) { + MADNESS_CHECK(get_a().size()<3); + for (int i=0; i(world()) + .g12(get_operator().get_kernel()) + .particle1(get_a()[i]) + .particle2(get_b()[i]); + } else if (is_decomposed_no_op()) { + tmp= CompositeFactory(world()) + .particle1(get_a()[i]) + .particle2(get_b()[i]); + } + tmp.fill_tree(); + result+=tmp; + } + } + component.reset(new TwoBodyFunctionPureComponent(result)); +}; + +std::vector consolidate(const std::vector& other) { + + std::vector result; + std::vector all_pure; + for (auto& c : other) { + if (c.is_pure_no_op()) all_pure.push_back(c.get_function()); + else result.push_back(c); + } + if (not all_pure.empty()) { + for (std::size_t i=1; i& v1) { + auto a012=std::array{0,1,2}; + auto a345=std::array{3,4,5}; + int particle=-1; + if (v1== a012) particle=0; + if (v1== a345) particle=1; + MADNESS_CHECK(particle==0 or particle==1); + World& world=other.world(); + + if (other.is_decomposed()) { + if (particle == 0) { + return CCPairFunction(other.get_operator_ptr(), f * other.get_a(), copy(world, other.get_b())); + } else { + return CCPairFunction(other.get_operator_ptr(), copy(world, other.get_a()), f * other.get_b()); + } + } else if (other.is_pure()) { + auto tmp=multiply(other.get_function(),f,particle+1); + return CCPairFunction(other.get_operator_ptr(),tmp); + } else { + MADNESS_EXCEPTION("confused CCPairFunction in multiply",1); + } +}; + + +/// multiplication with a 2-particle function +CCPairFunction& CCPairFunction::multiply_with_op_inplace(const std::shared_ptr op) { + if (has_operator()) { + auto newop=combine(get_operator_ptr(),op); + reset_operator(newop); + } else { + reset_operator(op); + } + return *this; +} + double CCPairFunction::make_xy_u(const CCFunction& xx, const CCFunction& yy) const { double result = 0.0; @@ -38,8 +132,11 @@ CCPairFunction::make_xy_u(const CCFunction& xx, const CCFunction& yy) const { real_function_3d CCPairFunction::project_out(const CCFunction& f, const size_t particle) const { MADNESS_ASSERT(particle == 1 or particle == 2); real_function_3d result; - if (is_pure()) { - result = pure().get_function().project_out(f.function, particle - 1); // this needs 0 or 1 for particle but we give 1 or 2 + if (is_pure_no_op()) { + result = pure().get_function().project_out(f.function, + particle - 1); // this needs 0 or 1 for particle but we give 1 or 2 + } else if (is_op_pure()) { + MADNESS_EXCEPTION("implement CCPairFunction::project_out for op_pure",1); } else if (is_decomposed_no_op()) { result = project_out_decomposed(f.function, particle); } else if (is_op_decomposed()) { @@ -90,6 +187,127 @@ real_function_3d inner(const CCPairFunction& c, const real_function_3d& f, return c.partial_inner(f,v11,v22); } +CCPairFunction inner(const CCPairFunction& c1, const CCPairFunction& c2, + const std::tuple v1, const std::tuple v2) { + auto v11=std::array({std::get<0>(v1),std::get<1>(v1),std::get<2>(v1)}); + auto v22=std::array({std::get<0>(v2),std::get<1>(v2),std::get<2>(v2)}); + + return c1.partial_inner(c2,v11,v22); +} + +CCPairFunction CCPairFunction::partial_inner(const CCPairFunction& other, + const std::array& v1, + const std::array& v2) const { + auto a012=std::array{0,1,2}; + auto a345=std::array{3,4,5}; + MADNESS_CHECK(v1==a012 or v1== a345); + MADNESS_CHECK(v2==a012 or v2== a345); + MADNESS_CHECK(not this->is_op_pure()); // not implemented yet + MADNESS_CHECK(not other.is_op_pure()); // not implemented yet + + auto integration_index=[&a012,&a345](auto v) {return (v==a012) ? 0l : 1l;}; + auto remaining_index=[&integration_index](auto v) {return (integration_index(v)+1)%2;}; + + CCPairFunction result; + if (this->is_pure()) { + if (other.is_pure()) { + real_function_6d tmp=madness::innerXX<6>(this->get_function(),other.get_function(),v1,v2); + return CCPairFunction(tmp); + + } else if (other.is_decomposed_no_op()) { + // \int \sum_i f(1,2) a_i(1) b_i(3) d1 = \sum_i b_i(3) \int a_i(1) f(1,2) d1 + vector_real_function_3d tmp; + for (auto& a : other.get_vector(integration_index(v2))) { + tmp.push_back(innerXX<3>(this->get_function(),a,v1,a012)); // a012 is correct, referring to 3D function + } + return CCPairFunction(tmp,other.get_vector(remaining_index(v2))); + + } else if (other.is_op_decomposed()) { + + // \int \sum_i h(1,2) f(1,3) c_i(1) d_i(3) d1 + // = \sum_i d_i(3) \int h_c_i(1,2) f(1,3) d1 + // = \sum_i d_i(3) H_i(3,2) + const auto& h=this->pure().get_function(); + const auto& c=other.get_vector(integration_index(v2)); + const auto& d=other.get_vector(remaining_index(v2)); + auto& op=*(other.get_operator().get_op()); + op.particle()=integration_index(v1)+1; + + const vector_real_function_6d tmp=partial_mul(h,c,integration_index(v1)+1); + auto H=apply(world(),op,tmp); + real_function_6d result=real_factory_6d(world()); +// const vector_real_function_6d result=partial_mul(H,d,integration_index(v1)+1); + for (int i=0; iis_decomposed_no_op()) { + if (other.is_pure()) { + return other.partial_inner(*this,v2,v1); + } else if (other.is_decomposed_no_op()) { + // \int \sum_i a_i(1) b_i(2) \sum_j c_j(1) d_j(3) d1 + // = \sum_ij b_i(2) d_j(3) + // = \sum_i b~_i(2) d~_i(3) // SVD decomposition of S_ac + Tensor ovlp=matrix_inner(world(),this->get_vector(integration_index(v1)),other.get_vector(integration_index(v2))); + Tensor< typename Tensor::scalar_type > s; + Tensor U,VT; + svd(ovlp,U,s,VT); + for (int i=0; iget_vector(remaining_index(v1)),U); + auto right=transform(world(),other.get_vector(remaining_index(v2)),transpose(VT)); + return CCPairFunction(left,right); + + } else if (other.is_op_decomposed()) { + // \int \sum_ij a_i(1) b_i(2) f(1,3) c_j(1) d_j(3) d1 + // = \sum_ij b_i(2) d_j(3) \int ac_ij(1) f(1,3) d1 + // = \sum_i b_i(2) \sum_j d_j(3) g_ij(3) + // = \sum_i b_i(2) h_i(3) + const auto& a=this->get_vector(integration_index(v1)); + const auto& b=this->get_vector(remaining_index(v1)); + const auto& c=other.get_vector(integration_index(v2)); + const auto& d=other.get_vector(remaining_index(v2)); + const auto& op=*(other.get_operator().get_op()); + std::decay_t h(a.size()); // /same type as a, without reference& + for (int i=0; iis_op_decomposed()) { + if (other.is_pure()) { + return other.partial_inner(*this,v2,v1); + } else if (other.is_decomposed_no_op()) { + return other.partial_inner(*this,v2,v1); + } else if (other.is_op_decomposed()) { + if (this->is_convertible_to_pure_no_op()) { + CCPairFunction tmp=copy(*this); + tmp.convert_to_pure_no_op_inplace(); + return tmp.partial_inner(other,v1,v2); + } else if (other.is_convertible_to_pure_no_op()) { + CCPairFunction tmp=copy(other); + tmp.convert_to_pure_no_op_inplace(); + return this->partial_inner(tmp,v1,v2); + } else { + MADNESS_EXCEPTION("no partial_inner for this combination: ",1); + } + } else { + MADNESS_EXCEPTION("confused CCPairfunction",1); + } + } else { + MADNESS_EXCEPTION("confused CCPairfunction",1); + } + return result; + +} real_function_3d CCPairFunction::partial_inner(const real_function_3d& f, const std::array& v1, @@ -98,6 +316,7 @@ real_function_3d CCPairFunction::partial_inner(const real_function_3d& f, auto a345=std::array{3,4,5}; MADNESS_CHECK(v2==a012 ); // only 3 dimension in f MADNESS_CHECK(v1==a012 or v1== a345); // 6 dimension in f + MADNESS_CHECK(not this->is_op_pure()); // not implemented yet int particle=-1; if (v1== a012) particle=0; if (v1== a345) particle=1; @@ -165,91 +384,118 @@ CCPairFunction::assign_particles(const size_t particle) const { } } +/// compute the inner product of this and other +/// there are 4 possible components: pure/decomposed with and without operator, gives us 16 pair combinations.. double CCPairFunction::inner_internal(const CCPairFunction& other, const real_function_3d& R2) const { const CCPairFunction& f1=*this; const CCPairFunction& f2=other; - double thresh=FunctionDefaults<6>::get_thresh(); + double thresh=FunctionDefaults<6>::get_thresh()*0.1; double result = 0.0; - if (f1.is_pure() and f2.is_pure()) { - CCTimer tmp(world(), "making R1R2|u>"); - if (R2.is_initialized()) { - real_function_6d R1u = multiply(::copy(f1.pure().get_function()), ::copy(R2), 1); - real_function_6d R1R2u = multiply(R1u, ::copy(R2), 2); // R1u function now broken -// tmp.info(); - result = f2.pure().get_function().inner(R1R2u); + if (f1.is_pure() and f2.is_pure()) { // these are 4 combinations pure/pure + pureT bra=f1.get_function(); + pureT ket=f2.get_function(); +// if (R2.is_initialized()) { +// real_function_6d R1u = multiply(::copy(f1.pure().get_function()), ::copy(R2), 1); +// real_function_6d R1R2u = multiply(R1u, ::copy(R2), 2); // R1u function now broken +// bra = R1R2u; +// } + // include the operator(s), if any + if (f1.has_operator() or f2.has_operator()) { + MADNESS_EXCEPTION("still to debug",1); +// auto ops=combine(f1.get_operator_ptr(),f2.get_operator_ptr()); +// for (const auto& single_op : ops) { +// auto fac=single_op.first; +// auto op=single_op.second; +// double bla=0.0; +// if (op.get_op()) { +// real_function_6d tmp1; +// if (R2.is_initialized()) { +// tmp1= CompositeFactory(world()).g12(op.get_kernel()).ket(ket).particle1(R2).particle2(R2); +// } else { +// tmp1= CompositeFactory(world()).g12(op.get_kernel()).ket(ket); +// } +// bla=fac*inner(bra,tmp1); +// } else { +// bla=fac*inner(bra,ket); +// } +// result+=bla; +// } } else { - result = f2.pure().get_function().inner(f1.pure().get_function()); + // no operators + result=inner(bra,ket); } } else if (f1.is_pure() and f2.is_decomposed()) { // with or without operator - vector_real_function_3d a = R2.is_initialized() ? R2 * f2.get_a() : copy(world(), f2.get_a()); - vector_real_function_3d b = R2.is_initialized() ? R2 * f2.get_b() : copy(world(), f2.get_b()); - real_function_6d op; - if (f2.has_operator()) { - if (f2.decomposed().get_operator_ptr()->type() == OT_F12) { - op = TwoElectronFactory(world()).dcut(f2.get_operator().parameters.lo).gamma( - f2.get_operator().parameters.gamma).f12().thresh(thresh); - } else if (f2.get_operator().type() == OT_G12) { - op = TwoElectronFactory(world()).dcut(f2.get_operator().parameters.lo).thresh(thresh); - } else { - MADNESS_EXCEPTION(("6D Overlap with operatortype " + assign_name(f2.get_operator().type()) + " not supported").c_str(), 1); + const vector_real_function_3d a = R2.is_initialized() ? R2 * f2.get_a() : copy(world(), f2.get_a()); + const vector_real_function_3d b = R2.is_initialized() ? R2 * f2.get_b() : copy(world(), f2.get_b()); + const pureT& bra=f1.get_function(); + + auto ops=combine(f1.get_operator_ptr(),f2.get_operator_ptr()); + MADNESS_EXCEPTION("still to debug",1); +// if (ops.size()>0) { +// for (const auto& single_op : ops) { +// auto fac = single_op.first; +// auto op = single_op.second; +// double bla=0.0; +// for (int i=0; i(world()).g12(op.get_kernel()).particle1(a[i]).particle2(b[i]); +// bla += fac * inner(bra, tmp); +// } else { +// real_function_6d tmp = CompositeFactory(world()).particle1(a[i]).particle2(b[i]); +// bla += fac * inner(bra,tmp); +// } +// } +// result+=bla; +// } +// } else { // no operators + for (int i=0; i(world()).particle1(a[i]).particle2(b[i]); + result+=inner(bra,tmp); } - } - for (int i=0; i(world()).g12(op).particle1(x).particle2(y); - else opxy= CompositeFactory(world()).particle1(x).particle2(y); - result+= f1.pure().get_function().inner(opxy); - } +// } } else if (f1.is_decomposed() and f2.is_pure()) { // with or without op result= f2.inner_internal(f1,R2); - } else if (f1.is_decomposed_no_op() and f2.is_decomposed_no_op()) { + + } else if (f1.is_decomposed() and f2.is_decomposed()) { MADNESS_ASSERT(f1.get_a().size() == f1.get_b().size()); MADNESS_ASSERT(f2.get_a().size() == f2.get_b().size()); - vector_real_function_3d a = R2.is_initialized() ? R2* f2.get_a() : f2.get_a(); - vector_real_function_3d b = R2.is_initialized() ? R2* f2.get_b() : f2.get_b(); - // = \sum_ij = \sum_ij - result = (matrix_inner(world(), a, f1.get_a()).emul(matrix_inner(world(), b, f1.get_b()))).sum(); - } else if (f1.is_decomposed_no_op() and f2.is_op_decomposed()) { - MADNESS_ASSERT(f1.get_a().size() == f1.get_b().size()); - // = const vector_real_function_3d& a1 = f1.get_a(); const vector_real_function_3d& b1 = f1.get_b(); - const vector_real_function_3d a2 = R2.is_initialized() ? R2 * f2.get_a() : copy(world(),f2.get_a()); - const vector_real_function_3d b2 = R2.is_initialized() ? R2 * f2.get_b() : copy(world(),f2.get_b()); - for (size_t i = 0; i < a1.size(); i++) { - vector_real_function_3d aa = truncate(a1[i] * a2); - vector_real_function_3d aopx = f2.get_operator()(aa); - vector_real_function_3d bb = truncate(b1[i] * b2); - result += inner(bb,aopx); - } - - } else if (f1.is_op_decomposed() and f2.is_decomposed_no_op()) { - return f2.inner_internal(f1,R2); - - } else if (f1.is_op_decomposed() and f2.is_op_decomposed()) { - MADNESS_CHECK(can_combine(f1.get_operator(),f2.get_operator())); - MADNESS_CHECK(f1.get_a().size()==1); - std::vector> ops=combine(f1.get_operator(),f2.get_operator()); - CCPairFunction bra(f1.get_a(),f1.get_b()); - for (const auto& op : ops) { - CCPairFunction ket; - if (op.second.get_op()) ket = CCPairFunction(&op.second,f2.get_a(),f2.get_b()); - else ket = CCPairFunction(f2.get_a(),f2.get_b()); - - double tmp=op.first * inner(ket,bra,R2); - print("inner",bra.name(true),ket.name()," : ",tmp); - result+=tmp; - } + const vector_real_function_3d a2 = R2.is_initialized() ? R2* f2.get_a() : f2.get_a(); + const vector_real_function_3d b2 = R2.is_initialized() ? R2* f2.get_b() : f2.get_b(); + + + MADNESS_EXCEPTION("still to debug",1); + auto ops=combine(f1.get_operator_ptr(),f2.get_operator_ptr()); +// if (ops.size()==0) { +// // = \sum_ij = \sum_ij +// result = (matrix_inner(world(), a1, a2)).trace(matrix_inner(world(),b1,b2)); +// } else { +// // = +// for (const auto& single_op : ops) { +// auto fac = single_op.first; +// auto op = single_op.second; +// +// double bla=0.0; +// if (op.get_op()) { +// for (size_t i = 0; i < a1.size(); i++) { +// vector_real_function_3d aa = truncate(a1[i] * a2); +// vector_real_function_3d bb = truncate(b1[i] * b2); +// vector_real_function_3d aopx = op(aa); +// bla += fac * inner(bb, aopx); +// } +// } else { +// bla += fac*(matrix_inner(world(), a1, a2)).trace(matrix_inner(world(),b1,b2)); +// } +// result+=bla; +// } +// } } else MADNESS_EXCEPTION( - ("CCPairFunction Overlap not supported for combination " + f1.name() + " and " + f2.name()).c_str(), 1) - - ; + ("CCPairFunction Overlap not supported for combination " + f1.name() + " and " + f2.name()).c_str(), 1) ; return result; } @@ -267,6 +513,7 @@ std::vector apply(const ProjectorBase& projector, const std::vec std::vector result; for (const auto& pf : argument) { if (pf.is_pure()) { + MADNESS_CHECK(not pf.has_operator()); // not yet implemented if (auto SO=dynamic_cast*>(&projector)) { auto tmp=(*SO)(pf.get_function()); auto tmp2=CCPairFunction(tmp); @@ -303,6 +550,7 @@ std::vector apply(const ProjectorBase& projector, const std::vec } } else if (pf.is_op_decomposed()) { if (auto SO=dynamic_cast*>(&projector)) { +// CCTimer t(world,"SO block"); // Q12 = 1 - O1 (1 - 1/2 O2) - O2 (1 - 1/2 O1) QProjector Q1(world,SO->bra1(),SO->ket1()); Q1.set_particle(0); @@ -310,8 +558,11 @@ std::vector apply(const ProjectorBase& projector, const std::vec Q2.set_particle(1); auto tmp=Q1(Q2(std::vector({pf}))); for (auto& t: tmp) result.push_back(t); +// t.print(); } else if (auto P=dynamic_cast*>(&projector)) { +// CCTimer t(world,"P block"); + // Q12 = 1 - O1 (1 - 1/2 O2) - O2 (1 - 1/2 O1) std::vector tmp= zero_functions_compressed(world,P->get_ket_vector().size()); // per term a_i b_k: @@ -330,8 +581,10 @@ std::vector apply(const ProjectorBase& projector, const std::vec truncate(world,tmp); if (P->get_particle()==0) result.push_back(CCPairFunction(P->get_ket_vector(),tmp)); if (P->get_particle()==1) result.push_back(CCPairFunction(tmp,P->get_ket_vector())); +// t.print(); } else if (auto Q=dynamic_cast*>(&projector)) { +// CCTimer t(world,"Q block"); // Q1 f12 |a_i b_i> = f12 |a_i b_i> - \sum_k |k(1) a_i(2)*f_(kb_i)(2) > result.push_back(pf); // reuse the projector code above @@ -340,6 +593,7 @@ std::vector apply(const ProjectorBase& projector, const std::vec t*=-1.0; result.push_back(t); } +// t.print(); } else { MADNESS_EXCEPTION("CCPairFunction: unknown projector type",1); @@ -354,5 +608,36 @@ std::vector apply(const ProjectorBase& projector, const std::vec return result; }; +template +std::vector apply(const SeparatedConvolution& op, const std::vector& argument) { + if (argument.size()==0) return argument; + World& world=argument.front().world(); + std::vector result; + for (const auto& arg : argument) { + if (arg.is_pure()) { + result.push_back(CCPairFunction(op(arg.get_function()))); + } else if (arg.is_op_pure()) { + auto tmp=arg.to_pure(); + result.push_back(apply(op,tmp)); + } else if (arg.is_decomposed_no_op()) { + MADNESS_CHECK(op.particle()==1 or op.particle()==2); + if (op.particle()==1) { + auto tmp= madness::apply(world,op,arg.get_a()); + result.push_back(CCPairFunction(tmp,arg.get_b())); + } else if (op.particle()==2) { + auto tmp= madness::apply(world,op,arg.get_b()); + result.push_back(CCPairFunction(arg.get_a(),tmp)); + } + } else if (arg.is_op_decomposed()) { // sucks.. + auto tmp=arg.to_pure(); + result.push_back(apply(op,tmp)); + } + } + + return result; +} + + +template std::vector apply(const SeparatedConvolution& op, const std::vector& argument); } // namespace madness \ No newline at end of file diff --git a/src/madness/chem/ccpairfunction.h b/src/madness/chem/ccpairfunction.h index 4ddc4af6e93..32d13fc2e60 100644 --- a/src/madness/chem/ccpairfunction.h +++ b/src/madness/chem/ccpairfunction.h @@ -23,11 +23,6 @@ class ProjectorBase; /// Types of Functions used by CC_function class enum FuncType { UNDEFINED, HOLE, PARTICLE, MIXED, RESPONSE }; -/// FuncTypes used by the CC_function_6d structure -enum PairFormat { PT_UNDEFINED, PT_FULL, PT_DECOMPOSED, PT_OP_DECOMPOSED }; - - - /// structure for a CC Function 3D which holds an index and a type // the type is defined by the enum FuncType (definition at the start of this file) struct CCFunction : public archive::ParallelSerializableObject { @@ -95,7 +90,9 @@ class TwoBodyFunctionComponentBase { virtual void swap_particles_inplace() = 0; virtual bool is_pure() const {return false;} virtual bool is_decomposed() const {return false;} - virtual bool has_operator() const {return false;} + virtual bool has_operator() const = 0; + virtual void set_operator(const std::shared_ptr op) = 0; + virtual const std::shared_ptr get_operator_ptr() const = 0; virtual void print_size() const = 0; virtual std::string name(const bool transpose=false) const = 0; virtual World& world() const =0; @@ -110,11 +107,12 @@ class TwoBodyFunctionPureComponent : public TwoBodyFunctionComponentBase { TwoBodyFunctionPureComponent() = default; explicit TwoBodyFunctionPureComponent(const Function& f) : u(f) {} + explicit TwoBodyFunctionPureComponent(const std::shared_ptr op, const Function& f) + : u(f), op(op) {} /// deep copy std::shared_ptr clone() override { - TwoBodyFunctionPureComponent result; - result.u=madness::copy(u); + TwoBodyFunctionPureComponent result(op,madness::copy(u)); return std::make_shared>(result); } @@ -126,6 +124,8 @@ class TwoBodyFunctionPureComponent : public TwoBodyFunctionComponentBase { bool is_pure() const override {return true;} + bool has_operator() const override {return op!=nullptr;} + World& world() const override {return u.world();}; void serialize() {} @@ -135,8 +135,12 @@ class TwoBodyFunctionPureComponent : public TwoBodyFunctionComponentBase { } std::string name(const bool transpose) const override { - if (transpose) return "< u |"; - return "|u>"; + if (transpose) { + if (has_operator()) return "< u |"+get_operator_ptr()->name(); + return "< u |"; + } + if (has_operator()) return get_operator_ptr()->name() + "| u >"; + return "| u >"; } @@ -151,6 +155,10 @@ class TwoBodyFunctionPureComponent : public TwoBodyFunctionComponentBase { u=swap_particles(u); } + const std::shared_ptr get_operator_ptr() const override {return op;}; + + void set_operator(const std::shared_ptr op1) override {op=op1;} + real_function_6d& get_function() { return u; } @@ -158,6 +166,7 @@ class TwoBodyFunctionPureComponent : public TwoBodyFunctionComponentBase { private: /// pure 6D function real_function_6d u; + std::shared_ptr op; }; @@ -174,7 +183,7 @@ class TwoBodyFunctionSeparatedComponent : public TwoBodyFunctionComponentBase { TwoBodyFunctionSeparatedComponent(const std::vector>& a, const std::vector>& b, - const CCConvolutionOperator* op) : a(a), b(b), op(op) {}; + const std::shared_ptr op) : a(a), b(b), op(op) {}; TwoBodyFunctionSeparatedComponent(const TwoBodyFunctionSeparatedComponent& other) = default; @@ -208,32 +217,50 @@ class TwoBodyFunctionSeparatedComponent : public TwoBodyFunctionComponentBase { std::string name(const bool transpose) const override { if (transpose) { - if (has_operator()) return "name(); - return "name(); + return "name() + "|xy>"; - return "|ab>"; + if (has_operator()) return get_operator_ptr()->name() + "| ab>"; + return "| ab>"; }; void serialize() {} template - TwoBodyFunctionPureComponent apply(const SeparatedConvolution* op, const int particle=0) {} + TwoBodyFunctionPureComponent apply(const SeparatedConvolution* op, const int particle=0) { + MADNESS_EXCEPTION("TwoBodyFunctionPureComponent apply not yet implemented",1); + } /// return f(2,1) void swap_particles_inplace() override { std::swap(a,b); } + long rank() const { + MADNESS_CHECK(a.size()==b.size()); + return a.size(); + } + std::vector> get_a() const {return a;} std::vector> get_b() const {return b;} - const CCConvolutionOperator* get_operator_ptr() const {return op;}; + std::vector> get_vector(const int i) const { + MADNESS_CHECK(i==0 or i==1); + if (i==0) return a; + else if (i==1) return b; + else { + MADNESS_EXCEPTION("confused index in TwoBodyFunctionSeparatedComponent",1); + } + } + + const std::shared_ptr get_operator_ptr() const override {return op;}; + + void set_operator(const std::shared_ptr op1) override {op=op1;} private: std::vector> a; std::vector> b; - const CCConvolutionOperator* op=nullptr; + std::shared_ptr op; }; @@ -263,9 +290,21 @@ class TwoBodyFunctionSeparatedComponent : public TwoBodyFunctionComponentBase { * - mul_partial */ +/// a 6D function, either in full or low rank form, possibly including an 2-particle function + +/** + * the function is stored as + * - pure: full rank form, 6D + * - op_pure: full rank form, 6D with an 2-particle function f(1,2) |u> + * - decomposed: sum of two vectors of 3D functions \sum_i |a_i(1) b_i(2)> + * - op_decomposed: as above, with an 2-particle function: f(1,2) \sum_i |a_i b_i> + * +**/ struct CCPairFunction { using T=double; +using pureT=Function; + public: /// empty ctor @@ -276,29 +315,45 @@ using T=double; component.reset(new TwoBodyFunctionPureComponent(copy(ket))); } + /// takes a deep copy of the argument function + explicit CCPairFunction(const std::shared_ptr op_, const real_function_6d& ket) { + component.reset(new TwoBodyFunctionPureComponent(op_,copy(ket))); + } + /// takes a deep copy of the argument functions explicit CCPairFunction(const vector_real_function_3d& f1, const vector_real_function_3d& f2) { World& world=f1.front().world(); component.reset(new TwoBodyFunctionSeparatedComponent(copy(world,f1),copy(world,f2))); } + /// takes a deep copy of the argument functions + explicit CCPairFunction(const real_function_3d& f1, const real_function_3d& f2) : + CCPairFunction(std::vector({f1}),std::vector({f2})) { + } + /// takes a deep copy of the argument functions explicit CCPairFunction(const std::pair& f) : CCPairFunction(f.first,f.second) { } /// takes a deep copy of the argument functions - explicit CCPairFunction(const CCConvolutionOperator *op_, const CCFunction& f1, const CCFunction& f2) : + explicit CCPairFunction(const std::shared_ptr op_, const CCFunction& f1, const CCFunction& f2) : CCPairFunction(op_,std::vector({f1.function}),std::vector({f2.function})) { } /// takes a deep copy of the argument functions - explicit CCPairFunction(const CCConvolutionOperator *op_, const std::vector& f1, + explicit CCPairFunction(const std::shared_ptr op_, const std::vector& f1, const std::vector& f2) { World& world=f1.front().world(); component.reset(new TwoBodyFunctionSeparatedComponent(copy(world,f1),copy(world,f2),op_)); } + /// takes a deep copy of the argument functions + explicit CCPairFunction(const std::shared_ptr op_, const real_function_3d& f1, + const real_function_3d& f2) : CCPairFunction(op_,std::vector({f1}), + std::vector({f2})) { + }; + /// shallow assignment operator CCPairFunction& operator()(const CCPairFunction& other) { component=other.component; @@ -315,6 +370,10 @@ using T=double; return result; } + /// add all like functions up, return *this for chaining + friend std::vector consolidate(const std::vector& other); + + void info() const { print_size(); } World& world() const { @@ -336,21 +395,58 @@ using T=double; /// deep copy necessary otherwise: shallow copy errors CCPairFunction invert_sign(); + /// scalar multiplication: f*fac CCPairFunction operator*(const double fac) const { CCPairFunction result=copy(*this); result*=fac; return result; } + /// scalar multiplication: fac*f friend CCPairFunction operator*(const double fac, const CCPairFunction& f) { return fac*f; } + /// multiplication with a 2-particle function + friend CCPairFunction operator*(const std::shared_ptr op, const CCPairFunction& f) { + CCPairFunction result=copy(f); + return result.multiply_with_op_inplace(op); + } + + /// multiplication with a 2-particle function + friend std::vector operator*(const std::shared_ptr op, + const std::vector& f) { + std::vector result; + for (auto& ff : f) { + result.push_back(copy(ff)); + result.back().multiply_with_op_inplace(op); + } + return result; + } + + friend std::vector multiply(const std::vector& other, const real_function_3d f, + const std::array& v1) { + std::vector result; + for (auto& o : other) result.push_back(multiply(o,f,v1)); + return result; + } + + /// multiplication with a 2-particle function + CCPairFunction operator*(const std::shared_ptr op) { + CCPairFunction result=copy(*this); + return result.multiply_with_op_inplace(op); + } + + CCPairFunction& multiply_with_op_inplace(const std::shared_ptr op); + + + bool has_operator() const {return component->has_operator();} bool is_pure() const {return component->is_pure();} + bool is_op_pure() const {return is_pure() and has_operator();} + bool is_pure_no_op() const {return is_pure() and (not has_operator());} bool is_decomposed() const {return component->is_decomposed();} - bool is_decomposed_no_op() const {return component->is_decomposed() and (not component->has_operator());} bool is_op_decomposed() const {return component->is_decomposed() and component->has_operator();} - bool has_operator() const {return component->has_operator();} + bool is_decomposed_no_op() const {return component->is_decomposed() and (not component->has_operator());} TwoBodyFunctionPureComponent& pure() const { if (auto ptr=dynamic_cast*>(component.get())) return *ptr; @@ -366,16 +462,46 @@ using T=double; MADNESS_CHECK(component->is_decomposed()); return decomposed().get_a(); } + vector_real_function_3d get_b() const { MADNESS_CHECK(component->is_decomposed()); return decomposed().get_b(); } + std::vector> get_vector(const int i) const { + MADNESS_CHECK(component->is_decomposed()); + return decomposed().get_vector(i); + } + const CCConvolutionOperator& get_operator() const { - MADNESS_CHECK(is_op_decomposed()); - return *decomposed().get_operator_ptr(); + MADNESS_CHECK(component and component->has_operator()); + return *(component->get_operator_ptr()); + } + + const std::shared_ptr get_operator_ptr() const { + MADNESS_CHECK(component); + return component->get_operator_ptr(); } + void reset_operator(const std::shared_ptr op) { + MADNESS_CHECK(component); + component->set_operator(op); + } + + /// can this be converted to a pure representation (depends on the operator, if present) + bool is_convertible_to_pure_no_op() const; + + /// out-of-place conversion to pure function + CCPairFunction to_pure() const { + auto tmp=copy(*this); + MADNESS_CHECK(tmp.is_convertible_to_pure_no_op()); + tmp.convert_to_pure_no_op_inplace(); + return tmp; + } + + /// convert this into a pure hi-dim function + void convert_to_pure_no_op_inplace(); + CCPairFunction& operator*=(const double fac) { if (component->is_pure()) pure()*=fac; if (component->is_decomposed()) decomposed()*=fac; @@ -392,12 +518,16 @@ using T=double; return component->name(transpose); } + /// multiply CCPairFunction with a 3D function of one of the two particles + friend CCPairFunction multiply(const CCPairFunction& other, const real_function_3d& f, const std::array& v1); + /// @param[in] f: a 3D-CC_function /// @param[in] particle: the particle on which the operation acts /// @param[out] _particle (projection from 6D to 3D) real_function_3d project_out(const CCFunction& f, const size_t particle) const; - // result is: _particle + /// result is: _particle + /// @param[in] x: a 3D-CC_function /// @param[in] op: a CC_convoltion_operator which is currently either f12 or g12 /// @param[in] particle: the particle on which the operation acts (can be 1 or 2) @@ -415,6 +545,7 @@ using T=double; double make_xy_u(const CCFunction& xx, const CCFunction& yy) const; + /// compute the inner product of this and other double inner_internal(const CCPairFunction& other, const real_function_3d& R2) const; friend double inner(const CCPairFunction& a, const CCPairFunction& b, const real_function_3d& R2) { @@ -426,13 +557,17 @@ using T=double; return a.inner_internal(b,R2); } - friend double inner(const std::vector& va, const std::vector& vb) { - real_function_3d R2; + friend double inner(const std::vector& va, const std::vector& vb, + const real_function_3d R2=real_function_3d()) { + double wall0=cpu_time(); +// real_function_3d R2; double result=0.0; for (auto& a : va) { for (auto& b : vb) { double tmp=a.inner_internal(b,R2); - print("result from inner",a.name(true),b.name(),tmp); + double wall1=cpu_time(); + print("result from inner",a.name(true),b.name(),tmp,wall1-wall0,"s"); + wall0=wall1; result+=tmp; } } @@ -444,22 +579,6 @@ using T=double; /// the 3 types of 6D-function that occur in the CC potential which coupled doubles to singles std::shared_ptr component; -// World& world; -// /// the type of the given 6D-function -// const PairFormat type; -// /// if type==decomposed this is the first particle -// vector_real_function_3d a; -// /// if type==decomposed this is the second particle -// vector_real_function_3d b; -// /// if type==op_decomposed_ this is the symmetric 6D-operator (g12 or f12) in u=op12|xy> -// const CCConvolutionOperator *op; -// /// if type==op_decomposed_ this is the first particle in u=op12|xy> -// CCFunction x; -// /// if type==op_decomposed_ this is the second particle in u=op12|xy> -// CCFunction y; -// /// if type=pure_ this is just the MRA 6D-function -// real_function_6d u; - /// @param[in] f: a 3D-CC_function /// @param[in] particle: the particle on which the operation acts /// @param[out] _particle (projection from 6D to 3D) for the case that u=|ab> so _particle = *|b> if particle==1 @@ -505,6 +624,10 @@ using T=double; const std::array& v1, const std::array& v2) const; + CCPairFunction partial_inner(const CCPairFunction& other, + const std::array& v1, + const std::array& v2) const; + }; /// apply the projector on the argument function, potentially yielding a vector of CCPairfunctions as result @@ -518,9 +641,19 @@ std::vector apply(const ProjectorBase& P, const std::vector +std::vector apply(const SeparatedConvolution& op, const std::vector& argument); + +/// convenience function +template +CCPairFunction apply(const SeparatedConvolution& op, const CCPairFunction& argument); + real_function_3d inner(const CCPairFunction& c, const real_function_3d& f, const std::tuple v1, const std::tuple v2={0,1,2}); +CCPairFunction inner(const CCPairFunction& c1, const CCPairFunction& c2, + const std::tuple v1, const std::tuple v2); } // namespace madness diff --git a/src/madness/chem/lowrankfunction.h b/src/madness/chem/lowrankfunction.h new file mode 100644 index 00000000000..5707f0d7688 --- /dev/null +++ b/src/madness/chem/lowrankfunction.h @@ -0,0 +1,949 @@ +// +// Created by Florian Bischoff on 8/10/23. +// + +#ifndef MADNESS_LOWRANKFUNCTION_H +#define MADNESS_LOWRANKFUNCTION_H + + +#include +#include +#include +#include +#include + + +namespace madness { + + struct LowRankFunctionParameters : QCCalculationParametersBase { + + LowRankFunctionParameters() : QCCalculationParametersBase() { + + // initialize with: key, value, comment (optional), allowed values (optional) + initialize("radius",2.0,"the radius"); + initialize("gamma",1.0,"the exponent of the correlation factor"); + initialize("volume_element",0.1,"volume covered by each grid point"); + initialize("hard_shell",true,"radius is hard"); + initialize("tol",1.e-8,"rank-reduced cholesky tolerance"); + initialize("f12type","Slater","correlation factor",{"Slater","SlaterF12"}); + initialize("orthomethod","cholesky","orthonormalization",{"cholesky","canonical","symmetric"}); + initialize("transpose","slater2","transpose of the matrix",{"slater1","slater2"}); + initialize("gridtype","random","the grid type",{"random","cartesian","spherical"}); + initialize("rhsfunctiontype","exponential","the type of function",{"delta","exponential"}); + initialize("optimize",1,"number of optimization iterations"); + } + + void read_and_set_derived_values(World& world, const commandlineparser& parser, std::string tag) { + read_input_and_commandline_options(world,parser,tag); + } + + double radius() const {return get("radius");} + double gamma() const {return get("gamma");} + double volume_element() const {return get("volume_element");} + double tol() const {return get("tol");} + bool hard_shell() const {return get("hard_shell");} + int optimize() const {return get("optimize");} + std::string gridtype() const {return get("gridtype");} + std::string orthomethod() const {return get("orthomethod");} + std::string rhsfunctiontype() const {return get("rhsfunctiontype");} + std::string f12type() const {return get("f12type");} + }; + + + class gridbase { + protected: + double volume_element=0.1; + double radius=3; + }; + + template + class randomgrid : public gridbase { + public: + randomgrid(const double volume_element, const double radius) : gridbase() { + this->volume_element=volume_element; + this->radius=radius; + } + + std::vector> get_grid() const { + std::vector> grid; + long npoint_within_volume=volume()/volume_element; + print("npoint_within_volume",npoint_within_volume); + + auto cell = FunctionDefaults::get_cell(); + auto is_in_cell = [&cell](const Vector& r) { + for (int d = 0; d < NDIM; ++d) if (r[d] < cell(d, 0) or r[d] > cell(d, 1)) return false; + return true; + }; + double rad=radius; + auto is_in_sphere = [&rad](const Vector& r) { + return (r.normf()0 and NDIM<4); + if (NDIM==1) return 2.0*radius; + if (NDIM==2) return constants::pi*radius*radius; + if (NDIM==3) return 4.0 / 3.0 * constants::pi * std::pow(radius, 3.0); + } + + static Vector gaussian_random_distribution(double mean, double variance) { + std::random_device rd{}; + std::mt19937 gen{rd()}; + std::normal_distribution<> d{mean, variance}; + Vector result; + for (int i = 0; i < NDIM; ++i) result[i]=d(gen); + return result; + } + +// double computed_volume_element() const { +// double volume = 4.0 / 3.0 * constants::pi * std::pow(radius, 3.0); +// print("volume element in random grid", volume / (0.67 * rank)); +// } + + }; + + template + class sphericalgrid : public gridbase { + public: + double volume_element=0.1; + double radius=3; + bool hard_shell=true; + + sphericalgrid(const double volume_element, const double radius, bool hard_shell) + :volume_element(volume_element), radius(radius), hard_shell(hard_shell) { + }; + + std::vector> get_coordinates() const { + // 1D grid + double volume_element_1d=std::pow(volume_element,1./NDIM); + long ngrid=std::ceil(radius/volume_element_1d); + double stepsize=radius/ngrid; + double scale=1.0; + if (not hard_shell) scale=std::pow(2.0,1.0/(ngrid+1)); + print("scale",scale); + + std::vector coord1d; + print("volume element, stepsize, ngrid" ,volume_element, std::pow(stepsize,NDIM),stepsize,ngrid); + for (int i=0; i> result; + for (int i=0; i c{coord1d[i],coord1d[j],coord1d[k]}; + double cutoff = hard_shell ? radius : 2.0*radius; + if (c.normf() + struct cartesian_grid { + Vector lovec,hivec; + std::vector stride; + long index=0; + long n_per_dim; + long total_n; + Vector increment; + + cartesian_grid(const double volume_per_gridpoint, const double radius) { + double length_per_gridpoint=std::pow(volume_per_gridpoint,1.0/NDIM); + n_per_dim=ceil(2.0*radius/length_per_gridpoint); + print("length per gridpoint, n_per_dim",length_per_gridpoint,n_per_dim); + print("volume_per_gridpoint",std::pow(length_per_gridpoint,NDIM)); + initialize(-radius,radius); + print("increment",increment); + } + + + cartesian_grid(const long n_per_dim, const double lo, const double hi) + : n_per_dim(n_per_dim) { + initialize(lo,hi); + } + + cartesian_grid(const cartesian_grid& other) : lovec(other.lovec), + hivec(other.hivec), stride(other.stride), index(0), n_per_dim(other.n_per_dim), + total_n(other.total_n), increment(other.increment) { + } + + cartesian_grid& operator=(const cartesian_grid& other) { + cartesian_grid tmp(other); + std::swap(*this,other); + return *this; + } + + void initialize(const double lo, const double hi) { + lovec.fill(lo); + hivec.fill(hi); + increment=(hivec-lovec)*(1.0/double(n_per_dim-1)); + stride=std::vector(NDIM,1l); + total_n=std::pow(n_per_dim,NDIM); + for (long i=NDIM-2; i>=0; --i) stride[i]=n_per_dim*stride[i+1]; + + } + + double volume_per_gridpoint() const{ + double volume=1.0; + for (int i=0; i get_coordinates() const { + Vector tmp(NDIM); + for (int idim=0; idim +struct particle { + std::array dims; + + /// default constructor + particle() = default; + + /// convenience for particle 1 (the left/first particle) + static particle particle1() { + particle p; + for (int i=0; i(1,p)) {} + particle(const int p1, const int p2) : particle(std::vector({p1,p2})) {} + particle(const int p1, const int p2,const int p3) : particle(std::vector({p1,p2,p3})) {} + particle(const std::vector p) { + for (int i=0; i + typename std::enable_if_t> + get_tuple() const {return std::tuple(dims[0]);} + + template + typename std::enable_if_t> + get_tuple() const {return std::tuple(dims[0],dims[1]);} + + template + typename std::enable_if_t> + get_tuple() const {return std::tuple(dims[0],dims[1],dims[2]);} +}; + +template +std::ostream& operator<<(std::ostream& os, const particle& p) { + os << "("; + for (auto i=0; i +struct LRFunctorBase { + + virtual ~LRFunctorBase() {}; + virtual std::vector> inner(const std::vector>& rhs, + const particle p1, const particle p2) const =0; + + virtual Function inner(const Function& rhs, const particle p1, const particle p2) const { + return inner(std::vector>({rhs}),p1,p2)[0]; + } + + virtual T operator()(const Vector& r) const =0; + virtual typename Tensor::scalar_type norm2() const { + MADNESS_EXCEPTION("L2 norm not implemented",1); + } + + virtual World& world() const =0; + friend std::vector> inner(const LRFunctorBase& functor, const std::vector>& rhs, + const particle p1, const particle p2) { + return functor.inner(rhs,p1,p2); + } + friend Function inner(const LRFunctorBase& functor, const Function& rhs, + const particle p1, const particle p2) { + return functor.inner(rhs,p1,p2); + } + +}; + +template +struct LRFunctorF12 : public LRFunctorBase { + LRFunctorF12() = default; + LRFunctorF12(const std::shared_ptr> f12, const Function& a, + const Function& b) : f12(f12), a(a), b(b) {} + + std::shared_ptr> f12; ///< a two-particle function + Function a,b; ///< the lo-dim functions + + World& world() const {return f12->get_world();} + std::vector> inner(const std::vector>& rhs, + const particle p1, const particle p2) const { + + std::vector> result; + // functor is now a(1) b(2) f12 + // result(1) = \int a(1) f(1,2) b(2) rhs(2) d2 + World& world=rhs.front().world(); + auto premultiply= p1.is_first() ? a : b; + auto postmultiply= p1.is_first() ? b : a; + + const int nbatch=30; + for (int i=0; i> tmp; + auto begin= rhs.begin()+i; + auto end= (i+nbatch)::scalar_type norm2() const { + const Function one=FunctionFactory(world()).f([](const Vector& r){return 1.0;}); + const Function pre=(a.is_initialized()) ? a : one; + const Function post=(b.is_initialized()) ? b : one; + const SeparatedConvolution& f12a=*(f12); + const SeparatedConvolution f12sq= SeparatedConvolution::combine(f12a,f12a); + + // \int f(1,2)^2 d1d2 = \int f(1,2)^2 pre(1)^2 post(2)^2 d1 d2 + typename Tensor::scalar_type term1 =madness::inner(post*post,f12sq(pre*pre)); + return sqrt(term1); + + } + + T operator()(const Vector& r) const { + + auto split = [](const Vector& r) { + Vector first, second; + for (int i=0; iinfo.mu; + auto [first,second]=split(r); + + double result=1.0; + if (a.is_initialized()) result*=a(first); + if (b.is_initialized()) result*=b(first); + if (f12->info.type==OT_SLATER) result*=exp(-gamma*(first-second).normf()); + else if (f12->info.type==OT_GAUSS) result*=exp(-gamma* madness::inner(first-second,first-second)); + else return 1.0; + return result; + + } +}; + +template +struct LRFunctorPure : public LRFunctorBase { + LRFunctorPure() = default; + LRFunctorPure(const Function& f) : f(f) {} + World& world() const {return f.world();} + + Function f; ///< a hi-dim function + + std::vector> inner(const std::vector>& rhs, + const particle p1, const particle p2) const { + std::vector> result; + for (const auto& r : rhs) result.push_back(madness::inner(f,r,p1.get_tuple(),p2.get_tuple())); + return result; + } + + T operator()(const Vector& r) const { + return f(r); + } + + typename Tensor::scalar_type norm2() const { + return f.norm2(); + } +}; + + + /// LowRankFunction represents a hi-dimensional (NDIM) function as a sum of products of low-dimensional (LDIM) functions + + /// f(1,2) = \sum_i g_i(1) h_i(2) + /// a LowRankFunction can be created from a hi-dim function directly, or from a composite like f(1,2) phi(1) psi(2), + /// where f(1,2) is a two-particle function (e.g. a Slater function) + template + class LowRankFunction { + public: + + World& world; + double rank_revealing_tol=1.e-8; // rrcd tol + std::string orthomethod="canonical"; + bool do_print=true; + std::vector> g,h; + const particle p1=particle::particle1(); + const particle p2=particle::particle2(); + + LowRankFunction(World& world) : world(world) {} + + LowRankFunction(std::vector> g, std::vector> h, + double tol, std::string orthomethod) : world(g.front().world()), + rank_revealing_tol(tol), orthomethod(orthomethod), g(g), h(h) { + + } + + /// shallow copy ctor + LowRankFunction(const LowRankFunction& other) : world(other.world), + g(other.g), h(other.h), + rank_revealing_tol(other.rank_revealing_tol), orthomethod(other.orthomethod) { + } + + /// deep copy + friend LowRankFunction copy(const LowRankFunction& other) { + return LowRankFunction(madness::copy(other.g),madness::copy(other.h),other.rank_revealing_tol,other.orthomethod); + } + + LowRankFunction& operator=(const LowRankFunction& f) { // Assignment required for storage in vector + LowRankFunction ff(f); + std::swap(ff.g,g); + std::swap(ff.h,h); + return *this; + } + + /// function evaluation + T operator()(const Vector& r) const { + Vector first, second; + for (int i=0; i result=copy(*this); + result+=b; + return result; + } + /// subtraction + LowRankFunction operator-(const LowRankFunction& b) const { + LowRankFunction result=copy(*this); + result-=b; + return result; + } + + /// in-place addition + LowRankFunction& operator+=(const LowRankFunction& b) { + + g=append(g,copy(b.g)); + h=append(h,copy(b.h)); + return *this; + } + + /// in-place subtraction + LowRankFunction& operator-=(const LowRankFunction& b) { + g=append(g,-1.0*b.g); // operator* implies deep copy of b.g + h=append(h,copy(b.h)); + return *this; + } + + /// scale by a scalar + template + LowRankFunction operator*(const Q a) const { + return LowRankFunction,NDIM>(g * a, Q(h),rank_revealing_tol,orthomethod); + } + + /// out-of-place scale by a scalar (no type conversion) + LowRankFunction operator*(const T a) const { + return LowRankFunction(g * a, h,rank_revealing_tol,orthomethod); + } + + /// multiplication with a scalar + friend LowRankFunction operator*(const T a, const LowRankFunction& other) { + return other*a; + } + + /// in-place scale by a scalar (no type conversion) + LowRankFunction& operator*=(const T a) { + g=g*a; + return *this; + } + + /// l2 norm + typename TensorTypeData::scalar_type norm2() const { + auto tmp1=matrix_inner(world,h,h); + auto tmp2=matrix_inner(world,g,g); + return sqrt(tmp1.trace(tmp2)); + } + + std::vector> get_functions(const particle& p) const { + MADNESS_CHECK(p.is_first() or p.is_last()); + if (p.is_first()) return g; + return h; + } + + std::vector> get_g() const {return g;} + std::vector> get_h() const {return h;} + + long rank() const {return g.size();} + + /// return the size in GByte + double size() const { + double sz=get_size(world,g); + sz+=get_size(world,h); + return sz; + } + + Function reconstruct() const { + auto fapprox=hartree_product(g[0],h[0]); + for (int i=1; i> orthonormalize(const std::vector>& g) const { + + double tol=rank_revealing_tol; + std::vector> g2; + auto ovlp=matrix_inner(world,g,g); + if (orthomethod=="canonical") { + tol*=0.01; + print("orthonormalizing with method/tol",orthomethod,tol); + g2=orthonormalize_canonical(g,ovlp,tol); + } else if (orthomethod=="cholesky") { + print("orthonormalizing with method/tol",orthomethod,tol); + g2=orthonormalize_rrcd(g,ovlp,tol); + } + else { + MADNESS_EXCEPTION("no such orthomethod",1); + } + double tight_thresh=FunctionDefaults<3>::get_thresh()*0.1; + return truncate(g2,tight_thresh); + } + + + /// optimize the lrf using the lrfunctor + + /// @param[in] nopt number of iterations (wrt to Alg. 4.3 in Halko) + void optimize(const LRFunctorBase& lrfunctor1, const long nopt=1) { + timer t(world); + t.do_print=do_print; + for (int i=0; i ovlp_g = matrix_inner(world, g, g); + Tensor ovlp_h = matrix_inner(world, h, h); + auto [eval_g, evec_g] = syev(ovlp_g); + auto [eval_h, evec_h] = syev(ovlp_h); + Tensor Xplus=copy(evec_g); + Tensor Xminus=copy(evec_g); + Tensor Yplus=copy(evec_h); + Tensor Yminus=copy(evec_h); + for (int i=0; i M=madness::inner(Xplus,Yplus,0,0); // (X+)^T Y+ + auto [U,s,VT]=svd(M); + + // truncate + typename Tensor::scalar_type s_accumulated=0.0; + int i=s.size()-1; + for (;i>=0; i--) { + s_accumulated+=s[i]; + if (s_accumulated>thresh) { + i++; + break; + } + } + for (int j=0; j XX=madness::inner(Xminus,U,1,0); + Tensor YY=madness::inner(Yminus,VT,1,1); + + g=truncate(transform(world,g,XX)); + h=truncate(transform(world,h,YY)); + } + + + double check_orthonormality(const std::vector>& v) const { + Tensor ovlp=matrix_inner(world,v,v); + return check_orthonormality(ovlp); + } + + double check_orthonormality(const Tensor& ovlp) const { + timer t(world); + t.do_print=do_print; + Tensor ovlp2=ovlp; + for (int i=0; i& lrfunctor1) const { + + timer t(world); + t.do_print=do_print; + + // \int f(1,2)^2 d1d2 + double term1 = lrfunctor1.norm2(); + t.tag("computing term1"); + + // \int f(1,2) pre(1) post(2) \sum_i g(1) h(2) d1d2 +// double term2=madness::inner(pre*g,f12(post*h)); + double term2=madness::inner(g,inner(lrfunctor1,h,p2,p1)); + t.tag("computing term2"); + + // g functions are orthonormal + // \int gh(1,2)^2 d1d2 = \int \sum_{ij} g_i(1) g_j(1) h_i(2) h_j(2) d1d2 + // = \sum_{ij} \int g_i(1) g_j(1) d1 \int h_i(2) h_j(2) d2 + // = \sum_{ij} delta_{ij} \int h_i(2) h_j(2) d2 + // = \sum_{i} \int h_i(2) h_i(2) d2 + double zero=check_orthonormality(g); + if (zero>1.e-10) print("g is not orthonormal",zero); + double term3a=madness::inner(h,h); + auto tmp1=matrix_inner(world,h,h); + auto tmp2=matrix_inner(world,g,g); + double term3=tmp1.trace(tmp2); + print("term3/a/diff",term3a,term3,term3-term3a); + t.tag("computing term3"); + + double arg=term1-2.0*term2+term3; + if (arg<0.0) throw std::runtime_error("negative argument in l2error"); + double error=sqrt(term1-2.0*term2+term3)/sqrt(term1); + if (world.rank()==0 and do_print) { + print("term1,2,3, error",term1, term2, term3, " --",error); + } + + return error; + } + + }; + + // This interface is necessary to compute inner products + template + double inner(const LowRankFunction& a, const LowRankFunction& b) { + World& world=a.world; + return (matrix_inner(world,a.g,b.g).emul(matrix_inner(world,a.h,b.h))).sum(); + } + + + +// template +// LowRankFunction inner(const Function& lhs, const LowRankFunction& rhs, +// const std::tuple v1, const std::tuple v2) { +// World& world=rhs.world; +// // int lhs(1,2) rhs(2,3) d2 = \sum \int lhs(1,2) g_i(2) h_i(3) d2 +// // = \sum \int lhs(1,2) g_i(2) d2 h_i(3) +// LowRankFunction result(world); +// result.h=rhs.h; +// decltype(rhs.g) g; +// for (int i=0; i + LowRankFunction inner(const Function& f1, const LowRankFunction& f2, + const particle p1, const particle p2) { + auto result=inner(f2,f1,p2,p1); + std::swap(result.g,result.h); + return result; + } + + /// lrf(1,3) = inner(lrf(1,2), full(2,3)) + + /// @param[in] f1 the first function + /// @param[in] f2 the second function + /// @param[in] p1 the integration variable of the first function + /// @param[in] p2 the integration variable of the second function + template + LowRankFunction inner(const LowRankFunction& f1, const Function& f2, + const particle p1, const particle p2) { + World& world=f1.world; + static_assert(2*PDIM==NDIM); + // int f(1,2) k(2,3) d2 = \sum \int g_i(1) h_i(2) k(2,3) d2 + // = \sum g_i(1) \int h_i(2) k(2,3) d2 + LowRankFunction result(world); + if (p1.is_last()) { // integrate over 2: result(1,3) = lrf(1,2) f(2,3) + result.g = f1.g; + decltype(f1.h) h; + for (int i=0; i::particle1().get_tuple()),p2.get_tuple())); + result.h=copy(h); + } else if (p1.is_first()) { // integrate over 1: result(2,3) = lrf(1,2) f(1,3) + result.g = f1.h; // correct! second variable of f1 becomes first variable of result + decltype(f1.g) h; + for (int i=0; i::particle1().get_tuple(),p2.get_tuple())); + result.h=copy(h); + } + return result; + } + + /// lrf(1,3) = inner(lrf(1,2), lrf(2,3)) + + /// @param[in] f1 the first function + /// @param[in] f2 the second function + /// @param[in] p1 the integration variable of the first function + /// @param[in] p2 the integration variable of the second function + template + LowRankFunction inner(const LowRankFunction& f1, const LowRankFunction& f2, + const particle p1, const particle p2) { + World& world=f1.world; + static_assert(2*PDIM==NDIM); + + // inner(lrf(1,2) ,lrf(2,3) ) = \sum_ij g1_i(1) h2_j(3) + auto matrix=matrix_inner(world,f2.get_functions(p2),f1.get_functions(p1)); + auto htilde=transform(world,f2.get_functions(p2.complement()),matrix); + auto gg=copy(world,f1.get_functions(p1.complement())); + return LowRankFunction(gg,htilde,f1.rank_revealing_tol,f1.orthomethod); + } + + /// f(1) = inner(lrf(1,2), f(2)) + + /// @param[in] f1 the first function + /// @param[in] vf vector of the second functions + /// @param[in] p1 the integration variable of the first function + /// @param[in] p2 the integration variable of the second function, dummy variable for consistent notation + template + std::vector> inner(const LowRankFunction& f1, const std::vector>& vf, + const particle p1, const particle p2=particle::particle1()) { + World& world=f1.world; + static_assert(2*PDIM==NDIM); + MADNESS_CHECK(p2.is_first()); + + // inner(lrf(1,2), f_k(2) ) = \sum_i g1_i(1) + auto matrix=matrix_inner(world,f1.get_functions(p1),vf); + return transform(world,f1.get_functions(p1.complement()),matrix); + } + + /// f(1) = inner(lrf(1,2), f(2)) + + /// @param[in] f1 the first function + /// @param[in] vf the second function + /// @param[in] p1 the integration variable of the first function + /// @param[in] p2 the integration variable of the second function, dummy variable for consistent notation + template + Function inner(const LowRankFunction& f1, const Function& f2, + const particle p1, const particle p2=particle::particle1()) { + return inner(f1,std::vector>({f2}),p1,p2)[0]; + } + + template + class LowRankFunctionFactory { + public: + + const particle p1=particle::particle1(); + const particle p2=particle::particle2(); + LowRankFunctionParameters parameters; + LowRankFunctionFactory() = default; + LowRankFunctionFactory(const LowRankFunctionParameters param) : parameters(param) {} + LowRankFunctionFactory(const LowRankFunctionFactory& other) = default; + + LowRankFunctionFactory& set_radius(const double radius) { + parameters.set_user_defined_value("radius",radius); + return *this; + } + LowRankFunctionFactory& set_volume_element(const double volume_element) { + parameters.set_user_defined_value("volume_element",volume_element); + return *this; + } + LowRankFunctionFactory& set_rank_revealing_tol(const double rrtol) { + parameters.set_user_defined_value("tol",rrtol); + return *this; + } + LowRankFunctionFactory& set_orthomethod(const std::string orthomethod) { + parameters.set_user_defined_value("orthomethod",orthomethod); + return *this; + } + + LowRankFunction project(const LRFunctorBase& lrfunctor) const { + World& world=lrfunctor.world(); + bool do_print=true; + timer t1(world); + t1.do_print=do_print; + auto orthomethod=parameters.orthomethod(); + auto rank_revealing_tol=parameters.tol(); + + // get sampling grid + randomgrid rgrid(parameters.volume_element(),parameters.radius()); + std::vector> grid=rgrid.get_grid(); + auto Y=Yformer(lrfunctor,grid,parameters.rhsfunctiontype()); + t1.tag("Yforming"); + + auto ovlp=matrix_inner(world,Y,Y); // error in symmetric matrix_inner, use non-symmetric form here! + t1.tag("compute ovlp"); + auto g=truncate(orthonormalize_rrcd(Y,ovlp,rank_revealing_tol)); + t1.tag("rrcd/truncate/thresh"); + auto sz=get_size(world,g); + if (world.rank()==0 and do_print) print("gsize",sz); +// check_orthonormality(g); + + if (world.rank()==0 and do_print) { + print("Y.size()",Y.size()); + print("g.size()",g.size()); + } + + auto h=truncate(inner(lrfunctor,g,p1,p1)); + t1.tag("Y backprojection with truncation"); + return LowRankFunction(g,h,parameters.tol(),parameters.orthomethod()); + + } + + /// apply a rhs (delta or exponential) on grid points to the hi-dim function and form Y = A_ij w_j (in Halko's language) + std::vector> Yformer(const LRFunctorBase& lrfunctor1, const std::vector>& grid, + const std::string rhsfunctiontype, const double exponent=30.0) const { + + World& world=lrfunctor1.world(); + std::vector> Y; + if (rhsfunctiontype=="exponential") { + std::vector> omega; + double coeff=std::pow(2.0*exponent/constants::pi,0.25*LDIM); + for (const auto& point : grid) { + omega.push_back(FunctionFactory(world) + .functor([&point,&exponent,&coeff](const Vector& r) + { + auto r_rel=r-point; + return coeff*exp(-exponent*madness::inner(r_rel,r_rel)); + })); + } + Y=inner(lrfunctor1,omega,p2,p1); + } else { + MADNESS_EXCEPTION("confused rhsfunctiontype",1); + } + auto norms=norm2s(world,Y); + std::vector> Ynormalized; + + for (int i=0; iparameters.tol()) Ynormalized.push_back(Y[i]); + normalize(world,Ynormalized); + return Ynormalized; + } + + }; + + +} // namespace madness + +#endif //MADNESS_LOWRANKFUNCTION_H diff --git a/src/madness/chem/mp2.cc b/src/madness/chem/mp2.cc index bf1cb96b4c1..d726e490f75 100644 --- a/src/madness/chem/mp2.cc +++ b/src/madness/chem/mp2.cc @@ -50,6 +50,8 @@ #include #include #include +#include +#include #include @@ -260,6 +262,7 @@ double MP2::value(const Tensor& x) { if (world.rank() == 0) { printf("current decoupled mp2 energy %12.8f\n", correlation_energy); } + if (param.no_compute()) return correlation_energy; correlation_energy = 0.0; if (hf->get_calc().param.do_localize()) { @@ -280,6 +283,186 @@ double MP2::value(const Tensor& x) { return correlation_energy; } +double MP2::mp3() const { + + typedef std::vector ClusterFunction; + Pairs clusterfunctions; + double mp3_energy=0.0; + CCTimer t1(world,"make pairs"); + // load converged MP1 wave functions + for (int i = param.freeze(); i < hf->nocc(); ++i) { + for (int j = i; j < hf->nocc(); ++j) { +// pairs(i, j) = make_pair(i, j); // initialize + ClusterFunction tmp; + tmp.push_back(CCPairFunction(pairs(i,j).function)); + CCPairFunction ij(hf->nemo(i),hf->nemo(j)); + + CCConvolutionOperator::Parameters cparam; + cparam.thresh_op*=0.1; + auto f12=CCConvolutionOperatorPtr(world,OpType::OT_F12,cparam); + auto vfij=Q12(std::vector({f12*ij})); + for (auto& p : vfij) tmp.push_back(p); + + clusterfunctions(i,j)=tmp; + print("prep pairs",t1.reset()); + } + } + t1.print(); + auto R2 = hf->nemo_ptr->ncf->square(); + + CCTimer t2(world,"recompute MP2"); + // recompute MP2 energy + CCConvolutionOperator::Parameters cparam; + auto g12=CCConvolutionOperatorPtr(world,OpType::OT_G12,cparam); + for (int i = param.freeze(); i < hf->nocc(); ++i) { + for (int j = i; j < hf->nocc(); ++j) { + auto bra=g12*CCPairFunction(hf->R2orbital(i),hf->R2orbital(j)); + double energy1=inner(bra,clusterfunctions(i,j).front()); + double energy=energy1+pairs(i,j).ij_gQf_ij; + print("MP2 energy: gQf, u",pairs(i,j).ij_gQf_ij,energy1,energy); + print("time compute ",t2.reset()); + + double energy2=inner({bra},clusterfunctions(i,j)); + printf("MP2 energy: cluster %12.8f\n",energy2); + print("time clusterfunction",t2.reset()); + } + } + + CCTimer t3(world,"MP3"); + + // compute the term (2) + madness::Pairs gij; + + for (int i = param.freeze(); i < hf->nocc(); ++i) { + for (int j = i; j < hf->nocc(); ++j) { + gij.insert(i,j,(*g12)(hf->nemo(i)*hf->R2orbital(j))); + } + } + print("\n compute term1 of the MP3 energy\n"); + // compute the MP3 energy + double term1=0.0; + if (0) { + for (int i = param.freeze(); i < hf->nocc(); ++i) { + for (int j = i; j < hf->nocc(); ++j) { + auto bra = clusterfunctions(i, j); + double tmp1 = inner(bra, g12 * clusterfunctions(i, j), R2); + double tmp2 = inner(bra, g12 * clusterfunctions(j, i), R2); + double fac = (i == j) ? 0.5 : 1.0; + double tmp = fac * (4.0 * tmp1 - 2.0 * tmp2); + print("mp3 energy: term1", i, j, tmp); + term1 += tmp; + } + } + } + printf("MP3 energy: term1 %12.8f\n",term1); + print("time term1",t3.reset()); + + + // compute intermediates for terms G, I, H, and J + + // \sum_j tau_ij(1,2) * phi_j(2) + std::vector tau_i_jj(hf->nocc()-param.freeze()); + // \sum_j tau_ij(1,2) * phi_j(1) + std::vector tau_ij_j(hf->nocc()-param.freeze()); + for (int i = param.freeze(); i < hf->nocc(); ++i) { + for (int j = param.freeze(); j < hf->nocc(); ++j) { + auto tmp2=multiply(clusterfunctions(i,j),hf->R2orbital(j),{3,4,5}); + for (auto& t : tmp2) tau_i_jj[i].push_back(t); + + auto tmp4=multiply(clusterfunctions(i,j),hf->R2orbital(j),{0,1,2}); + for (auto& t : tmp4) tau_ij_j[i].push_back(t); + } + } + print("info on tau_i_jj and tau_ij_j"); + for (int i = param.freeze(); i < hf->nocc(); ++i) { + for (auto& c: tau_i_jj[i]) c.info(); + } + for (int i = param.freeze(); i < hf->nocc(); ++i) { + for (auto& c: tau_ij_j[i]) c.info(); + } + + print("time GHIJ prep",t3.reset()); + // terms G, I, H, J of Bartlett/Silver 1975 + double term2a=0.0; + real_convolution_3d& g=*(g12->get_op()); + g.particle()=2; + for (int i = param.freeze(); i < hf->nocc(); ++i) { + auto gtau=g(tau_i_jj[i]); + print("info on gtau_i_jj"); + for (auto& g :gtau) g.info(); + double G=inner(gtau,tau_i_jj[i],R2); + print("G",G); +// double G=inner(tau_i_jj[i],g12*tau_i_jj[i],R2); + double I=inner(tau_i_jj[i],g12*tau_ij_j[i],R2); + print("I",I); + double H=inner(tau_ij_j[i],g12*tau_i_jj[i],R2); + print("H",H); + double J=inner(tau_ij_j[i],g12*tau_ij_j[i],R2); + print("J",J); + double tmp = (8.0 *G - 4.0*I + 2.0* H - 4.0*J); + print("mp3 energy: term2",i,tmp); + term2a+=tmp; + } + print("time GHIJ ",t3.reset()); + + double term2=0.0; + for (int i = param.freeze(); i < hf->nocc(); ++i) { + for (int j = param.freeze(); j < hf->nocc(); ++j) { + auto bra_ij=clusterfunctions(i,j); + + double tmp=0.0; + for (int k=param.freeze(); knocc(); ++k) { + auto tau_ij_gik = multiply(bra_ij,gij(i,k),{3,4,5}); + auto tau_ij_gjk = multiply(bra_ij,gij(j,k),{3,4,5}); + + double tmp1=2.0*inner(tau_ij_gik,clusterfunctions(k,j),R2) + - inner(tau_ij_gik,clusterfunctions(j,k),R2) + +2.0* inner(tau_ij_gjk,clusterfunctions(i,k),R2) + - inner(tau_ij_gjk,clusterfunctions(k,i),R2); + tmp-=2.0*tmp1; + } + print("mp3 energy: term2",i,j,tmp); + term2+=tmp; + } + } + print("time term2",t3.reset()); + + double term3=0.0; + for (int i = param.freeze(); i < hf->nocc(); ++i) { + for (int j = param.freeze(); j < hf->nocc(); ++j) { + double tmp=0.0; + for (int k=param.freeze(); knocc(); ++k) { + for (int l=param.freeze(); lnocc(); ++l) { + auto bra = clusterfunctions(i,k); + double ovlp1=inner(bra,clusterfunctions(j,l),R2); + double ovlp2=inner(bra,clusterfunctions(l,j),R2); + auto ket_i=hf->nemo(i); + auto ket_k=hf->nemo(k); + auto bra_j=hf->R2orbital(j); + auto bra_l=hf->R2orbital(l); + + double g_jlik=inner(bra_j*ket_i, (*g12)(bra_l*ket_k)); + double g_jlki=inner(bra_j*ket_k, (*g12)(bra_l*ket_i)); + tmp+=(2.0*ovlp1 - ovlp2)*g_jlik; + } + } + print("mp3 energy: term3",i,j,tmp); + term3+=tmp; + } + } + print("time term3",t3.reset()); + t3.print(); + printf("term1 %12.8f\n",term1); + printf("term2a %12.8f\n",term2a); + printf("term2 %12.8f\n",term2); + printf("term3 %12.8f\n",term3); + mp3_energy=term1+term2a+term2+term3; + printf("MP3 en %12.8f\n",mp3_energy); + + + return mp3_energy; +} + /// solve the residual equation for electron pair (i,j) void MP2::solve_residual_equations(ElectronPair& result, const double econv, const double dconv) const { @@ -525,13 +708,6 @@ double MP2::asymmetry(const real_function_6d& f, const std::string s) const { return diff; } -void MP2::test(const std::string filename) { - if (world.rank() == 0) - printf("starting coupling at time %8.1fs\n", wall_time()); - if (world.rank() == 0) - printf("ending coupling at time %8.1fs\n", wall_time()); -} - /// compute the matrix element /// scales quartically. I think I can get this down to cubically by @@ -540,6 +716,7 @@ void MP2::test(const std::string filename) { /// @return the energy double MP2::compute_gQf(const int i, const int j, ElectronPair& pair) const { + CCTimer t1(world,"gQf old"); // for clarity of notation const int k = pair.i; const int l = pair.j; @@ -619,6 +796,33 @@ double MP2::compute_gQf(const int i, const int j, ElectronPair& pair) const { if (world.rank() == 0) printf("<%d%d | g Q12 f | %d%d> %12.8f\n", i, j, k, l, e); + print("gQf old",t1.reset()); + + CCTimer timer(world,"gQf with ccpairfunction"); + CCConvolutionOperator::Parameters param; + auto f=CCConvolutionOperatorPtr(world,OpType::OT_F12,param); + auto g=CCConvolutionOperatorPtr(world,OpType::OT_G12,param); +// print("operator constructor",timer.reset()); + + CCPairFunction fij(f,ket_i,ket_j); + CCPairFunction gij(g,bra_k,bra_l); +// CCPairFunction ij({psi},{psi}); + std::vector vfij={fij}; + std::vector vgij={gij}; +// std::vector vij={ij}; +// print("g/f ij constructor",timer.reset()); + +// StrongOrthogonalityProjector SO(world); +// SO.set_spaces({psi},{psi},{psi},{psi}); + std::vector Qfij=Q12(vfij); +// std::vector Qgij=Q12(vgij); +// print("SO application",timer.reset()); + + double result1=inner(vgij,Qfij); + print("inner",timer.reset()); + + print(")",result1); + return e; } diff --git a/src/madness/chem/mp2.h b/src/madness/chem/mp2.h index 4d0d9472b15..b60844703f5 100644 --- a/src/madness/chem/mp2.h +++ b/src/madness/chem/mp2.h @@ -329,6 +329,7 @@ class MP2 : public OptimizationTargetInterface, public QCPropertyInterface { initialize < int > ("freeze", 0); initialize < int > ("maxsub", 2); initialize < bool > ("restart", true); + initialize < bool > ("no_compute", false); initialize < int > ("maxiter", 5); } @@ -362,6 +363,7 @@ class MP2 : public OptimizationTargetInterface, public QCPropertyInterface { int i() const { return this->get >("pair")[0]; } /// convenience function int j() const { return this->get >("pair")[1]; } /// convenience function int restart() const { return this->get("restart"); } /// convenience function + int no_compute() const { return this->get("no_compute"); } /// convenience function int maxiter() const { return this->get("maxiter"); } /// convenience function int maxsub() const { return this->get("maxsub"); } /// convenience function bool do_oep() const { return do_oep1;} @@ -460,6 +462,8 @@ class MP2 : public OptimizationTargetInterface, public QCPropertyInterface { return hf->orbital_energy(i) + hf->orbital_energy(j); } + double mp3() const; + /// solve the residual equation for electron pair (i,j) /// \todo Parameter documentation. Below are un-doxygenated comments that no longer seem relevant? @@ -503,8 +507,6 @@ class MP2 : public OptimizationTargetInterface, public QCPropertyInterface { double asymmetry(const real_function_6d& f, const std::string s) const; - void test(const std::string filename); - /// compute the matrix element /// scales quartically. I think I can get this down to cubically by diff --git a/src/madness/chem/test_ccpairfunction.cc b/src/madness/chem/test_ccpairfunction.cc index 505c0828c29..e035e5803b0 100644 --- a/src/madness/chem/test_ccpairfunction.cc +++ b/src/madness/chem/test_ccpairfunction.cc @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -15,15 +16,16 @@ using namespace madness; +bool longtest=false; struct data { real_function_3d f1,f2,f3,f4,f5; - real_function_6d f; + real_function_6d f12,f23; - std::shared_ptr f12; + std::shared_ptr f12_op; data() {} data(World& world, CCParameters& parameters) { - f12.reset(new CCConvolutionOperator(world,OT_F12,parameters)); + f12_op.reset(new CCConvolutionOperator(world,OT_F12,parameters)); if (not f1.is_initialized()) initialize(world); } @@ -44,38 +46,93 @@ struct data { double r2=r[3]*r[3] + r[4]*r[4] + r[5]*r[5]; return exp(-1.0*r1 - 2.0*r2); }; - f = real_factory_6d(world).f(g); + auto g23 = [](const coord_6d& r) { + double r1=r[0]*r[0] + r[1]*r[1] + r[2]*r[2]; + double r2=r[3]*r[3] + r[4]*r[4] + r[5]*r[5]; + return exp(-1.0*r1 - 2.0*r2) + exp(-2.0*r1 - 3.0*r2); + }; + f12 = real_factory_6d(world).f(g); + f23 = real_factory_6d(world).f(g23); } void clear() { - f12.reset(); + f12_op.reset(); f1.clear(); f2.clear(); f3.clear(); f4.clear(); f5.clear(); - f.clear(); + f12.clear(); + f23.clear(); } + /// get some standard functions + + /// f1: exp(-1.0 r^2) + /// f2: exp(-2.0 r^2) + /// f3: exp(-3.0 r^2) + /// f4: exp(-4.0 r^2) + /// f5: exp(-5.0 r^2) + /// f12: exp(-r_1^2 - 2 r_2^2) + /// f23: exp(-r_1^2 - 2 r_2^2) + exp(-2 r_1^2 - 3 r_2^2) auto get_functions() const { - return std::make_tuple(f1,f2,f3,f4,f5,f); + return std::make_tuple(f1,f2,f3,f4,f5,f12); } + /// get some standard ccpairfunctions + + /// p1: pure, corresponds to f12 + /// p2: dec, corresponds to f23 + /// p3: op_dec, corresponds to f23 + /// p4: pure, corresponds to f23 + /// p5: op_pure, corresponds to f23 auto get_ccpairfunctions() const { - CCPairFunction p1; - CCPairFunction p2(f); - CCPairFunction p3({f1,f2},{f1,f3}); - CCPairFunction p4(f12.get(),{f1,f2},{f1,f3}); - return std::make_tuple(p1,p2,p3); + CCPairFunction p1(f12); + CCPairFunction p2({f1,f2},{f2,f3}); + CCPairFunction p3(f12_op,{f1,f2},{f2,f3}); + CCPairFunction p4(f23); // two-term, corresponds to p2 + CCPairFunction p5(f12_op,f23); // two-term, corresponds to p2 + return std::make_tuple(p1,p2,p3,p4,p5); } }; data data1; +/// Computes the electrostatic potential due to a Gaussian charge distribution + +/// stolen from testsuite.cc +class GaussianPotential : public FunctionFunctorInterface { +public: + typedef Vector coordT; + const coordT center; + const double exponent; + const double coefficient; + + GaussianPotential(const coordT& center, double expnt, double coefficient) + : center(center) + , exponent(sqrt(expnt)) + , coefficient(coefficient*pow(constants::pi/exponent,1.5)*pow(expnt,-0.75)) {} + + double operator()(const coordT& x) const { + double sum = 00; + for (int i=0; i<3; ++i) { + double xx = center[i]-x[i]; + sum += xx*xx; + }; + double r = sqrt(sum); + if (r<1.e-4) { // correct thru order r^3 + const double sqrtpi=sqrt(constants::pi); + const double a=exponent; + return coefficient*(2.0*a/sqrtpi - 2.0*a*a*a*r*r/(3.0*sqrtpi)); + } else { + return coefficient*erf(exponent*r)/r; + } + } +}; + int test_constructor(World& world, std::shared_ptr ncf, const Molecule& molecule, - const CCParameters& parameter) { - int success=0; + const CCParameters& parameter) { test_output t1("CCPairFunction::constructor"); real_function_6d f=real_factory_6d(world); @@ -83,13 +140,13 @@ int test_constructor(World& world, std::shared_ptr ncf vector_real_function_3d a= zero_functions(world,3); vector_real_function_3d b= zero_functions(world,3); - CCConvolutionOperator f12(world, OT_F12, parameter); + auto f12=CCConvolutionOperatorPtr(world, OT_F12, parameter); t1.checkpoint(true,"preparation"); CCPairFunction p1; CCPairFunction p2(f); CCPairFunction p3({f1,f2},{f1,f3}); - CCPairFunction p4(&f12,{f1,f2},{f2,f3}); + CCPairFunction p4(f12,{f1,f2},{f2,f3}); t1.checkpoint(true,"construction"); { @@ -133,203 +190,372 @@ int test_constructor(World& world, std::shared_ptr ncf } t1.checkpoint(true,"checks on assignment"); - t1.end(); - return (t1.get_final_success()) ? 0 : 1; + return t1.end(); +} + +int test_operator_apply(World& world, std::shared_ptr ncf, const Molecule& molecule, + const CCParameters& parameter) { + test_output t1("CCPairFunction::test_operator_apply"); +// t1.set_cout_to_terminal(); + + double exponent=1.0; // corresponds to the exponent of data::f1 and data::ff +// double coefficient=pow(1.0/constants::pi*exponent,0.5*3); + double coefficient=1.0; + const coord_3d center={0.0,0.0,-0.0}; + + auto Gaussian = [¢er, &exponent, &coefficient](const coord_3d& r) { + return coefficient * exp(-exponent*inner(r-center,r-center)); + }; + + // this won't work, the simulation cell is only [-1,1].. + // Normalized Gaussian exponent a produces potential erf(sqrt(a)*r)/r +// GaussianPotential gpot(center, exponent, coefficient); +// real_function_3d op_a=real_factory_3d(world).functor(gpot); + + real_function_3d a=real_factory_3d(world).functor(Gaussian); + exponent=2.0; + real_function_3d b=real_factory_3d(world).functor(Gaussian); + + auto [f1,f2,f3,f4,f5,ff]=data1.get_functions(); + CCPairFunction c1(a,b); + CCPairFunction c2(f1,f2); +// CCPairFunction ref(op_a,b); + auto op= CoulombOperator(world,1.e-5,FunctionDefaults<3>::get_thresh()); + op.print_timings=false; + op.particle()=1; + + auto op_c1=op(c1); + auto op_c2=op(c2); +// double a1=inner(ref,ref); + double norm1=inner(op_c1,op_c1); + double norm2=inner(op_c2,op_c2); + print("norm1,norm2",norm1,norm2); + bool good=fabs(norm1-norm2)::get_thresh(); + t1.checkpoint(good,"op(xx)"); + return t1.end(); +} + + +int test_transformations(World& world, std::shared_ptr ncf, const Molecule& molecule, + const CCParameters& parameter) { + test_output t1("CCPairFunction::test_transformations"); + + real_function_6d f=real_factory_6d(world); + auto [f1,f2,f3,f4,f5,ff]=data1.get_functions(); + auto f12=CCConvolutionOperatorPtr(world, OT_F12, parameter); + auto g12=CCConvolutionOperatorPtr(world, OT_G12, parameter); + + + CCPairFunction p1(ff); + t1.checkpoint(p1.is_pure(),"is_pure"); + t1.checkpoint(p1.is_pure_no_op(),"is_pure_no_op"); + + CCPairFunction p2(f12,ff); + t1.checkpoint(p2.is_pure(),"is_pure"); + t1.checkpoint(p2.is_op_pure(),"is_op_pure"); + t1.checkpoint(p2.is_convertible_to_pure_no_op(),"is_convertible_to_pure_no_op"); + CCPairFunction p3=copy(p2); + p3.convert_to_pure_no_op_inplace(); + t1.checkpoint(p2.is_op_pure(),"is_op_pure"); + t1.checkpoint(p3.is_pure_no_op(),"is_pure_no_op"); + + CCPairFunction p4(g12,ff); + t1.checkpoint(p4.is_pure(),"is_pure"); + t1.checkpoint(p4.is_op_pure(),"is_op_pure"); + t1.checkpoint(not p4.is_convertible_to_pure_no_op(),"not is_convertible_to_pure_no_op"); + + CCPairFunction p5(f12,f1,f2); + t1.checkpoint(not p5.is_pure(),"is_pure"); + t1.checkpoint(p5.is_op_decomposed(),"is_op_decomposed"); + t1.checkpoint(p5.is_convertible_to_pure_no_op(),"is_convertible_to_pure_no_op"); + CCPairFunction p6=copy(p5); + p6.convert_to_pure_no_op_inplace(); + t1.checkpoint(p6.is_pure_no_op(),"is_pure_no_op"); + + return t1.end(); } -int test_overlap(World& world, std::shared_ptr ncf, const Molecule& molecule, +int test_multiply_with_f12(World& world, std::shared_ptr ncf, const Molecule& molecule, + const CCParameters& parameters) { + test_output t1("CCPairFunction::test_multiply_with_f12"); + + // p1: pure, corresponds to f12 + // p2: dec, corresponds to f23 + // p3: op_dec, corresponds to f23 + // p4: pure, corresponds to f23 + // p5: op_pure, corresponds to f23 + auto [p1,p2,p3,p4,p5]=data1.get_ccpairfunctions(); // p2-p5 correspond to f23 + auto f12=data1.f12_op; + + double thresh=FunctionDefaults<3>::get_thresh(); + + // decomposed + CCPairFunction tmp1=f12*p2; // should now be identical to p3 + CCPairFunction tmp2=p2*f12; // should now be identical to p3 + double ref=inner(p2,p3); + + double r1=inner(p2,tmp1); + bool good=(fabs(ref-r1) ncf, const Molecule& molecule, const CCParameters& parameters) { - CCTimer timer(world, "testing"); - auto R2 = ncf->square(); + test_output t1("CCPairFunction::test_multiply"); - CCConvolutionOperator g12(world, OT_G12, parameters); - CCConvolutionOperator f12(world, OT_F12, parameters); + // consistency check, relies on CCPairFunction::inner to work correctly - CorrelationFactor corrfac(world, 1.0, 1.e-7, molecule); + double thresh=FunctionDefaults<3>::get_thresh(); + auto [p1,p2,p3,p4,p5]=data1.get_ccpairfunctions(); // p2-p5 correspond to f23 + auto [f1,f2,f3,f4,f5,f]=data1.get_functions(); - auto g = [](const coord_3d& r) { return exp(-1.0 * r.normf()); }; - real_function_3d f1 = real_factory_3d(world).f(g); - real_function_3d R2f = R2 * f1; - double norm = inner(f1, R2f); - f1.scale(1.0 / sqrt(norm)); - R2f.scale(1.0 / sqrt(norm)); - CC_vecfunction mo_ket_(std::vector(1, f1), HOLE); - CC_vecfunction mo_bra_(std::vector(1, R2f), HOLE); + // reference value is = + CCPairFunction bra(f1,f2); + CCPairFunction bra1(f1*f2,f2); + CCPairFunction bra2(f1,f2*f2); + for (auto& p : {p2,p3,p4,p5}) { - bool passed_lo = true; - bool passed_hi = true; - bool success=true; - int isuccess=0; - const double lo = parameters.thresh_6D(); - const double hi = parameters.thresh_3D(); - const double hi_loose = parameters.thresh_3D()*5.0; + auto tmp1=multiply(p,f2,{0,1,2}); + double ovlp1=inner(bra,tmp1); + double ref1=p.has_operator() ? inner(bra1,p3) : inner(bra1,p2); - print("threshold 3D / 6D", hi,lo); - test_output t1("CCPairFunction::inner"); + bool good=(fabs(ovlp1-ref1) ncf, const Molecule& molecule, + const CCParameters& parameters) { + test_output t1("CCPairFunction::test_inner"); t1.set_cout_to_terminal(); - // f^2 = 1/(4y^2)(1 - 2*f2(y) + f2(2y)) , f2(2y) =f2(y)^2 + + /// f1: exp(-1.0 r^2) + /// f2: exp(-2.0 r^2) + /// f3: exp(-3.0 r^2) + /// f4: exp(-4.0 r^2) + /// f5: exp(-5.0 r^2) + /// f: exp(-r_1^2 - 2 r_2^2) + /// f23: exp(-r_1^2 - 2 r_2^2) + exp(-2 r_1^2 - 3 r_2^2) + auto [f1,f2,f3,f4,f5,f]=data1.get_functions(); + /// p1: pure, corresponds to f12 + /// p2: dec, corresponds to f23 + /// p3: op_dec, corresponds to f23 + /// p4: pure, corresponds to f23 + /// p5: op_pure, corresponds to f23 + auto [p1,p2,p3,p4,p5]=data1.get_ccpairfunctions(); + auto f12 = *(data1.f12_op); + + /// results + auto a=std::vector({f1,f2}); + auto b=std::vector({f2,f3}); + std::vector a_ij_functions, b_ij_functions; + for (int i=0; i = \sum_{ij} + double ab_ab=aij.trace(bij); + + // = \sum_{ij} < _1(2) | b_ib_j(2) > + double ab_f_ab=dot(world,f12(a_ij_functions),b_ij_functions).trace(); + + // = \sum_{ij} < _2 | b_ib_j(2) > + // f^2 = 1/(4y^2)(1 - 2*f(y) + f2(2y)) , f2(2y) =f2(y)^2 // operator apply of SlaterF12Operator includes a factor of 1/(2 gamma) and the identity // operator apply of SlaterOperator has no further terms const double y = parameters.gamma(); - SeparatedConvolution f = SlaterF12Operator(world, y, parameters.lo(), parameters.thresh_bsh_3D()); - SeparatedConvolution f2 = SlaterOperator(world, y, parameters.lo(), parameters.thresh_bsh_3D()); - SeparatedConvolution ff = SlaterOperator(world, 2.0 * y, parameters.lo(), parameters.thresh_bsh_3D()); - - real_function_3d a = mo_ket_(0).function; - real_function_3d b = mo_ket_(0).function; - real_function_6d fab_6d = CompositeFactory(world).g12(corrfac.f()).particle1(copy(a)).particle2( - copy(b)); - fab_6d.fill_tree().truncate().reduce_rank(); - fab_6d.print_size("fab_6d"); - - real_function_3d aR = mo_bra_(0).function; - real_function_3d bR = mo_bra_(0).function; - - const real_function_3d aa = (aR * a).truncate(); - const real_function_3d bb = (bR * b).truncate(); - const real_function_3d af2a = f2(aa); - const real_function_3d affa = ff(aa); - const real_function_3d afa = f(aa); + SeparatedConvolution fop= SlaterOperator(world, y, parameters.lo(), parameters.thresh_bsh_3D()); + SeparatedConvolution fsq = SlaterOperator(world, 2.0 * y, parameters.lo(), parameters.thresh_bsh_3D()); const double prefactor = 1.0 / (4 * y * y); - const double term1= prefactor * (aR.inner(a) * bR.inner(b)); - const double term2= prefactor * 2.0 * bb.inner(af2a) ; - const double term3= prefactor * bb.inner(affa); - const double ab_f2_ab = prefactor * (aR.inner(a) * bR.inner(b) - 2.0 * bb.inner(af2a) + bb.inner(affa)); - const double ab_f_ab = inner(f(aa), bb); - - const long rank = world.rank(); - auto printer = [&rank](std::string msg, const double result, const double reference, const double time) { - if (rank == 0) { - long len = std::max(0l, long(40 - msg.size())); - msg += std::string(len, ' '); - std::cout << msg << std::fixed << std::setprecision(8) << "result, ref, diff " - << result << " " << reference << " " << std::setw(9) << std::setprecision(2) << std::scientific - << result - reference - << ", elapsed time " << std::fixed << std::setw(5) << std::setprecision(2) << time << std::endl; + const double ab_f2_ab = prefactor*( ab_ab + -2.0*dot(world,apply(world,fop,a_ij_functions),b_ij_functions).trace() + +dot(world,apply(world,fsq,a_ij_functions),b_ij_functions).trace() ); + + + for (auto& ket : {p2, p3, p4, p5}) { + for (auto& bra : {p2, p3, p4, p5}) { + double ref=0.0; + if (bra.has_operator() and ket.has_operator()) ref=ab_f2_ab; + if (bra.has_operator() and (not ket.has_operator())) ref=ab_f_ab; + if ((not bra.has_operator()) and ket.has_operator()) ref=ab_f_ab; + if ((not bra.has_operator()) and (not ket.has_operator())) ref=ab_ab; + double result=inner(bra,ket); + + print(bra.name(true)+ket.name(),"ref, result, diff", ref, result, ref-result); + double thresh=FunctionDefaults<3>::get_thresh(); + bool good=(fabs(result-ref) lo) passed_lo = false; - if (fabs(diff) > hi) passed_hi = false; - - printer(" op_dec/op_dec : ", test1, ab_f2_ab, timer.reset()); - success=(fabs(diff) < hi); - } - if (not success) isuccess++; - t1.checkpoint(success, "op_dec/op_dec"); - { - CCPairFunction ab(mo_ket_.get_vecfunction(), mo_ket_.get_vecfunction()); - const double test1 = inner(ab, ab, R2); - const double test2 = double(mo_ket_.size()); // mos are normed - const double diff = test1 - test2; - printer(" dec/dec : ", test1, double(mo_ket_.size()), timer.reset()); - if (fabs(diff) > lo) passed_lo = false; - if (fabs(diff) > hi) passed_hi = false; - success=(fabs(diff) < hi); } - if (not success) isuccess++; - t1.checkpoint(success, "dec/dec"); - { - CCPairFunction fab(fab_6d); - const double test1 = inner(fab, fab, R2); - const double diff = test1 - ab_f2_ab; - printer(" pure/pure : ", test1, ab_f2_ab, timer.reset()); - if (fabs(diff) > lo) passed_lo = false; - if (fabs(diff) > hi) passed_hi = false; - success=(fabs(diff) < hi_loose); - } - if (not success) isuccess++; - t1.checkpoint(success, "pure/pure"); + return t1.end(); +} - { - CCPairFunction ab1(mo_ket_.get_vecfunction(), mo_ket_.get_vecfunction()); - CCPairFunction fab(&f12, a, b); - timer.reset(); - const double test1 = inner(ab1, fab, R2); - const double test2 = ab_f_ab; - const double diff = test1 - test2; - printer(" dec/op_dec : ", test1, ab_f_ab, timer.reset()); - if (fabs(diff) > lo) passed_lo = false; - if (fabs(diff) > hi) passed_hi = false; - success=(fabs(diff) < hi); - } - if (not success) isuccess++; - t1.checkpoint(success, "dec/op_dec"); - { - // the next tests evaulate in different ways - CCPairFunction fab(fab_6d); - CCPairFunction ab2(mo_ket_.get_vecfunction(), mo_ket_.get_vecfunction()); - const double test1 = inner(fab, ab2, R2); - const double test2 = bb.inner(afa); - const double diff = test1 - test2; - printer(" dec/pure", test1, test2, timer.reset()); - if (fabs(diff) > lo) passed_lo = false; - if (fabs(diff) > hi) passed_hi = false; - - success=(fabs(diff) < hi_loose); - } - if (not success) isuccess++; - t1.checkpoint(success, "dec/pure"); -// { -// CCPairFunction fab(fab_6d); -// CCPairFunction ab2(ab_6d); -// const double test1 = overlap(fab, ab2); -// const double test2 = bb.inner(afa); -// const double diff = test1 - test2; -// if (world.rank() == 0) -// std::cout << "Overlap Test 6 : " << std::fixed << std::setprecision(10) << "result=" << test1 -// << ", test=" << test2 << ", diff=" << diff << "\n"; -// if (fabs(diff) > lo) passed_lo = false; -// if (fabs(diff) > hi) passed_hi = false; -// } - { - CCPairFunction fab(&f12, a, b); - CCPairFunction ab2(fab_6d); - const double test1 = inner(fab, ab2, R2); - const double test2 = ab_f2_ab; - const double diff = test1 - test2; - printer(" op_dec/pure : ", test1, test2, timer.reset()); - if (fabs(diff) > lo) passed_lo = false; - if (fabs(diff) > hi) passed_hi = false; - success=(fabs(diff) < hi_loose); // be a bit loose here .. +int test_partial_inner_6d(World& world, std::shared_ptr ncf, const Molecule& molecule, + const CCParameters& parameter) { + + CCTimer timer(world, "testing"); + test_output t1("CCPairFunction::test_partial_inner_6d"); + t1.set_cout_to_terminal(); + auto [f1,f2,f3,f4,f5,f] = data1.get_functions(); + + std::vector a = {f1, f2}; + std::vector b = {f2, f3}; + + auto f12=CCConvolutionOperatorPtr(world, OT_F12, parameter); + + auto [p1,p2,p3,p4,nil]=data1.get_ccpairfunctions(); + CCPairFunction p11({f1},{f1}); + CCPairFunction p12({f1},{f2}); + + double g11=inner(f1,f1); + double g22=inner(f2,f2); + double g12=inner(f1,f2); + double g13=inner(f1,f3); + double g23=inner(f2,f3); + double g33=inner(f3,f3); + real_function_3d gf11=(*f12)(f1*f1); + real_function_3d gf12=(*f12)(f1*f2); + real_function_3d gf22=(*f12)(f2*f2); + real_function_3d gf13=(*f12)(f1*f3); + real_function_3d gf23=(*f12)(f2*f3); + + t1.checkpoint(true,"prep",timer.reset()); + + // p1 = p12 = e(-r1 -2r2) + // p2 = e(-r1) * e(-2r2) + e(-2r1) * e(-3e2) separated + // p4 = e(-r1) * e(-2r2) + e(-2r1) * e(-3e2) 6d + for (auto test_p1 : {p2,p4}) { + for (auto test_p2 : {p2,p4}) { + CCPairFunction r1=inner(test_p1,test_p2,{0,1,2},{0,1,2}); + CCPairFunction r2=inner(test_p1,test_p2,{0,1,2},{3,4,5}); + CCPairFunction r3=inner(test_p1,test_p2,{3,4,5},{0,1,2}); + CCPairFunction r4=inner(test_p1,test_p2,{3,4,5},{3,4,5}); + + double n1=inner(r1,p11); + double n2=inner(r2,p11); + double n3=inner(r3,p11); + double n4=inner(r4,p11); + + double ref_n1=g11*g12*g12 + g12*g12*g13 + g12*g13*g12 + g22*g13*g13; + double ref_n2=g12*g12*g11 + g13*g12*g12 + g22*g13*g11 + g23*g13*g12; + double ref_n3=ref_n2; + double ref_n4=g22*g11*g11 + g23*g11*g12 + g23*g12*g11 + g33*g12*g12; + + bool good=fabs(n1-ref_n1)::get_thresh(); + t1.checkpoint(good,test_p1.name(true)+test_p2.name()+" -- 1"); + good=fabs(n2-ref_n2)::get_thresh(); + t1.checkpoint(good,test_p1.name(true)+test_p2.name()+" -- 2"); + good=fabs(n3-ref_n3)::get_thresh(); + t1.checkpoint(good,test_p1.name(true)+test_p2.name()+" -- 3"); + good=fabs(n4-ref_n4)::get_thresh(); + t1.checkpoint(good,test_p1.name(true)+test_p2.name()+" -- 4",timer.reset()); + + } } - if (not success) isuccess++; - t1.checkpoint(success, "op_dec/pure"); - { - CCPairFunction fab(&f12, a, b); - CCPairFunction ab2(mo_ket_.get_vecfunction(), mo_ket_.get_vecfunction()); - timer.reset(); - const double test1 = inner(fab, ab2, R2); -// const double test2 = bb.inner(afa); - const double test2 = ab_f_ab; - const double diff = test1 - test2; - printer(" op_dec/dec : ", test1, test2, timer.reset()); - if (fabs(diff) > lo) passed_lo = false; - if (fabs(diff) > hi) passed_hi = false; - success=(fabs(diff) < hi); + // test < sth | f(1,2) a(1)b(2) > + // CCPairFunction p2({f1,f2},{f2,f3}); + // CCPairFunction p4(f23); // two-term, corresponds to p2 + CCPairFunction p5(f12,std::vector({f1}),std::vector({f2})); + for (auto& test_p1 : {p2, p4}) { + CCPairFunction r1=inner(test_p1,p5,{0,1,2},{0,1,2}); + CCPairFunction r2=inner(test_p1,p5,{0,1,2},{3,4,5}); + CCPairFunction r3=inner(test_p1,p5,{3,4,5},{0,1,2}); + CCPairFunction r4=inner(test_p1,p5,{3,4,5},{3,4,5}); + + double ref_n1=inner(gf11,f2*f1) * g12 + inner(gf12,f1*f2) * g13; + double ref_n2=inner(gf12,f1*f1) * g12 + inner(gf22,f1*f1) * g13; + double ref_n3=inner(gf12,f2*f1) * g11 + inner(gf13,f1*f2) * g12; + double ref_n4=inner(gf22,f1*f1) * g11 + inner(gf23,f1*f1) * g12; + + double n1=inner(r1,p11); + double n2=inner(r2,p11); + double n3=inner(r3,p11); + double n4=inner(r4,p11); +// print("n1, ref_n1",n1,ref_n1, n1-ref_n1); +// print("n2, ref_n2",n2,ref_n2, n2-ref_n2); +// print("n3, ref_n3",n3,ref_n3, n3-ref_n3); +// print("n4, ref_n4",n4,ref_n4, n4-ref_n4); + + bool good=fabs(n1-ref_n1)::get_thresh(); + t1.checkpoint(good,test_p1.name(true)+p5.name()+" -- 1"); + good=fabs(n2-ref_n2)::get_thresh(); + t1.checkpoint(good,test_p1.name(true)+p5.name()+" -- 2"); + good=fabs(n3-ref_n3)::get_thresh(); + t1.checkpoint(good,test_p1.name(true)+p5.name()+" -- 3"); + good=fabs(n4-ref_n4)::get_thresh(); + t1.checkpoint(good,test_p1.name(true)+p5.name()+" -- 4",timer.reset()); + } - if (not success) isuccess++; - t1.checkpoint(success, "op_dec/dec"); - return (t1.get_final_success()) ? 0 : 1; -} + // test < a(1) b(2) f(1,2) | f(1,3) c(1)d(3) > + // CCPairFunction p3(f12_op.get(),{f1,f2},{f2,f3}); + // CCPairFunction p5(&f12,{f1},{f2}); + if (longtest) { + CCPairFunction r1=inner(p3,p5,{0,1,2},{0,1,2}); + print("time after r1 ", timer.reset()); + double n1=inner(r1,p11); + p3.convert_to_pure_no_op_inplace(); + p5.convert_to_pure_no_op_inplace(); + print("n1",n1); + print("time after r2a ", timer.reset()); + CCPairFunction r1a=inner(p3,p5,{0,1,2},{0,1,2}); + double n1a=inner(r1a,p11); + print("n1a",n1a); + print("diff",n1-n1a); + CCPairFunction r2=inner(p3,p5,{0,1,2},{3,4,5}); + print("time after r2 ", timer.reset()); + CCPairFunction r3=inner(p3,p5,{3,4,5},{0,1,2}); + print("time after r3 ", timer.reset()); + CCPairFunction r4=inner(p3,p5,{3,4,5},{3,4,5}); + print("time after r4 ", timer.reset()); + } -int test_partial_inner(World& world, std::shared_ptr ncf, const Molecule& molecule, - const CCParameters& parameter) { - int success=0; + return (t1.get_final_success()) ? 0 : 1; +} +int test_partial_inner_3d(World& world, std::shared_ptr ncf, const Molecule& molecule, + const CCParameters& parameter) { CCTimer timer(world, "testing"); test_output t1("CCPairFunction::test_partial_inner"); @@ -338,41 +564,43 @@ int test_partial_inner(World& world, std::shared_ptr n std::vector a = {f1, f2}; std::vector b = {f2, f3}; - CCConvolutionOperator f12(world, OT_F12, parameter); + auto f12=CCConvolutionOperatorPtr(world, OT_F12, parameter); - CCPairFunction p1(f); + CCPairFunction p1(f); // e(-r1 - 2r2) CCPairFunction p2(a,b); - CCPairFunction p3(&f12,a,b); + CCPairFunction p3(f12,a,b); + CCPairFunction p11({f1},{f1}); + CCPairFunction p12({f1},{f2}); double g11=inner(f1,f1); double g22=inner(f2,f2); double g12=inner(f1,f2); double g13=inner(f1,f3); double g23=inner(f2,f3); - real_function_3d gf11=f12(f1*f1); - real_function_3d gf12=f12(f1*f2); - real_function_3d gf13=f12(f1*f3); + real_function_3d gf11=(*f12)(f1*f1); + real_function_3d gf12=(*f12)(f1*f2); + real_function_3d gf13=(*f12)(f1*f3); print("g11, g22",g11,g22); print("time in preparation",timer.reset()); t1.checkpoint(true,"prep"); - // test pure + // test pure/3d { real_function_3d r=inner(p1,f1,{0,1,2},{0,1,2}); double norm=inner(r,f2); double ref_norm=g11 * g22; print("norm, ref_norm",norm,ref_norm); - bool good=fabs(norm-ref_norm::get_thresh()); + bool good=fabs(norm-ref_norm)::get_thresh(); t1.checkpoint(good,"pure -- 1"); } - // test pure + // test pure/3d { real_function_3d r=inner(p1,f1,{3,4,5},{0,1,2}); double norm=inner(r,f2); double ref_norm=g12 * g12; print("norm, ref_norm",norm,ref_norm); - bool good=fabs(norm-ref_norm::get_thresh()); + bool good=fabs(norm-ref_norm)::get_thresh(); t1.checkpoint(good,"pure -- 2"); } // test decomposed @@ -381,7 +609,7 @@ int test_partial_inner(World& world, std::shared_ptr n double norm=inner(r,f2); double ref_norm=g11 * g22 + g12 * g23; print("norm, ref_norm",norm,ref_norm); - bool good=fabs(norm-ref_norm::get_thresh()); + bool good=fabs(norm-ref_norm)::get_thresh(); t1.checkpoint(good,"decomposed -- 1"); } // test decomposed @@ -390,7 +618,7 @@ int test_partial_inner(World& world, std::shared_ptr n double norm=inner(r,f2); double ref_norm=g12 * g12 + g13 * g22; print("norm, ref_norm",norm,ref_norm); - bool good=fabs(norm-ref_norm::get_thresh()); + bool good=fabs(norm-ref_norm)::get_thresh(); t1.checkpoint(good,"decomposed -- 2"); } // test op_decomposed @@ -400,7 +628,7 @@ int test_partial_inner(World& world, std::shared_ptr n double norm=inner(r,f2); double ref_norm=inner(gf11*f2,f2) + inner(gf12*f2,f3); print("norm, ref_norm",norm,ref_norm); - bool good=fabs(norm-ref_norm::get_thresh()); + bool good=fabs(norm-ref_norm)::get_thresh(); t1.checkpoint(good,"op_decomposed -- 1"); } // test op_decomposed @@ -410,7 +638,7 @@ int test_partial_inner(World& world, std::shared_ptr n double norm=inner(r,f2); double ref_norm=inner(gf12*f1,f2) + inner(gf13*f2,f2); print("norm, ref_norm",norm,ref_norm); - bool good=fabs(norm-ref_norm::get_thresh()); + bool good=fabs(norm-ref_norm)::get_thresh(); t1.checkpoint(good,"op_decomposed -- 2"); } @@ -420,15 +648,13 @@ int test_partial_inner(World& world, std::shared_ptr n int test_apply(World& world, std::shared_ptr ncf, const Molecule& molecule, const CCParameters& parameter) { - int success=0; test_output t1("CCPairFunction::test_apply"); return (t1.get_final_success()) ? 0 : 1; } int test_scalar_multiplication(World& world, std::shared_ptr ncf, const Molecule& molecule, - const CCParameters& parameter) { - int success=0; + const CCParameters& parameter) { CCTimer timer(world, "testing"); test_output t1("CCPairFunction::test_scalar_multiplication"); @@ -464,7 +690,6 @@ int test_scalar_multiplication(World& world, std::shared_ptr::get_thresh(); t1.checkpoint(bsuccess,"scaling"); - if (bsuccess) success++; t1.end(); return (t1.get_final_success()) ? 0 : 1; @@ -472,7 +697,6 @@ int test_scalar_multiplication(World& world, std::shared_ptr ncf, const Molecule& molecule, const CCParameters& parameters) { - int success = 0; test_output t1("CCPairFunction::swap_particles"); CCTimer timer(world, "testing swap_particles"); @@ -535,7 +759,7 @@ int test_swap_particles(World& world, std::shared_ptr print("difference norms in exp(-21 - 21):",norm21_21,ref_12_12); double total_error= fabs(norm12_12-ref_12_12)+ fabs(norm12_21-ref_12_21) - + fabs(norm21_12-ref_12_21)+ fabs(norm21_21-ref_12_12); + + fabs(norm21_12-ref_12_21)+ fabs(norm21_21-ref_12_12); t1.checkpoint(total_error < FunctionDefaults<3>::get_thresh(), "swap_particles u"); @@ -546,8 +770,7 @@ int test_swap_particles(World& world, std::shared_ptr } int test_projector(World& world, std::shared_ptr ncf, const Molecule& molecule, - const CCParameters& parameter) { - int success=0; + const CCParameters& parameter) { test_output t1("CCPairFunction::test_projector"); CCTimer timer(world, "testing"); @@ -561,7 +784,7 @@ int test_projector(World& world, std::shared_ptr ncf, std::vector b = {f3, f1, f2}; std::vector o = orthonormalize_cd({f1-f3, f5}); // projects on an orthonormal basis a = orthonormalize_cd(a); // projects on an orthonormal basis - CCConvolutionOperator f12(world, OT_F12, parameter); + auto f12=CCConvolutionOperatorPtr(world, OT_F12, parameter); { double n1=inner(f1,o[0]); @@ -571,7 +794,7 @@ int test_projector(World& world, std::shared_ptr ncf, } CCPairFunction p1(a,b); - CCPairFunction p2(&f12,a,b); + CCPairFunction p2(f12,a,b); CCPairFunction p3(f); // outer (f1,f2) std::vector vp1({p1}); @@ -660,12 +883,88 @@ int test_projector(World& world, std::shared_ptr ncf, int test_dirac_convolution(World& world, std::shared_ptr ncf, const Molecule& molecule, const CCParameters& parameter) { - int success=0; test_output t1("CCPairFunction::test_dirac_convolution"); return (t1.get_final_success()) ? 0 : 1; } +/// testing = 0.032 mEh +int test_helium(World& world, std::shared_ptr ncf, const Molecule& molecule, + const CCParameters& parameters) { + + test_output t1("CCPairFunction::test_helium"); + CCTimer timer(world, "testing"); + + real_function_3d Vnuc = real_factory_3d(world).f([](const coord_3d& r) {return -2.0/(r.normf()+1.e-8);}); + real_function_3d psi = real_factory_3d(world).f([](const coord_3d& r) {return exp(-r.normf());}); + + auto iterate=[&world](const real_function_3d& potential, real_function_3d& psi, double& eps) { + real_convolution_3d op = BSHOperator3D(world, sqrt(-2*eps), 0.001, 1e-6); + real_function_3d Vpsi = (potential*psi); + Vpsi.scale(-2.0).truncate(); + real_function_3d tmp = op(Vpsi).truncate(); + double norm = tmp.norm2(); + real_function_3d r = tmp-psi; + double rnorm = r.norm2(); + double eps_new = eps - 0.5*inner(Vpsi,r)/(norm*norm); + if (world.rank() == 0) { + print("norm=",norm," eps=",eps," err(psi)=",rnorm," err(eps)=",eps_new-eps); + } + psi = tmp.scale(1.0/norm); + eps = eps_new; + }; + psi.truncate(); + psi.scale(1.0/psi.norm2()); + double eps = -0.6; + real_convolution_3d op = CoulombOperator(world, 0.001, 1e-6); + for (int iter=0; iter<10; iter++) { + real_function_3d rho = square(psi).truncate(); + real_function_3d potential = Vnuc + op(rho).truncate(); + iterate(potential, psi, eps); + } + double kinetic_energy = 0.0; + for (int axis=0; axis<3; axis++) { + real_derivative_3d D = free_space_derivative(world, axis); + real_function_3d dpsi = D(psi); + kinetic_energy += inner(dpsi,dpsi); + } + real_function_3d rho = square(psi); + double two_electron_energy = inner(op(rho),rho); + double nuclear_attraction_energy = 2.0*inner(Vnuc,rho); + double total_energy = kinetic_energy + two_electron_energy + nuclear_attraction_energy; + + t1.checkpoint(fabs(total_energy+2.8616533)<1.e-4,"helium iterations",timer.reset()); + print("ke, total", kinetic_energy, total_energy); + + + CCConvolutionOperator::Parameters param; + auto f=CCConvolutionOperatorPtr(world,OT_F12,param); + auto g=CCConvolutionOperatorPtr(world,OT_G12,param); + + CCPairFunction fij(f,psi,psi); + CCPairFunction gij(g,psi,psi); + CCPairFunction ij({psi},{psi}); + std::vector vfij={fij}; + std::vector vgij={gij}; + std::vector vij={ij}; + + StrongOrthogonalityProjector SO(world); + SO.set_spaces({psi},{psi},{psi},{psi}); + std::vector Qfij=SO(vfij); + std::vector Qgij=SO(vgij); + + double result1=inner(vgij,Qfij); + print(")",result1); + double result2=inner(vfij,Qgij); + print("(",result2); + bool good=fabs(result1-result2)<1.e-5; + good=good and fabs(result1+3.2624783e-02)<1.e-5; + t1.checkpoint(good,"V matrix element",timer.reset()); + + return (t1.get_final_success()) ? 0 : 1; + +} + /** functionality * * - ctor OK @@ -697,7 +996,7 @@ int main(int argc, char **argv) { #ifdef USE_GENTENSOR try { - parser.set_keyval("geometry", "source=library,he"); + parser.set_keyval("geometry", "he"); parser.print_map(); Molecule mol(world, parser); mol.print(); @@ -709,14 +1008,23 @@ int main(int argc, char **argv) { mol, nullptr, std::make_pair("slater", 2.0)); isuccess+=test_constructor(world, ncf, mol, ccparam); - isuccess+=test_overlap(world, ncf, mol, ccparam); + isuccess+=test_operator_apply(world, ncf, mol, ccparam); + isuccess+=test_transformations(world, ncf, mol, ccparam); + isuccess+=test_inner(world, ncf, mol, ccparam); + isuccess+=test_multiply(world, ncf, mol, ccparam); + isuccess+=test_multiply_with_f12(world, ncf, mol, ccparam); isuccess+=test_swap_particles(world, ncf, mol, ccparam); isuccess+=test_scalar_multiplication(world, ncf, mol, ccparam); - isuccess+=test_partial_inner(world, ncf, mol, ccparam); + isuccess+=test_partial_inner_3d(world, ncf, mol, ccparam); + isuccess+=test_partial_inner_6d(world, ncf, mol, ccparam); isuccess+=test_projector(world, ncf, mol, ccparam); + FunctionDefaults<3>::set_cubic_cell(-10,10); + isuccess+=test_helium(world,ncf,mol,ccparam); + data1.clear(); + } catch (std::exception& e) { + madness::print("an error occured"); + madness::print(e.what()); data1.clear(); - } catch (...) { - } #else print("could not run test_ccpairfunction: U need to compile with ENABLE_GENTENSOR=1"); diff --git a/src/madness/chem/test_low_rank_function.cc b/src/madness/chem/test_low_rank_function.cc new file mode 100644 index 00000000000..1ebbd506383 --- /dev/null +++ b/src/madness/chem/test_low_rank_function.cc @@ -0,0 +1,700 @@ +// +// Created by Florian Bischoff on 4/24/23. +// + + +#include +#include + + +#include +#include +#include + + +using namespace madness; + + + +int test_lowrank_function(World& world, LowRankFunctionParameters& parameters) { + test_output t1("CCPairFunction::low rank function"); + t1.set_cout_to_terminal(); + madness::default_random_generator.setstate(int(cpu_time())%4149); + madness::default_random_generator.setstate(int(cpu_time())%4149); + std::string id=unique_fileid(); + + constexpr std::size_t LDIM=3; + constexpr std::size_t NDIM=2*LDIM; + print("eps, k, NDIM, id",FunctionDefaults::get_thresh(),FunctionDefaults::get_k(),NDIM,id); + + parameters.print("grid","end"); + std::string f12type=parameters.f12type(); + std::string transpose=parameters.get("transpose"); + + json j; + std::string jsonfilename="test_low_rank_function."+id+".json"; + j["radius"]=parameters.radius(); + j["f12type"]=parameters.f12type(); + j["gamma"]=parameters.gamma(); + j["volume_element"]=parameters.volume_element(); + j["tol"]=parameters.tol(); + j["hard_shell"]=parameters.hard_shell(); + j["transpose"]=transpose; + j["orthomethod"]=parameters.orthomethod(); + j["gridtype"]=parameters.gridtype(); + j["rhsfunctiontype"]=parameters.rhsfunctiontype(); + j["optimize"]=parameters.optimize(); + std::ofstream of(jsonfilename,std::ios::out); + of< offset; + offset.fill(0.0); + Function phi1=FunctionFactory(world).functor([](const Vector& r) + { return 1.0;}); + Function phi2=FunctionFactory(world).functor([&offset](const Vector& r) + { return exp(-1.0*(r-offset).normf());}); + if (transpose=="slater1") std::swap(phi1,phi2); + { + double n1 = phi1.norm2(); + double n2 = phi2.norm2(); + bool first_one = (fabs(phi1({1, 1, 1}) - 1.0) < 1.e-6); + if (world.rank() == 0) { + if (first_one) print("1(1) phi(2)"); + else print("phi(1) 1(2)"); + print("norms", n1, n2); + } + } + + Function one=FunctionFactory(world) + .functor([](const Vector& r) { return 1.0;}); + + + std::shared_ptr f12; + if (f12type=="slaterf12") f12.reset(SlaterF12OperatorPtr(world,parameters.gamma(),1.e-6,FunctionDefaults::get_thresh())); + else if (f12type=="slater") f12.reset(SlaterOperatorPtr(world,parameters.gamma(),1.e-6,FunctionDefaults::get_thresh())); + else { + MADNESS_EXCEPTION(std::string("unknown f12type"+f12type).c_str(),1); + } + + + auto compute_result = [&world, &one](const auto& lrf) { + real_function_3d result=real_factory_3d(world); + for (int r=0; r(world,reference,"reference."+id,PlotParameters(world).set_plane({"x1","x2"})); + double n2=reference.norm2(); + print("reference.norm2() = int f12 phi2 d2",n2); + output(0.0,0.0,0.0,0.0,0.0,0.0); + + LRFunctorF12 lrfunctor(f12,phi1,phi1); + double cpu0=cpu_time(); + auto lrf=LowRankFunctionFactory(parameters).project(lrfunctor); +// plot_plane<6>(world,lrfunctor,"plot_original."+id,PlotParameters(world).set_plane({"x1","x4"})); + double cpu1=cpu_time(); + double error1=lrf.l2error(lrfunctor); + print("l2error projection",error1); +// plot_plane<6>(world,lrf,"plot_lrf_projection."+id,PlotParameters(world).set_plane({"x1","x4"})); + + // compare + // \phi(1) \bar \phi(1) = \int phi(1) \phi(2) f(1,2) d2 + // = \int \sum_r g_r(1) h_r(2) d2 + // = \sum_r g_r(1) <\phi|h_r> + real_function_3d result=compute_result(lrf); + double projection_error=compute_relative_error(reference,result,lrf); + auto diff=reference-result; +// plot_plane<3>(world,diff,"plot_diff_int_projection."+id,PlotParameters(world).set_plane({"x1","x2"})); +// plot_plane<3>(world,result,"plot_lrf_int_projection."+id,PlotParameters(world).set_plane({"x1","x2"})); + output(projection_error,lrf.rank(),cpu1-cpu0,0.0,0.0,0.0); + j["projection_error_integrated"]=projection_error; + j["projection_error_l2"]=error1; + j["projection_time"]=cpu1-cpu0; + std::ofstream of1(jsonfilename,std::ios::out); + of1<(world,diff,"plot_diff_int_optimization."+id,PlotParameters(world).set_plane({"x1","x2"})); + plot_plane<3>(world,result,"plot_lrf_int_optimization."+id,PlotParameters(world).set_plane({"x1","x2"})); + double optimization_error=compute_relative_error(reference,result,lrf); + output(projection_error,lrf.rank(),cpu1-cpu0,optimization_error,lrf.rank(),cpu3-cpu2); + bool success=(projection_error<5.e-2) and (optimization_error<1.e-2); + + j["optimization_error_integrated"]=optimization_error; + j["optimization_error_l2"]=error2; + j["optimization_time"]=cpu3-cpu2; + + std::ofstream of2(jsonfilename,std::ios::out); + of2< +int test_Kcommutator(World& world, LowRankFunctionParameters& parameters) { + test_output t1("CCPairFunction::low exchange commutator"); + t1.set_cout_to_terminal(); + madness::default_random_generator.setstate(int(cpu_time())%4149); + std::string id=unique_fileid(); + + constexpr std::size_t LDIM=3; + constexpr std::size_t NDIM=2*LDIM; + print("eps, k, NDIM, id",FunctionDefaults::get_thresh(),FunctionDefaults::get_k(),NDIM,id); + parameters.print("grid"); + + real_convolution_3d g12=(CoulombOperator(world,1.e-6,FunctionDefaults::get_thresh())); + std::shared_ptr f12ptr; + std::string f12type=parameters.f12type(); + if (f12type=="slaterf12") f12ptr.reset(SlaterF12OperatorPtr(world,parameters.gamma(),1.e-6,FunctionDefaults::get_thresh())); + else if (f12type=="slater") f12ptr.reset(SlaterOperatorPtr(world,parameters.gamma(),1.e-6,FunctionDefaults::get_thresh())); + else { + MADNESS_EXCEPTION(std::string("unknown f12type"+f12type).c_str(),1); + } + real_convolution_3d& f12=*f12ptr; + + real_function_3d phi=real_factory_3d(world).f([](const coord_3d& r){return exp(-r.normf());}); + double n=phi.norm2(); + phi.scale(1/n); + real_function_3d phi_k=phi; // lookys silly, helps reading. + + + // reference term ( < ij | K(1) ) = = < Ki(1) j(2) | f12 | i(1) j(2) > = ",reference); + + json j; + std::string jsonfilename="test_kcommuntator."+id+".json"; + j["radius"]=parameters.radius(); + j["f12type"]=parameters.f12type(); + j["gamma"]=parameters.gamma(); + j["thresh"]=FunctionDefaults<3>::get_thresh(); + j["volume_element"]=parameters.volume_element(); + j["tol"]=parameters.tol(); + j["hard_shell"]=parameters.hard_shell(); + j["orthomethod"]=parameters.orthomethod(); + j["gridtype"]=parameters.gridtype(); + j["rhsfunctiontype"]=parameters.rhsfunctiontype(); + j["optimize"]=parameters.optimize(); + j["reference"]=reference; + + auto json2file= [](const json& j, const std::string& jsonfilename) { + std::ofstream of(jsonfilename, std::ios::out); + of << j; + of.close(); + }; + + json2file(j,jsonfilename); + timer t(world); + auto compute_error = [&](const std::string msg, const LowRankFunction& lrf) { + auto gk = mul(world, phi_k, g12(lrf.g * phi_k)); // function of 1 + auto hj = lrf.h * phi; // function of 2 + Tensor j_hj = inner(world, phi, hj); + Tensor i_gk = inner(world, phi, gk); + double result_right = j_hj.trace(i_gk); + print(msg, result_right); + j[msg]=result_right-reference; + j[msg+"_rank"]=lrf.rank(); + j[msg+"_compute_time"]=t.tag(msg+"_compute_time"); + json2file(j,jsonfilename); + }; + + + { + // lowrankfunction left phi: lrf(1',2) = f12(1',2) i(1') + // K f12 ij = \sum_k k(1) \int g(1,1') f12(1'2) i(1') j(2) k(1') d1' + // = \sum_kr k(1) j(2) \int g(1,1') g_r(1') h_r(2) k(1') d1' + // = \sum_r j(2) h_r(2) \sum_k k(1) \int g(1,1') g_r(1') k(1') d1' + real_function_3d one = real_factory_3d(world).f([](const coord_3d& r) { return 1.0; }); + LRFunctorF12 lrfunctor(f12ptr,phi,one); +// LowRankFunction fi_one(f12ptr, copy(phi), copy(one)); + auto fi_one=LowRankFunctionFactory(parameters).project(lrfunctor); +// fi_one.project(parameters); + double l2error=fi_one.l2error(lrfunctor); + print("left_project_l2error",l2error); + + j["left_project_time"]=t.tag("left_project_time"); + json2file(j,jsonfilename); + compute_error("left_project",fi_one); + + fi_one.optimize(lrfunctor); + l2error=fi_one.l2error(lrfunctor); + print("left_optimize_l2error",l2error); + j["left_optimize_time"]=t.tag("left_optimize_time"); + json2file(j,jsonfilename); + compute_error("left_optimize",fi_one); + + fi_one.reorthonormalize(); + j["left_reorthonormalize"]=t.tag("left_reorthonormalize"); + json2file(j,jsonfilename); + compute_error("left_reorthonormalize",fi_one); + } + +// // lowrankfunction right phi: lrf(1',2) = f12(1',2) i(1') +// { +// real_function_3d one = real_factory_3d(world).f([](const coord_3d &r) { return 1.0; }); +// LowRankFunction fi_one(f12ptr, copy(one), copy(phi)); +// fi_one.project(parameters); +// std::swap(fi_one.g,fi_one.h); +// j["right_project_time"]=t.tag("right_project_time"); +// json2file(j,jsonfilename); +// +// +// { +// auto gk = mul(world, phi_k, g12(fi_one.g * phi_k)); // function of 1 +// auto hj = fi_one.h * phi; // function of 2 +// Tensor j_hj = inner(world, phi, hj); +// Tensor i_gk = inner(world, phi, gk); +// double result_right = j_hj.trace(i_gk); +// print("result_right, project only", result_right); +// j["right_project"]=result_right-reference; +// j["right_project_rank"]=fi_one.rank(); +// j["left_optimize_compute_time"]=t.tag("left_optimize_compute_time"); +// j["right_project_compute_time"]=t.tag("right_project_compute_time"); +// } +// json2file(j,jsonfilename); +// std::swap(fi_one.g,fi_one.h); +// fi_one.optimize(); +// std::swap(fi_one.g,fi_one.h); +// { +// auto gk = mul(world, phi_k, g12(fi_one.g * phi_k)); // function of 1 +// auto hj = fi_one.h * phi; // function of 2 +// Tensor j_hj = inner(world, phi, hj); +// Tensor i_gk = inner(world, phi, gk); +// double result_right = j_hj.trace(i_gk); +// print("result_right, optimize", result_right); +// j["right_optimize"]=result_right-reference; +// j["right_optimize_rank"]=fi_one.rank(); +// j["right_optimize_compute_time"]=t.tag("right_optimize_compute_time"); +// } +// json2file(j,jsonfilename); +// +// } + + return 0; + +} + +template +int test_full_rank_functor(World& world, LowRankFunctionParameters& parameters) { + + test_output t1("test_full_rank_functor"); +// t1.set_cout_to_terminal(); + print_header2("entering test_full_rank_functor"); + constexpr int NDIM=2*LDIM; + FunctionDefaults::set_thresh(1.e-6); + FunctionDefaults::set_thresh(1.e-6); + double tol=1.e-3; + double gaussexponent=2.0; + + const particle p1=particle::particle1(); + const particle p2=particle::particle2(); + + LRFunctorPure functorpure; + Function gauss=FunctionFactory(world) + .functor([&gaussexponent](const Vector& r){ + Vector a,b; + for (int i=0; i functorf12; + Function b=FunctionFactory(world).functor([](const Vector& r){return exp(-inner(r,r));}); + functorf12.f12.reset(GaussOperatorPtr(world,gaussexponent)); + + auto builder= LowRankFunctionFactory(parameters).set_radius(8) + .set_volume_element(0.1).set_rank_revealing_tol(1.e-10).set_orthomethod("canonical"); + + auto lrfunction1=builder.project(functorf12); + t1.checkpoint(true,"construction f12 functor"); + auto lrfunction2=builder.project(functorpure); + t1.checkpoint(true,"construction full rank functor"); + lrfunction1.optimize(functorf12); + t1.checkpoint(true,"optimization f12 functor"); + lrfunction2.optimize(functorpure); + t1.checkpoint(true,"optimization full rank functor"); + + + try { + double error1=lrfunction1.l2error(functorf12); + t1.checkpoint(error1 +int test_arithmetic(World& world, LowRankFunctionParameters parameters) { + constexpr std::size_t NDIM = 2 * LDIM; + test_output t1("LowRankFunction::arithmetic in dimension " + std::to_string(NDIM)); + t1.set_cout_to_terminal(); + double thresh=FunctionDefaults::get_thresh(); + double thresh_ndim=FunctionDefaults::get_thresh(); + print("thresh ldim/ndim",thresh,thresh_ndim); + Function phi=FunctionFactory(world) + .functor([](const Vector& r){return exp(-4.0*inner(r,r));}); + + LRFunctorF12 functor1; + functor1.f12.reset(GaussOperatorPtr(world,1.0)); + functor1.a=phi; + LRFunctorF12 functor2; + functor2.f12.reset(GaussOperatorPtr(world,2.0)); + functor2.a=phi; + + auto p1=particle::particle1(); + auto p2=particle::particle2(); + + auto builder= LowRankFunctionFactory(parameters).set_radius(4) + .set_volume_element(0.1).set_rank_revealing_tol(1.e-10).set_orthomethod("canonical"); + auto lrf1=builder.project(functor1); + auto lrf2=builder.project(functor2); + + Vector r; + r.fill(0.2); + + // addition/subtraction + { + auto l1=lrf1+lrf1; + t1.checkpoint(fabs(l1(r)-2.0*lrf1(r)) +int test_inner(World& world, LowRankFunctionParameters parameters) { + + static_assert(LDIM==1 or LDIM==2); + constexpr std::size_t NDIM = 2 * LDIM; + test_output t1("LowRankFunction::test_inner in dimension " + std::to_string(NDIM)); + t1.set_cout_to_terminal(); + double thresh=FunctionDefaults::get_thresh(); + double thresh_ndim=FunctionDefaults::get_thresh(); + print("thresh ldim/ndim",thresh,thresh_ndim); + Function phi=FunctionFactory(world) + .functor([](const Vector& r){return exp(-4.0*inner(r,r));}); + + LRFunctorF12 functor1; +// functor1.f12.reset(GaussOperatorPtr(world,1.0)); + functor1.f12.reset(SlaterOperatorPtr_ND(world,1.0,1.e-4,thresh)); + functor1.a=phi; + LRFunctorF12 functor2; +// functor2.f12.reset(GaussOperatorPtr(world,2.0)); + functor2.f12.reset(SlaterOperatorPtr_ND(world,2.0,1.e-4,thresh)); + functor2.a=phi; + + auto p1=particle::particle1(); + auto p2=particle::particle2(); + + auto builder= LowRankFunctionFactory(parameters).set_radius(4) + .set_volume_element(0.1).set_rank_revealing_tol(1.e-10).set_orthomethod("canonical"); + auto lrf1=builder.project(functor1); + auto lrf2=builder.project(functor2); + + // reference numbers: (by mathematica) + // f1(x,y) = exp(-a*x^2) * exp(-(x-y)^2) + // f2(x,y) = exp(-a*x^2) * exp(-g (x-y)^2) + // with a=4, g=2 + // int f1(x,y),f2(x,z) dx = inner(f1,f2,0,0) : norm^2 = Pi^2/(2 Sqrt[2] Sqrt[a gamma] Sqrt[1 + 2 a + gamma]) = 0.37197471167788324677 + // int f1(x,y),f2(z,x) dx = inner(f1,f2,0,1) : norm^2 = 0.32972034117743393239 + // int f1(y,x),f2(x,z) dx = inner(f1,f2,1,0) : norm^2 = 0.26921553123369812300 + // int f1(y,x),f2(z,x) dx = inner(f1,f2,1,1) : norm^2 = 0.35613867236025352322 + + // inner f(1,2) f(2,3) + auto fullrank1=lrf1.reconstruct(); + auto fullrank2=lrf2.reconstruct(); + t1.checkpoint(true,"prep inner"); + { + std::vector reference={ 0.37197471167788324677, 0.32972034117743393239, 0.26921553123369812300, 0.35613867236025352322}; + int counter=0; + for (auto p11 : {p1,p2}) { + for (auto p22 : {p1,p2}) { + double ref=reference[counter]; + if (LDIM==2) ref*=ref; + + // full/full + auto lhs1=inner(fullrank1,fullrank2,p11.get_tuple(),p22.get_tuple()); + double l1=lhs1.norm2(); + t1.checkpoint(fabs(l1*l1-ref)::get_thresh()*50.0; + // fresh start + lrf1=builder.project(functor1); + fullrank1=FunctionFactory(world).functor(functor1); + + std::vector> arg(3); + for (int i=0; i<3; ++i) arg[i]=FunctionFactory(world) + .functor([&i](const Vector& r) + {return exp(-r.normf());}); + + std::vector> lhs_full1, lhs_full2,lhs_func1,lhs_func2; + for (auto& a : arg) { + lhs_full1.push_back(inner(fullrank1,a,p1.get_tuple(),p1.get_tuple())); + lhs_full2.push_back(inner(fullrank1,a,p2.get_tuple(),p1.get_tuple())); + + lhs_func1.push_back(inner(functor1,a,p1,p1)); + lhs_func2.push_back(inner(functor1,a,p2,p1)); + } + auto lhs_lrf1=inner(lrf1,arg,p1,p1); + auto lhs_lrf2=inner(lrf1,arg,p2,p1); + + double norm_func1=norm2(world,lhs_func1); + double norm_func2=norm2(world,lhs_func2); + double norm_full1=norm2(world,lhs_full1); + double norm_full2=norm2(world,lhs_full2); + double norm_lrf1=norm2(world,lhs_lrf1); + double norm_lrf2=norm2(world,lhs_lrf2); + print("norms 1",norm_func1,norm_full1,norm_lrf1); + print("norms 2",norm_func2,norm_full2,norm_lrf2); + + double error1=norm2(world,lhs_full1-lhs_func1); + double error2=norm2(world,lhs_full2-lhs_func2); + double error3=norm2(world,lhs_lrf1 -lhs_func1); + double error4=norm2(world,lhs_lrf2 -lhs_func2); + + print("error1/2",error1,error2,error3,error4); +// t1.checkpoint(error1 +int test_construction_optimization(World& world, LowRankFunctionParameters parameters) { + parameters.set_user_defined_value("volume_element",0.05); + constexpr std::size_t NDIM=2*LDIM; + test_output t1("LowRankFunction::construction/optimization in dimension "+std::to_string(NDIM)); + t1.set_cout_to_terminal(); + OperatorInfo info(1.0,1.e-6,FunctionDefaults::get_thresh(),OT_SLATER); + auto slater=std::shared_ptr >(new SeparatedConvolution(world,info)); + Function one=FunctionFactory(world).functor([](const Vector& r){return exp(-0.2*inner(r,r));}); + + LRFunctorF12 lrfunctor(slater,one,one); + LowRankFunctionFactory builder(parameters); + auto lrf=builder.project(lrfunctor); + t1.checkpoint(lrf.rank()>0,"construction"); + + double error=lrf.l2error(lrfunctor); + print("l2 error project ",error); + t1.checkpoint(error<2.e-2,"l2 error in projection "+std::to_string(error)); + + auto lrf2(lrf); + error=lrf2.l2error(lrfunctor); + print("l2 error copy ctor ",error); + MADNESS_CHECK(lrf.rank()==lrf2.rank()); + MADNESS_CHECK(&(lrf.g[0]) != &(lrf2.g[0])); // deep copy + t1.checkpoint(error<2.e-2,"l2 error in copy ctor "+std::to_string(error)); + + lrf.optimize(lrfunctor); + error=lrf.l2error(lrfunctor); + print("l2 error optimize",error); + t1.checkpoint(error<1.e-2,"l2 error in optimization "+std::to_string(error)); + + lrf.reorthonormalize(); + error=lrf.l2error(lrfunctor); + print("l2 error reorthonormalize",error); + t1.checkpoint(error<1.e-2,"l2 error in reorthonormalization "+std::to_string(error)); + return t1.end(); +} + +int main(int argc, char **argv) { + + madness::World& world = madness::initialize(argc, argv); + startup(world, argc, argv); + commandlineparser parser(argc, argv); + int k = parser.key_exists("k") ? std::atoi(parser.value("k").c_str()) : 6; + double thresh = parser.key_exists("thresh") ? std::stod(parser.value("thresh")) : 1.e-4; + FunctionDefaults<6>::set_tensor_type(TT_2D); + + FunctionDefaults<1>::set_thresh(thresh); + FunctionDefaults<2>::set_thresh(thresh); + FunctionDefaults<3>::set_thresh(thresh); + FunctionDefaults<4>::set_thresh(thresh); + FunctionDefaults<5>::set_thresh(thresh); + FunctionDefaults<6>::set_thresh(thresh); + + FunctionDefaults<1>::set_k(k); + FunctionDefaults<2>::set_k(k); + FunctionDefaults<3>::set_k(k); + FunctionDefaults<4>::set_k(k); + FunctionDefaults<5>::set_k(k); + FunctionDefaults<6>::set_k(k); + + FunctionDefaults<1>::set_cubic_cell(-10.,10.); + FunctionDefaults<2>::set_cubic_cell(-10.,10.); + FunctionDefaults<3>::set_cubic_cell(-10.,10.); + FunctionDefaults<4>::set_cubic_cell(-10.,10.); + FunctionDefaults<5>::set_cubic_cell(-10.,10.); + FunctionDefaults<6>::set_cubic_cell(-10.,10.); + + FunctionDefaults<2>::set_tensor_type(TT_FULL); + print("numerical parameters: k, eps(3D), eps(6D)", FunctionDefaults<3>::get_k(), FunctionDefaults<3>::get_thresh(), + FunctionDefaults<6>::get_thresh()); + LowRankFunctionParameters parameters; + parameters.read_and_set_derived_values(world,parser,"grid"); + parameters.print("grid"); + int isuccess=0; +#ifdef USE_GENTENSOR + + + try { + +// isuccess+=test_grids<1>(world,parameters); +// isuccess+=test_grids<2>(world,parameters); +// isuccess+=test_grids<3>(world,parameters); +// isuccess+=test_full_rank_functor<1>(world, parameters); +// isuccess+=test_construction_optimization<1>(world,parameters); +// isuccess+=test_construction_optimization<2>(world,parameters); + isuccess+=test_arithmetic<1>(world,parameters); + isuccess+=test_arithmetic<2>(world,parameters); + isuccess+=test_inner<1>(world,parameters); + isuccess+=test_inner<2>(world,parameters); + +// isuccess+=test_lowrank_function(world,parameters); +// isuccess+=test_Kcommutator(world,parameters); + } catch (std::exception& e) { + madness::print("an error occured"); + madness::print(e.what()); + } +#else + print("could not run test_ccpairfunction: U need to compile with ENABLE_GENTENSOR=1"); +#endif + finalize(); + + return isuccess; +} \ No newline at end of file diff --git a/src/madness/misc/CMakeLists.txt b/src/madness/misc/CMakeLists.txt index 81f5842b3d3..790a0d1e09a 100644 --- a/src/madness/misc/CMakeLists.txt +++ b/src/madness/misc/CMakeLists.txt @@ -2,7 +2,7 @@ set(MADMISC_HEADERS misc.h ran.h phandler.h interpolation_1d.h cfft.h info.h) set(MADMISC_SOURCES - checksum_file.cc position_stream.cc gprofexit.cc ran.cc cfft.cc info.cc) + checksum_file.cc position_stream.cc gprofexit.cc ran.cc cfft.cc info.cc unique_filename.cc) # retrieve git metadata include(GetGitMetadata) vgkit_cmake_git_metadata() diff --git a/src/madness/misc/misc.h b/src/madness/misc/misc.h index 313c53dcdeb..f8ace2438f1 100644 --- a/src/madness/misc/misc.h +++ b/src/madness/misc/misc.h @@ -46,6 +46,8 @@ namespace madness { const char comment='#', bool rewind=true, bool silent=false); std::string lowercase(const std::string& s); void gprofexit(int id, int nproc); + /// creates a unique filename, using PBS ID if available + std::string unique_fileid(); } #endif // MADNESS_MISC_MISC_H__INCLUDED diff --git a/src/madness/misc/unique_filename.cc b/src/madness/misc/unique_filename.cc new file mode 100644 index 00000000000..7788ef328e4 --- /dev/null +++ b/src/madness/misc/unique_filename.cc @@ -0,0 +1,25 @@ +// +// Created by Florian Bischoff/chatgpt on 9/7/23. +// + +#include +#include + +namespace madness { + +std::string unique_fileid() { + std::string uniqueFileName; + + // Check if the PBS_JOBID environment variable is set + const char* pbsId = std::getenv("PBS_JOBID"); + if (pbsId != nullptr) { + uniqueFileName =std::string(pbsId); + } else { + // If PBS_ID is not available, generate a random number + int randomNumber = RandomValue(); + uniqueFileName = std::to_string(randomNumber); + } + + return uniqueFileName; +} +} diff --git a/src/madness/mra/funcimpl.h b/src/madness/mra/funcimpl.h index 2c89e977652..17ce65612cf 100644 --- a/src/madness/mra/funcimpl.h +++ b/src/madness/mra/funcimpl.h @@ -519,7 +519,7 @@ namespace madness { std::vector s0(NDIM/2,s); const double tol=f->get_thresh(); - const double thresh=f->truncate_tol(tol, key); + const double thresh=f->truncate_tol(tol, key)*0.3; // custom factor to "ensure" accuracy // include the wavelets in the norm, makes it much more accurate const double fnorm=fcoeff.normf(); const double gnorm=gcoeff.normf(); @@ -913,6 +913,7 @@ namespace madness { private: typedef WorldObject< FunctionImpl > woT; ///< Base class world object type public: + typedef T typeT; typedef FunctionImpl implT; ///< Type of this class (implementation) typedef std::shared_ptr< FunctionImpl > pimplT; ///< pointer to this class typedef Tensor tensorT; ///< Type of tensor for anything but to hold coeffs @@ -3263,11 +3264,11 @@ namespace madness { /// /// k=number of wavelets, so k=5 means max order is 4, so max exactly /// representable squarable polynomial is of order 2. - void tnorm(const tensorT& t, double* lo, double* hi) const; + void static tnorm(const tensorT& t, double* lo, double* hi); - void tnorm(const GenTensor& t, double* lo, double* hi) const; + void static tnorm(const GenTensor& t, double* lo, double* hi); - void tnorm(const SVDTensor& t, double* lo, double* hi, const int particle) const; + void static tnorm(const SVDTensor& t, double* lo, double* hi, const int particle); // This invoked if node has not been autorefined void do_square_inplace(const keyT& key); @@ -3764,69 +3765,63 @@ namespace madness { + /// pointwise multiplication of two tensors, returns result and estimates error + + /// provide one of the two factors upon construction, the other factor upon operator() call. + /// + /// error is estimated by oversampling: pointwise multiplication will result in a coefficient + /// tensor of order 2k, estimate the error through the norm of the k+1 contribution + /// + /// error does not account for inaccurate representation of the input tensors! + /// U need to compute that somewhere else! template struct pointwise_multiplier { - const implT* impl; - const FunctionImpl* gimpl; coeffT val_lhs, coeff_lhs; + long oversampling=1; double error=0.0; double lo=0.0, hi=0.0, lo1=0.0, hi1=0.0, lo2=0.0, hi2=0.0; - pointwise_multiplier() :gimpl(0), impl(0) {} - pointwise_multiplier(const Key key, const coeffT& clhs, implT* i, const FunctionImpl* gimpl) - : impl(i), gimpl(gimpl), coeff_lhs(clhs) { - val_lhs=impl->coeffs2values(key,coeff_lhs); + pointwise_multiplier() {} + pointwise_multiplier(const Key key, const coeffT& clhs) : coeff_lhs(clhs) { + auto fcf=FunctionCommonFunctionality(coeff_lhs.dim(0)); + val_lhs=fcf.coeffs2values(key,coeff_lhs); error=0.0; - impl->tnorm(coeff_lhs,&lo,&hi); - gimpl->tnorm(coeff_lhs.get_svdtensor(),&lo1,&hi1,1); - gimpl->tnorm(coeff_lhs.get_svdtensor(),&lo2,&hi2,2); - } - - pointwise_multiplier(const Key key, const coeffT& clhs, implT* i) - : impl(i), gimpl(0), coeff_lhs(clhs) { - val_lhs=impl->coeffs2values(key,coeff_lhs); - error=0.0; - impl->tnorm(coeff_lhs,&lo,&hi); - } - - /// multiply values of rhs and lhs, result on rhs, rhs and lhs are of the same dimensions - coeffT operator()(const Key key, const coeffT& coeff_rhs) { - double rlo, rhi; - impl->tnorm(coeff_rhs,&rlo,&rhi); - error = hi*rlo + rhi*lo + rhi*hi; - coeffT val_rhs=impl->coeffs2values(key, coeff_rhs); - val_rhs.emul(val_lhs); - return impl->values2coeffs(key,val_rhs); - + implT::tnorm(coeff_lhs,&lo,&hi); + if (coeff_lhs.is_svd_tensor()) { + FunctionImpl::tnorm(coeff_lhs.get_svdtensor(),&lo1,&hi1,1); + FunctionImpl::tnorm(coeff_lhs.get_svdtensor(),&lo2,&hi2,2); + } } /// multiply values of rhs and lhs, result on rhs, rhs and lhs are of the same dimensions tensorT operator()(const Key key, const tensorT& coeff_rhs) { MADNESS_ASSERT(coeff_rhs.dim(0)==coeff_lhs.dim(0)); + auto fcf=FunctionCommonFunctionality(coeff_lhs.dim(0)); // the tnorm estimate is not tight enough to be efficient, better use oversampling bool use_tnorm=false; if (use_tnorm) { double rlo, rhi; - impl->tnorm(coeff_rhs,&rlo,&rhi); + implT::tnorm(coeff_rhs,&rlo,&rhi); error = hi*rlo + rhi*lo + rhi*hi; - tensorT val_rhs=impl->coeffs2values(key, coeff_rhs); + tensorT val_rhs=fcf.coeffs2values(key, coeff_rhs); val_rhs.emul(val_lhs.full_tensor_copy()); - return impl->values2coeffs(key,val_rhs); + return fcf.values2coeffs(key,val_rhs); } else { // use quadrature of order k+1 - auto cdata=FunctionCommonData::get(impl->get_k()+1); // npt=k+1 - FunctionCommonFunctionality fcf_hi_npt(cdata); + auto& cdata=FunctionCommonData::get(coeff_rhs.dim(0)); // npt=k+1 + auto& cdata_npt=FunctionCommonData::get(coeff_rhs.dim(0)+oversampling); // npt=k+1 + FunctionCommonFunctionality fcf_hi_npt(cdata_npt); // coeffs2values for rhs: k -> npt=k+1 - tensorT coeff1(cdata.vk); - coeff1(impl->cdata.s0)=coeff_rhs; // s0 is smaller than vk! + tensorT coeff1(cdata_npt.vk); + coeff1(cdata.s0)=coeff_rhs; // s0 is smaller than vk! tensorT val_rhs_k1=fcf_hi_npt.coeffs2values(key,coeff1); // coeffs2values for lhs: k -> npt=k+1 - tensorT coeff_lhs_k1(cdata.vk); - coeff_lhs_k1(impl->cdata.s0)=coeff_lhs.full_tensor_copy(); + tensorT coeff_lhs_k1(cdata_npt.vk); + coeff_lhs_k1(cdata.s0)=coeff_lhs.full_tensor_copy(); tensorT val_lhs_k1=fcf_hi_npt.coeffs2values(key,coeff_lhs_k1); // multiply @@ -3836,24 +3831,25 @@ namespace madness { tensorT result1=fcf_hi_npt.values2coeffs(key,val_lhs_k1); // extract coeffs up to k - tensorT result=copy(result1(impl->cdata.s0)); - result1(impl->cdata.s0)=0.0; + tensorT result=copy(result1(cdata.s0)); + result1(cdata.s0)=0.0; error=result1.normf(); return result; - } - } /// multiply values of rhs and lhs, result on rhs, rhs and lhs are of differnet dimensions coeffT operator()(const Key key, const tensorT& coeff_rhs, const int particle) { Key key1, key2; key.break_apart(key1,key2); - MADNESS_ASSERT(gimpl); - FunctionCommonFunctionality fcf_lo(gimpl->cdata); - FunctionCommonFunctionality fcf_hi(impl->cdata); - FunctionCommonFunctionality fcf_lo_npt(gimpl->get_k()+1); - FunctionCommonFunctionality fcf_hi_npt(impl->get_k()+1); + const long k=coeff_rhs.dim(0); + auto& cdata=FunctionCommonData::get(k); + auto& cdata_lowdim=FunctionCommonData::get(k); + FunctionCommonFunctionality fcf_lo(cdata_lowdim); + FunctionCommonFunctionality fcf_hi(cdata); + FunctionCommonFunctionality fcf_lo_npt(k+oversampling); + FunctionCommonFunctionality fcf_hi_npt(k+oversampling); + // make hi-dim values from lo-dim coeff_rhs on npt grid points tensorT ones=tensorT(fcf_lo_npt.cdata.vk); @@ -3880,16 +3876,14 @@ namespace madness { coeffT result1=fcf_hi_npt.values2coeffs(key,val_lhs_npt); // extract coeffs up to k - coeffT result=copy(result1(impl->cdata.s0)); - result1(impl->cdata.s0)=0.0; + coeffT result=copy(result1(cdata.s0)); + result1(cdata.s0)=0.0; error=result1.normf(); - result.reduce_rank(impl->get_tensor_args().thresh); return result; - } template void serialize(const Archive& ar) { - ar & error & lo & lo1 & lo2 & hi & hi1& hi2 & gimpl & impl & val_lhs & coeff_lhs; + ar & error & lo & lo1 & lo2 & hi & hi1& hi2 & val_lhs & coeff_lhs; } @@ -4109,20 +4103,14 @@ namespace madness { double error=refine_error; // prepare the multiplication - pointwise_multiplier pm; - if (have_v1()) pm=pointwise_multiplier(key,coeff_ket,result,iav1.get_impl()); - else if (have_v2()) { - pm=pointwise_multiplier(key,coeff_ket,result,iav2.get_impl()); - } else { - pm=pointwise_multiplier(key,coeff_ket,result); - } + pointwise_multiplier pm(key,coeff_ket); // perform the multiplication, compute tnorm part of the total error coeffT cresult(result->cdata.vk,result->get_tensor_args()); if (have_v1()) { cresult+=pm(key,cpot1.get_tensor(),1); error+=pm.error; - } + } if (have_v2()) { cresult+=pm(key,cpot2.get_tensor(),2); error+=pm.error; @@ -5342,6 +5330,129 @@ namespace madness { (rangeT(coeffs.begin(),coeffs.end()),do_inner_local(&g, leaves_only)); } + + /// compute the inner product of this range with other + template + struct do_inner_local_on_demand { + const FunctionImpl* ket; + const FunctionImpl* bra; + bool leaves_only=true; + typedef TENSOR_RESULT_TYPE(T,R) resultT; + + do_inner_local_on_demand(const FunctionImpl* bra, const FunctionImpl* ket, + const bool leaves_only=true) + : bra(bra), ket(ket), leaves_only(leaves_only) {} + resultT operator()(typename dcT::const_iterator& it) const { + + constexpr std::size_t LDIM=std::max(NDIM/2,std::size_t(1)); + + const keyT& key=it->first; + const nodeT& fnode = it->second; + if (not fnode.has_coeff()) return resultT(0.0); // probably internal nodes + + // assuming all boxes (esp the low-dim ones) are local, i.e. the functions are replicated + auto find_valid_parent = [](auto& key, auto& impl, auto&& find_valid_parent) { + MADNESS_CHECK(impl->get_coeffs().owner(key)==impl->world.rank()); // make sure everything is local! + if (impl->get_coeffs().probe(key)) return key; + auto parentkey=key.parent(); + return find_valid_parent(parentkey, impl, find_valid_parent); + }; + + // returns coefficients, empty if no functor present + auto get_coeff = [&find_valid_parent](const auto& key, const auto& impl) { + bool have_impl=impl.get(); + if (have_impl) { + auto parentkey = find_valid_parent(key, impl, find_valid_parent); + MADNESS_CHECK(impl->get_coeffs().probe(parentkey)); + typename decltype(impl->coeffs)::accessor acc; + impl->get_coeffs().find(acc,parentkey); + auto parentcoeff=acc->second.coeff(); + auto coeff=impl->parent_to_child(parentcoeff, parentkey, key); + return coeff; + } else { + return GenTensor::typeT>(); + } + }; + + Key key1,key2; + key.break_apart(key1,key2); + + auto func=dynamic_cast* >(ket->functor.get()); + MADNESS_ASSERT(func); + + auto coeff_bra=fnode.coeff(); + auto coeff_ket=get_coeff(key,func->impl_ket); + auto coeff_v1=get_coeff(key1,func->impl_m1); + auto coeff_v2=get_coeff(key2,func->impl_m2); + auto coeff_p1=get_coeff(key1,func->impl_p1); + auto coeff_p2=get_coeff(key2,func->impl_p2); + + // construct |ket(1,2)> or |p(1)p(2)> or |p(1)p(2) ket(1,2)> + double error=0.0; + if (coeff_ket.has_data() and coeff_p1.has_data()) { + pointwise_multiplier pm(key,coeff_ket); + coeff_ket=pm(key,outer(coeff_p1,coeff_p2,TensorArgs(TT_FULL,-1.0)).full_tensor()); + error+=pm.error; + } else if (coeff_ket.has_data() or coeff_p1.has_data()) { + coeff_ket = (coeff_ket.has_data()) ? coeff_ket : outer(coeff_p1,coeff_p2); + } else { // not ket and no p1p2 + MADNESS_EXCEPTION("confused ket/p1p2 in do_inner_local_on_demand",1); + } + + // construct (v(1) + v(2)) |ket(1,2)> + coeffT v1v2ket; + if (coeff_v1.has_data()) { + pointwise_multiplier pm(key,coeff_ket); + v1v2ket = pm(key,coeff_v1.full_tensor(), 1); + error+=pm.error; + v1v2ket+= pm(key,coeff_v2.full_tensor(), 2); + error+=pm.error; + } else { + v1v2ket = coeff_ket; + } + + resultT result; + if (func->impl_eri) { // project bra*ket onto eri, avoid multiplication with eri + MADNESS_CHECK(func->impl_eri->get_functor()->provides_coeff()); + coeffT coeff_eri=func->impl_eri->get_functor()->coeff(key).full_tensor(); + pointwise_multiplier pm(key,v1v2ket); + tensorT braket=pm(key,coeff_bra.full_tensor_copy().conj()); + error+=pm.error; + if (error>1.e-3) print("error in key",key,error); + result=coeff_eri.full_tensor().trace(braket); + + } else { // no eri, project ket onto bra + result=coeff_bra.full_tensor_copy().trace_conj(v1v2ket.full_tensor_copy()); + } + return result; + } + + resultT operator()(resultT a, resultT b) const { + return (a+b); + } + + template void serialize(const Archive& ar) { + throw "NOT IMPLEMENTED"; + } + }; + + /// Returns the inner product of this with function g constructed on-the-fly + + /// the leaf boxes of this' MRA tree defines the inner product + template + TENSOR_RESULT_TYPE(T,R) inner_local_on_demand(const FunctionImpl& gimpl) const { + PROFILE_MEMBER_FUNC(FunctionImpl); + typedef TENSOR_RESULT_TYPE(T,R) resultT; + typedef FunctionImpl implR; + + MADNESS_CHECK(this->is_reconstructed()); + + typedef Range rangeT; + rangeT range(coeffs.begin(), coeffs.end()); + return world.taskq.reduce>(range, + do_inner_local_on_demand(this, &gimpl)); + } + /// Type of the entry in the map returned by make_key_vec_map typedef std::vector< std::pair > mapvecT; @@ -5904,9 +6015,9 @@ namespace madness { int gparticle= v1[0]==0 ? 0 : 1; // which particle to integrate over int hparticle= v2[0]==0 ? 0 : 1; // which particle to integrate over // merge multiple contraction dimensions into one - Tensor gtensor = gcoeff1.get_svdtensor().make_vector_with_weights(gparticle); + Tensor gtensor = gcoeff1.get_svdtensor().flat_vector_with_weights(gparticle); Tensor gtensor_other = gcoeff1.get_svdtensor().flat_vector((gparticle+1)%2); - Tensor htensor = hcoeff1.get_svdtensor().make_vector_with_weights(hparticle); + Tensor htensor = hcoeff1.get_svdtensor().flat_vector_with_weights(hparticle); Tensor htensor_other = hcoeff1.get_svdtensor().flat_vector((hparticle+1)%2); Tensor tmp1=inner(gtensor,htensor,1,1); // tmp1(r,r') = sum_j b(r,j) a(r',j) Tensor tmp2=inner(tmp1,gtensor_other,0,0); // tmp2(r',i) = sum_r tmp1(r,r') a(r,i) @@ -5918,9 +6029,9 @@ namespace madness { if (key.level() > 0) { GenTensor gcoeff2 = copy(gcoeff1(g->get_cdata().s0)); GenTensor hcoeff2 = copy(hcoeff1(h->get_cdata().s0)); - Tensor gtensor = gcoeff2.get_svdtensor().make_vector_with_weights(gparticle); + Tensor gtensor = gcoeff2.get_svdtensor().flat_vector_with_weights(gparticle); Tensor gtensor_other = gcoeff2.get_svdtensor().flat_vector((gparticle+1)%2); - Tensor htensor = hcoeff2.get_svdtensor().make_vector_with_weights(hparticle); + Tensor htensor = hcoeff2.get_svdtensor().flat_vector_with_weights(hparticle); Tensor htensor_other = hcoeff2.get_svdtensor().flat_vector((hparticle+1)%2); Tensor tmp1=inner(gtensor,htensor,1,1); // tmp1(r,r') = sum_j b(r,j) a(r',j) Tensor tmp2=inner(tmp1,gtensor_other,0,0); // tmp2(r',i) = sum_r tmp1(r,r') a(r,i) diff --git a/src/madness/mra/funcplot.h b/src/madness/mra/funcplot.h index f73a41f1428..995d9e9b692 100644 --- a/src/madness/mra/funcplot.h +++ b/src/madness/mra/funcplot.h @@ -55,7 +55,7 @@ namespace madness { // initialize with: key, value, comment (optional), allowed values (optional) initialize("zoom",2,"zoom into the simulation cell"); initialize("npoints",151,"number of plot points per dimension"); - initialize>("origin",{0.0,0.0,0.0},"origin of the plot"); + initialize>("origin",{},"origin of the plot"); initialize>("plane",{"x1","x2"},"plot plane: x1, x2, .., x6"); } @@ -83,6 +83,9 @@ namespace madness { template Vector origin() const { auto origin_vec=get>("origin"); + // fill in zeros if the default origin has fewer dimensions than the actual origin + int missing=NDIM-origin_vec.size(); + for (auto i=0; i o(origin_vec); return o; } @@ -1021,9 +1024,9 @@ void plot_plane(World& world, const std::vector >& vfuncti World& world=vf.front().world(); for(size_t i=0;i(world,vf[i],namei); - plot_cubefile(world,vf[i],namei+".cube",header); +// plot_cubefile(world,vf[i],namei+".cube",header); } } diff --git a/src/madness/mra/function_interface.h b/src/madness/mra/function_interface.h index b253cdf4283..b0c498730c1 100644 --- a/src/madness/mra/function_interface.h +++ b/src/madness/mra/function_interface.h @@ -189,6 +189,16 @@ namespace madness { } + /// replicate low-dimensional functions over all ranks of this world + void replicate_low_dim_functions(const bool fence) { + if (impl_m1 and (not impl_m1->is_on_demand())) impl_m1->replicate(false); + if (impl_m2 and (not impl_m2->is_on_demand())) impl_m2->replicate(false); + + if (impl_p1 and (not impl_p1->is_on_demand())) impl_p1->replicate(false); + if (impl_p2 and (not impl_p2->is_on_demand())) impl_p2->replicate(false); + if (fence) world.gop.fence(); + } + void make_redundant(const bool fence) { // prepare base functions that make this function if (impl_ket and (not impl_ket->is_on_demand())) impl_ket->make_redundant(false); diff --git a/src/madness/mra/gfit.h b/src/madness/mra/gfit.h index bb13810e930..33a89ba81f1 100644 --- a/src/madness/mra/gfit.h +++ b/src/madness/mra/gfit.h @@ -57,7 +57,17 @@ class GFit { public: /// default ctor does nothing - GFit() {} + GFit() = default; + + /// copy constructor + GFit(const GFit& other) = default; + + /// assignment operator + GFit& operator=(const GFit& other) { + coeffs_ = other.coeffs_; + exponents_ = other.exponents_; + return *this; + } /// return a fit for the Coulomb function @@ -79,8 +89,15 @@ class GFit { /// @parma[in] prnt print level static GFit BSHFit(double mu, double lo, double hi, double eps, bool prnt=false) { GFit fit; - if (NDIM==3) bsh_fit(mu,lo,hi,eps,fit.coeffs_,fit.exponents_,prnt); + bool fix_interval=false; + if (NDIM==3) bsh_fit(mu,lo,hi,eps,fit.coeffs_,fit.exponents_,prnt,fix_interval); else bsh_fit_ndim(NDIM,mu,lo,hi,eps,fit.coeffs_,fit.exponents_,prnt); + + if (prnt) { + print("bsh fit"); + auto exact = [&mu](const double r) -> double { return 1.0/(4.0 * constants::pi) * exp(-mu * r)/r; }; + fit.print_accuracy(exact, lo, hi); + } return fit; } @@ -95,10 +112,161 @@ class GFit { /// @parma[in] prnt print level static GFit SlaterFit(double gamma, double lo, double hi, double eps, bool prnt=false) { GFit fit; - slater_fit(gamma,lo,hi,eps,fit.coeffs_,fit.exponents_,prnt); + slater_fit(gamma,lo,hi,eps,fit.coeffs_,fit.exponents_,false); + if (prnt) { + print("Slater fit"); + auto exact = [&gamma](const double r) -> double { return exp(-gamma * r); }; + fit.print_accuracy(exact, lo, hi); + } return fit; } + /// return a (trivial) fit for a single Gauss function + + /// the Gauss function is defined by + /// f(r) = exp(-\gamma r^2) + /// @param[in] gamma the exponent of the Gauss function + /// @param[in] lo the smallest length scale that needs to be precisely represented + /// @param[in] hi the largest length scale that needs to be precisely represented + /// @param[in] eps the precision threshold + /// @parma[in] prnt print level + static GFit GaussFit(double gamma, double lo, double hi, double eps, bool prnt=false) { + GFit fit; + fit.coeffs_=Tensor(1); + fit.exponents_=Tensor(1); + fit.coeffs_=1.0; + fit.exponents_=gamma; + if (prnt) { + print("Gauss fit"); + auto exact = [&gamma](const double r) -> double { return exp(-gamma * r*r); }; + fit.print_accuracy(exact, lo, hi); + } + return fit; + } + + /// return a fit for the F12 correlation factor + + /// the Slater function is defined by + /// f(r) = 1/(2 gamma) * (1 - exp(-\gamma r)) + /// @param[in] gamma the exponent of the Slater function + /// @param[in] lo the smallest length scale that needs to be precisely represented + /// @param[in] hi the largest length scale that needs to be precisely represented + /// @param[in] eps the precision threshold + /// @parma[in] prnt print level + static GFit F12Fit(double gamma, double lo, double hi, double eps, bool prnt=false) { + GFit fit; + f12_fit(gamma,lo*0.1,hi,eps*0.01,fit.coeffs_,fit.exponents_,false); + fit.coeffs_*=(0.5/gamma); + if (prnt) { + print("f12 fit"); + auto exact=[&gamma](const double r) -> double {return 0.5/gamma*(1.0-exp(-gamma*r));}; + fit.print_accuracy(exact,lo,hi); + } + return fit; + } + + /// return a fit for the F12^2 correlation factor + + /// the Slater function square is defined by + /// f(r) = [ 1/(2 gamma) * (1 - exp(-\gamma r)) ] ^2 + /// @param[in] gamma the exponent of the Slater function + /// @param[in] lo the smallest length scale that needs to be precisely represented + /// @param[in] hi the largest length scale that needs to be precisely represented + /// @param[in] eps the precision threshold + /// @parma[in] prnt print level + static GFit F12sqFit(double gamma, double lo, double hi, double eps, bool prnt=false) { + GFit fit; + f12sq_fit(gamma,lo*0.1,hi,eps*0.01,fit.coeffs_,fit.exponents_,false); + fit.coeffs_*=(0.25/(gamma*gamma)); + if (prnt) { + print("f12sq fit"); + auto exact=[&gamma](const double r) -> double {return std::pow(0.5/gamma*(1.0-exp(-gamma*r)),2.0);}; + fit.print_accuracy(exact,lo,hi); + } + return fit; + } + + + /// return a fit for the FG function + + /// fg = 1/(2 mu) * (1 - exp(-gamma r12)) / r12 + /// = 1/(2 mu) *( 1/r12 - exp(-gamma r12)/r12) + /// = 1/(2 mu) * (coulomb - bsh) + /// @param[in] gamma the exponent of the Slater function + /// @param[in] lo the smallest length scale that needs to be precisely represented + /// @param[in] hi the largest length scale that needs to be precisely represented + /// @param[in] eps the precision threshold + /// @parma[in] prnt print level + static GFit FGFit(double gamma, double lo, double hi, double eps, bool prnt=false) { + GFit bshfit,coulombfit; + eps*=0.1; + lo*=0.1; +// bool restrict_interval=false; + bool fix_interval=true; + bsh_fit(gamma,lo,hi,eps,bshfit.coeffs_,bshfit.exponents_,false,fix_interval); + bsh_fit(0.0,lo,hi,eps,coulombfit.coeffs_,coulombfit.exponents_,false,fix_interval); + // check the exponents are identical + auto diffexponents=(coulombfit.exponents() - bshfit.exponents()); + MADNESS_CHECK(diffexponents.normf()/coulombfit.exponents().size()<1.e-12); + auto diffcoefficients=(coulombfit.coeffs() - bshfit.coeffs()); + GFit fgfit; + fgfit.exponents_=bshfit.exponents_; + fgfit.coeffs_=4.0*constants::pi*0.5/gamma*diffcoefficients; + GFit::prune_small_coefficients(eps,lo,hi,fgfit.coeffs_,fgfit.exponents_); + + if (prnt) { + print("fg fit"); + auto exact=[&gamma](const double r) -> double {return 0.5/gamma*(1.0-exp(-gamma*r))/r;}; + fgfit.print_accuracy(exact,lo,hi); + } + return fgfit; + } + + /// return a fit for the F2G function + + /// f2g = (1/(2 mu) * (1 - exp(-gamma r12)))^2 / r12 + /// = 1/(4 mu^2) * [ 1/r12 - 2 exp(-gamma r12)/r12) + exp(-2 gamma r12)/r12 ] + /// @param[in] gamma the exponent of the Slater function + /// @param[in] lo the smallest length scale that needs to be precisely represented + /// @param[in] hi the largest length scale that needs to be precisely represented + /// @param[in] eps the precision threshold + /// @parma[in] prnt print level + static GFit F2GFit(double gamma, double lo, double hi, double eps, bool prnt=false) { + GFit bshfit,coulombfit,bsh2fit; + eps*=0.1; + lo*=0.1; +// bool restrict_interval=false; + bool fix_interval=true; + bsh_fit(gamma,lo,hi,eps,bshfit.coeffs_,bshfit.exponents_,false,fix_interval); + bsh_fit(2.0*gamma,lo,hi,eps,bsh2fit.coeffs_,bsh2fit.exponents_,false,fix_interval); + bsh_fit(0.0,lo,hi,eps,coulombfit.coeffs_,coulombfit.exponents_,false,fix_interval); + + // check the exponents are identical + auto diffexponents=(coulombfit.exponents() - bshfit.exponents()); + MADNESS_CHECK(diffexponents.normf()/coulombfit.exponents().size()<1.e-12); + auto diffexponents1=(coulombfit.exponents() - bsh2fit.exponents()); + MADNESS_CHECK(diffexponents1.normf()/coulombfit.exponents().size()<1.e-12); + + auto coefficients=(coulombfit.coeffs() - 2.0* bshfit.coeffs() + bsh2fit.coeffs()); + GFit f2gfit; + f2gfit.exponents_=bshfit.exponents_; + // additional factor 4 pi due to implementation of bsh_fit + double fourpi=4.0*constants::pi; + double fourmu2=4.0*gamma*gamma; + f2gfit.coeffs_=fourpi/fourmu2*coefficients; + GFit::prune_small_coefficients(eps,lo,hi,f2gfit.coeffs_,f2gfit.exponents_); + + if (prnt) { + print("fg fit"); + auto exact=[&gamma](const double r) -> double { + return 0.25/(gamma*gamma)*(1.0-2.0*exp(-gamma*r)+exp(-2.0*gamma*r))/r; + }; + f2gfit.print_accuracy(exact,lo,hi); + } + return f2gfit; + } + + /// return a fit for a general isotropic function /// note that the error is controlled over a uniform grid, the boundaries @@ -114,7 +282,26 @@ class GFit { /// return the exponents of the fit Tensor exponents() const {return exponents_;} - void truncate_periodic_expansion(Tensor& c, Tensor& e, + void static prune_small_coefficients(const double eps, const double lo, const double hi, + Tensor& coeff, Tensor& expnt) { + double mid = lo + (hi-lo)*0.5; + long npt=coeff.size(); + long i; + for (i=npt-1; i>0; --i) { + double cnew = coeff[i]*exp(-(expnt[i]-expnt[i-1])*mid*mid); + double errlo = coeff[i]*exp(-expnt[i]*lo*lo) - + cnew*exp(-expnt[i-1]*lo*lo); + double errhi = coeff[i]*exp(-expnt[i]*hi*hi) - + cnew*exp(-expnt[i-1]*hi*hi); + if (std::max(std::abs(errlo),std::abs(errhi)) > 0.03*eps) break; + npt--; + coeff[i-1] = coeff[i-1] + cnew; + } + coeff = coeff(Slice(0,npt-1)); + expnt = expnt(Slice(0,npt-1)); + } + + void truncate_periodic_expansion(Tensor& c, Tensor& e, double L, bool discardG0) const { double tcut = 0.25/L/L; @@ -157,6 +344,35 @@ class GFit { } + /// print coefficients and exponents, and values and errors + + /// @param[in] op the exact function, e.g. 1/r, exp(-mu r), etc + /// @param[in] lo lower bound for the range r + /// @param[in] hi higher bound for the range r + template + void print_accuracy(opT op, const double lo, const double hi) const { + + std::cout << "weights and roots" << std::endl; + for (int i=0; i coeffs_; @@ -173,11 +389,14 @@ class GFit { /// Multiresolution Quantum Chemistry in Multiwavelet Bases, /// Lecture Notes in Computer Science, vol. 2660, p. 707, 2003. static void bsh_fit(double mu, double lo, double hi, double eps, - Tensor& pcoeff, Tensor& pexpnt, bool prnt) { + Tensor& pcoeff, Tensor& pexpnt, bool prnt, bool fix_interval) { - if (mu < 0.0) throw "cannot handle negative mu in bsh_fit"; + if (mu < 0.0) throw "cannot handle negative mu in bsh_fit"; +// bool restrict_interval=(mu>0) and use_mu_for_restricting_interval; - if (mu > 0) { + + if ((mu > 0) and (not fix_interval)) { +// if (restrict_interval) { // Restrict hi according to the exponential decay double r = -log(4*constants::pi*0.01*eps); r = -log(r * 4*constants::pi*0.01*eps); @@ -195,7 +414,8 @@ class GFit { else if (eps >= 1e-12) TT = 26; else TT = 30; - if (mu > 0) { + if ((mu > 0) and (not fix_interval)) { +// if (restrict_interval) { slo = -0.5*log(4.0*TT/(mu*mu)); } else { @@ -253,21 +473,9 @@ class GFit { // end points ... if this error is less than the desired // precision, can discard the diffuse gaussian. - if (mu == 0.0) { - double mid = lo + (hi-lo)*0.5; - long i; - for (i=npt-1; i>0; --i) { - double cnew = coeff[i]*exp(-(expnt[i]-expnt[i-1])*mid*mid); - double errlo = coeff[i]*exp(-expnt[i]*lo*lo) - - cnew*exp(-expnt[i-1]*lo*lo); - double errhi = coeff[i]*exp(-expnt[i]*hi*hi) - - cnew*exp(-expnt[i-1]*hi*hi); - if (std::max(std::abs(errlo),std::abs(errhi)) > 0.03*eps) break; - npt--; - coeff[i-1] = coeff[i-1] + cnew; - } - coeff = coeff(Slice(0,npt-1)); - expnt = expnt(Slice(0,npt-1)); + if ((mu == 0.0) and (not fix_interval)) { +// if (restrict_interval) { + GFit::prune_small_coefficients(eps,lo,hi,coeff,expnt); } // Modify the coeffs of the largest exponents to satisfy the moment conditions @@ -328,27 +536,6 @@ class GFit { coeff(Slice(0,nmom-1)) = ncoeff; } - if (prnt) { - for (int i=0; i& pcoeff, Tensor& pexpnt, bool prnt) { + MADNESS_CHECK(gamma >0.0); // empirical number TT for the upper integration limit double TT; if (eps >= 1e-2) TT = 5; @@ -373,6 +561,7 @@ class GFit { else TT = 30; // integration limits for quadrature over s: slo and shi + // slo and shi must not depend on gamma!!! double slo=0.5 * log(eps) - 1.0; double shi=log(TT/(lo*lo))*0.5; @@ -401,32 +590,6 @@ class GFit { expnt[i] = 0.25*exp(-2.0*s); } - // Prune large exponents from the fit ... never necessary due to construction - - // Prune small exponents from Coulomb fit. Evaluate a gaussian at - // the range midpoint, and replace it there with the next most - // diffuse gaussian. Then examine the resulting error at the two - // end points ... if this error is less than the desired - // precision, can discard the diffuse gaussian. - - if (gamma == 0.0) { - double mid = lo + (hi-lo)*0.5; - long i; - for (i=npt-1; i>0; --i) { - double cnew = coeff[i]*exp(-(expnt[i]-expnt[i-1])*mid*mid); - double errlo = coeff[i]*exp(-expnt[i]*lo*lo) - - cnew*exp(-expnt[i-1]*lo*lo); - double errhi = coeff[i]*exp(-expnt[i]*hi*hi) - - cnew*exp(-expnt[i-1]*hi*hi); - if (std::max(std::abs(errlo),std::abs(errhi)) > 0.03*eps) break; - npt--; - coeff[i-1] = coeff[i-1] + cnew; - } - coeff = coeff(Slice(0,npt-1)); - expnt = expnt(Slice(0,npt-1)); - } - - if (prnt) { std::cout << "weights and roots for a Slater function with gamma=" << gamma << std::endl; for (int i=0; i& pcoeff, Tensor& pexpnt, bool prnt) { + Tensor coeff,expnt; + slater_fit(gamma, lo, hi, eps, coeff, expnt, prnt); + + pcoeff=Tensor(coeff.size()+1); + pcoeff(Slice(1,-1,1))=-coeff(_); + pexpnt=Tensor(expnt.size()+1); + pexpnt(Slice(1,-1,1))=expnt(_); + + pcoeff(0l)=1.0; + pexpnt(0l)=1.e-10; + } + + + /// fit a correlation factor f12^2 = (1- exp(-mu r))^2 = 1 - 2 exp(-mu r) + exp(-2 mu r) + + /// no factor 1/(2 mu) or square of it included! + /// use the Slater fit with an additional term: 1*exp(-0 r^2) + static void f12sq_fit(double gamma, double lo, double hi, double eps, + Tensor& pcoeff, Tensor& pexpnt, bool prnt) { + Tensor coeff1,expnt1, coeff2, expnt2; + slater_fit(gamma, lo, hi, eps, coeff1, expnt1, prnt); + slater_fit(2.0*gamma, lo, hi, eps, coeff2, expnt2, prnt); + + // check exponents are the same + MADNESS_CHECK((expnt1-expnt2).normf()/expnt1.size()<1.e-12); + + // add exponential terms + auto coeff=coeff2-2.0*coeff1; + pcoeff=Tensor(coeff1.size()+1); + pcoeff(Slice(1,-1,1))=coeff(_); + pexpnt=Tensor(expnt1.size()+1); + pexpnt(Slice(1,-1,1))=expnt1(_); + + // add constant term + pcoeff(0l)=1.0; + pexpnt(0l)=1.e-10; + } + + + void static bsh_fit_ndim(int ndim, double mu, double lo, double hi, double eps, Tensor& pcoeff, Tensor& pexpnt, bool prnt) { if (mu > 0) { diff --git a/src/madness/mra/mra.h b/src/madness/mra/mra.h index 4c3b2fb3875..0de6ba925d2 100644 --- a/src/madness/mra/mra.h +++ b/src/madness/mra/mra.h @@ -1252,6 +1252,7 @@ namespace madness { // if this and g are the same, use norm2() if (this->get_impl()==g.get_impl()) { + if (this->get_impl()->is_redundant()) this->get_impl()->undo_redundant(true); double norm=this->norm2(); return norm*norm; } @@ -1359,29 +1360,24 @@ namespace madness { /// g is constructed with an implicit multiplication, e.g. /// result = , with g = 1/r12 | gg> /// @param[in] g on-demand function - template - TENSOR_RESULT_TYPE(T,R) inner_on_demand(const Function& g) const { - MADNESS_ASSERT(g.is_on_demand() and (not this->is_on_demand())); - - this->reconstruct(); + template + TENSOR_RESULT_TYPE(T, R) inner_on_demand(const Function& g) const { + MADNESS_ASSERT(g.is_on_demand() and (not this->is_on_demand())); - // save for later, will be removed by make_Vphi - std::shared_ptr< FunctionFunctorInterface > func=g.get_impl()->get_functor(); - //leaf_op fnode_is_leaf(this->get_impl().get()); - Leaf_op_other fnode_is_leaf(this->get_impl().get()); - g.get_impl()->make_Vphi(fnode_is_leaf,true); // fence here + constexpr std::size_t LDIM=std::max(NDIM/2,std::size_t(1)); + auto func=dynamic_cast* >(g.get_impl()->get_functor().get()); + MADNESS_ASSERT(func); + func->make_redundant(true); + func->replicate_low_dim_functions(true); + this->reconstruct(); // if this == &g we don't need g to be redundant - if (VERIFY_TREE) verify_tree(); - TENSOR_RESULT_TYPE(T,R) local = impl->inner_local(*g.get_impl()); - impl->world.gop.sum(local); - impl->world.gop.fence(); + if (VERIFY_TREE) verify_tree(); - // restore original state - g.get_impl()->set_functor(func); - g.get_impl()->get_coeffs().clear(); - g.get_impl()->set_tree_state(on_demand); + TENSOR_RESULT_TYPE(T, R) local = impl->inner_local_on_demand(*g.get_impl()); + impl->world.gop.sum(local); + impl->world.gop.fence(); - return local; + return local; } /// project this on the low-dim function g: h(x) = @@ -2149,7 +2145,6 @@ namespace madness { if (op.modified()) { - MADNESS_ASSERT(not op.is_slaterf12); ff.get_impl()->make_redundant(true); result = apply_only(op, ff, fence); ff.get_impl()->undo_redundant(false); @@ -2157,14 +2152,6 @@ namespace madness { } else { - // the slaterf12 function is - // 1/(2 mu) \int d1 (1 - exp(- mu r12)) f(1) - // = 1/(2 mu) (f.trace() - \int d1 exp(-mu r12) f(1) ) - // f.trace() is just a number - R ftrace=0.0; - if (op.is_slaterf12) ftrace=f.trace(); -// print("ftrace",ftrace); - // saves the standard() step, which is very expensive in 6D // Function fff=copy(ff); Function fff=(ff); @@ -2196,7 +2183,6 @@ namespace madness { } else { ff.standard(); } - if (op.is_slaterf12) result=(result-ftrace).scale(-0.5/op.mu()); } if (print_timings) result.print_size("result after reconstruction"); diff --git a/src/madness/mra/mraimpl.h b/src/madness/mra/mraimpl.h index c748a640c91..dd3b155f038 100644 --- a/src/madness/mra/mraimpl.h +++ b/src/madness/mra/mraimpl.h @@ -960,6 +960,11 @@ namespace madness { /// Returns true if this block of coeffs needs autorefining template bool FunctionImpl::autorefine_square_test(const keyT& key, const nodeT& t) const { + // Chosen approach looks stupid but it is more accurate + // than the simple approach of summing everything and + // subtracting off the low-order stuff to get the high + // order (assuming the high-order stuff is small relative + // to the low-order) double lo, hi; tnorm(t.coeff().full_tensor_copy(), &lo, &hi); double test = 2*lo*hi + hi*hi; @@ -3035,13 +3040,9 @@ namespace madness { template - void FunctionImpl::tnorm(const tensorT& t, double* lo, double* hi) const { + void FunctionImpl::tnorm(const tensorT& t, double* lo, double* hi) { //PROFILE_MEMBER_FUNC(FunctionImpl); // Too fine grain for routine profiling - // Chosen approach looks stupid but it is more accurate - // than the simple approach of summing everything and - // subtracting off the low-order stuff to get the high - // order (assuming the high-order stuff is small relative - // to the low-order) + auto& cdata=FunctionCommonData::get(t.dim(0)); tensorT work = copy(t); tensorT tlo = work(cdata.sh); *lo = tlo.normf(); @@ -3050,7 +3051,8 @@ namespace madness { } template - void FunctionImpl::tnorm(const GenTensor& t, double* lo, double* hi) const { + void FunctionImpl::tnorm(const GenTensor& t, double* lo, double* hi) { + auto& cdata=FunctionCommonData::get(t.dim(0)); coeffT shalf=t(cdata.sh); *lo=shalf.normf(); coeffT sfull=copy(t); @@ -3060,9 +3062,10 @@ namespace madness { template void FunctionImpl::tnorm(const SVDTensor& t, double* lo, double* hi, - const int particle) const { + const int particle) { *lo=0.0; *hi=0.0; + auto& cdata=FunctionCommonData::get(t.dim(0)); if (t.rank()==0) return; const tensorT vec=t.flat_vector(particle-1); for (long i=0; i + class SeparatedConvolution; + + class CCPairFunction; + template + std::vector apply(const SeparatedConvolution& op, const std::vector& argument); + + template + std::vector< Function > + apply(const SeparatedConvolution& op, const std::vector< Function > f); + + + template + CCPairFunction apply(const SeparatedConvolution& op, const CCPairFunction& argument); + /// SeparatedConvolutionInternal keeps data for 1 term and all dimensions and 1 displacement /// Why is this here?? Why don't you just use ConvolutionND in SeparatedConvolutionData?? template @@ -110,17 +125,46 @@ namespace madness { */ + /// operator types + enum OpType { + OT_UNDEFINED, + OT_ONE, /// indicates the identity + OT_G12, /// 1/r + OT_SLATER, /// exp(-r) + OT_GAUSS, /// exp(-r2) + OT_F12, /// 1-exp(-r) + OT_FG12, /// (1-exp(-r))/r + OT_F212, /// (1-exp(-r))^2 + OT_F2G12, /// (1-exp(-r))^2/r = 1/r + exp(-2r)/r - 2 exp(-r)/r + OT_BSH /// exp(-r)/r + }; + + struct OperatorInfo { + OperatorInfo() = default; +// template +// OperatorInfo(double mu, double lo=1.e-5, double thresh=FunctionDefaults::get_thresh()) : mu(mu), thresh(thresh), lo(lo) {} + OperatorInfo(double mu, double lo, double thresh, OpType type) : mu(mu), thresh(thresh), lo(lo), type(type) {} + OpType type=OT_UNDEFINED; ///< introspection + double mu=0.0; ///< some introspection + double thresh=1.e-4; + double lo=1.e-5; + }; template class SeparatedConvolution : public WorldObject< SeparatedConvolution > { public: + typedef Q opT; ///< The apply function uses this to infer resultT=opT*inputT + + OperatorInfo info; + bool doleaves; ///< If should be applied to leaf coefficients ... false by default bool isperiodicsum;///< If true the operator 1D kernels have been summed over lattice translations ///< and may be non-zero at both ends of the unit cell - bool modified_; ///< use modified NS form - int particle_; - bool destructive_; ///< destroy the argument or restore it (expensive for 6d functions) + bool modified_=false; ///< use modified NS form + int particle_=1; ///< must only be 1 or 2 + bool destructive_=false; ///< destroy the argument or restore it (expensive for 6d functions) + bool print_timings=false; typedef Key keyT; const static size_t opdim=NDIM; @@ -129,11 +173,6 @@ namespace madness { Timer timer_low_accumulate; Timer timer_stats_accumulate; - // if this is a Slater-type convolution kernel: 1-exp(-mu r12)/(2 mu) - bool is_slaterf12; - bool print_timings=false; - double mu_; - private: @@ -141,7 +180,7 @@ namespace madness { const BoundaryConditions bc; const int k; const FunctionCommonData& cdata; - const int rank; + int rank; const std::vector vk; const std::vector v2k; const std::vector s0; @@ -161,8 +200,8 @@ namespace madness { bool& destructive() {return destructive_;} const bool& destructive() const {return destructive_;} - const double& gamma() const {return mu_;} - const double& mu() const {return mu_;} + const double& gamma() const {return info.mu;} + const double& mu() const {return info.mu;} private: @@ -181,6 +220,43 @@ namespace madness { const Q* VT; }; + static inline std::pair,Tensor> + make_coeff_for_operator(World& world, double mu, double lo, double eps, OpType type, + const BoundaryConditions& bc=FunctionDefaults::get_bc()) { + + const Tensor& cell_width = FunctionDefaults<3>::get_cell_width(); + double hi = cell_width.normf(); // Diagonal width of cell + if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation + + GFit fit; + if (type==OT_G12) {fit=GFit::CoulombFit(lo,hi,eps,false); + } else if (type==OT_SLATER) {fit=GFit::SlaterFit(mu,lo,hi,eps,false); + } else if (type==OT_GAUSS) {fit=GFit::GaussFit(mu,lo,hi,eps,false); + } else if (type==OT_F12) {fit=GFit::F12Fit(mu,lo,hi,eps,false); + } else if (type==OT_FG12) {fit=GFit::FGFit(mu,lo,hi,eps,false); + } else if (type==OT_F212) {fit=GFit::F12sqFit(mu,lo,hi,eps,false); + } else if (type==OT_F2G12) {fit=GFit::F2GFit(mu,lo,hi,eps,false); + } else if (type==OT_BSH) {fit=GFit::BSHFit(mu,lo,hi,eps,false); + } else { + MADNESS_EXCEPTION("Operator type not implemented",1); + } + + Tensor coeff=fit.coeffs(); + Tensor expnt=fit.exponents(); + + if (bc(0,0) == BC_PERIODIC) { + fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); + } + + return std::make_pair(coeff,expnt); + } + + static inline std::pair,Tensor> + make_coeff_for_operator(World& world, OperatorInfo info, + const BoundaryConditions& bc=FunctionDefaults::get_bc()) { + return make_coeff_for_operator(world,info.mu,info.lo,info.thresh,info.type,bc); + } + // /// return the right block of the upsampled operator (modified NS only) // // /// unlike the operator matrices on the natural level the upsampled operator @@ -882,8 +958,7 @@ namespace madness { , modified_(false) , particle_(1) , destructive_(false) - , is_slaterf12(false) - , mu_(0.0) + , info() , bc(bc) , k(k) , cdata(FunctionCommonData::get(k)) @@ -917,8 +992,7 @@ namespace madness { , modified_(false) , particle_(1) , destructive_(false) - , is_slaterf12(false) - , mu_(0.0) + , info() , ops(argops) , bc(bc) , k(k) @@ -935,9 +1009,23 @@ namespace madness { this->process_pending(); } + /// Constructor for Gaussian Convolutions (mostly for backward compatability) + SeparatedConvolution(World& world, const OperatorInfo info1, + const BoundaryConditions& bc = FunctionDefaults::get_bc(), + int k=FunctionDefaults::get_k(), + bool doleaves = false) + : SeparatedConvolution(world,Tensor(0l),Tensor(0l),info1.lo,info1.thresh,bc,k,doleaves,info1.mu) { + info.type=info1.type; + auto [coeff, expnt] =make_coeff_for_operator(world, info1, bc); + rank=coeff.dim(0); + ops.resize(rank); + initialize(coeff,expnt); + } + /// Constructor for Gaussian Convolutions (mostly for backward compatability) SeparatedConvolution(World& world, const Tensor& coeff, const Tensor& expnt, + double lo, double thresh, const BoundaryConditions& bc = FunctionDefaults::get_bc(), int k=FunctionDefaults::get_k(), bool doleaves = false, @@ -945,11 +1033,7 @@ namespace madness { : WorldObject< SeparatedConvolution >(world) , doleaves(doleaves) , isperiodicsum(bc(0,0)==BC_PERIODIC) - , modified_(false) - , particle_(1) - , destructive_(false) - , is_slaterf12(mu>0.0) - , mu_(mu) + , info(mu,lo,thresh,OT_UNDEFINED) , ops(coeff.dim(0)) , bc(bc) , k(k) @@ -957,8 +1041,11 @@ namespace madness { , rank(coeff.dim(0)) , vk(NDIM,k) , v2k(NDIM,2*k) - , s0(std::max(2,NDIM),Slice(0,k-1)) - { + , s0(std::max(2,NDIM),Slice(0,k-1)) { + initialize(coeff,expnt); + } + + void initialize(const Tensor& coeff, const Tensor& expnt) { // Presently we must have periodic or non-periodic in all dimensions. for (std::size_t d=1; d + std::vector> operator()(const std::vector>& f) const { + return madness::apply(*this, f); + } + /// apply this operator on a separable function f(1,2) = f(1) f(2) /// @param[in] f1 a function of dim LDIM @@ -1105,7 +1197,6 @@ namespace madness { template Function operator()(const Function& f1, const Function& f2) const { - MADNESS_ASSERT(not is_slaterf12); return madness::apply(*this, std::vector>({f1}), std::vector>({f2})); } @@ -1118,10 +1209,17 @@ namespace madness { template Function operator()(const std::vector>& f1, const std::vector>& f2) const { - MADNESS_ASSERT(not is_slaterf12); return madness::apply(*this, f1, f2); } + /// apply this onto another suitable argument, returning the same type + + /// argT must implement argT::apply(const SeparatedConvolution& op, const argT& arg) + template + argT operator()(const argT& argument) const { + return madness::apply(*this,argument); + } + /// apply this operator on coefficients in full rank form @@ -1551,6 +1649,73 @@ namespace madness { return tt; } + + static bool can_combine(const SeparatedConvolution& left, const SeparatedConvolution& right) { + return (combine_OT(left,right).type!=OT_UNDEFINED); + } + + /// return operator type and other info of the combined operator (e.g. fg = f(1,2)* g(1,2) + static OperatorInfo combine_OT(const SeparatedConvolution& left, const SeparatedConvolution& right) { + OperatorInfo info=left.info; + if ((left.info.type==OT_F12) and (right.info.type==OT_G12)) { + info.type=OT_FG12; + } else if ((left.info.type==OT_GAUSS) and (right.info.type==OT_GAUSS)) { + info=right.info; + info.type=OT_GAUSS; + info.mu=2.0*right.info.mu; + } else if ((left.info.type==OT_SLATER) and (right.info.type==OT_SLATER)) { + info=right.info; + info.type=OT_SLATER; + info.mu=2.0*right.info.mu; + } else if ((left.info.type==OT_G12) and (right.info.type==OT_F12)) { + info=right.info; + info.type=OT_FG12; + } else if ((left.info.type==OT_G12) and (right.info.type==OT_F212)) { + info=right.info; + info.type=OT_F2G12; + } else if (((left.info.type==OT_F212) and (right.info.type==OT_G12)) or + ((left.info.type==OT_F12) and (right.info.type==OT_FG12)) or + ((left.info.type==OT_FG12) and (right.info.type==OT_F12))) { + info=right.info; + info.type=OT_F2G12; + if (right.info.type!=OT_G12) MADNESS_CHECK(right.info.mu == left.info.mu); + } else if ((left.info.type==OT_F12) and (right.info.type==OT_F12)) { + info.type=OT_F212; + // keep the original gamma + // (f12)^2 = (1- slater12)^2 = 1/(4 gamma) (1 - 2 exp(-gamma) + exp(-2 gamma)) + MADNESS_CHECK(right.info.mu == left.info.mu); + } else { + MADNESS_EXCEPTION("unknown combination of SeparatedConvolutions: feel free to extend in operator.h",1); + } + return info; + } + + + /// combine 2 convolution operators to one + static SeparatedConvolution combine(const SeparatedConvolution& left, + const SeparatedConvolution& right) { + MADNESS_CHECK(can_combine(left,right)); + MADNESS_CHECK(left.get_world().id()==right.get_world().id()); + + auto info=combine_OT(left,right); + return SeparatedConvolution(left.get_world(),info,left.bc,left.k); + } + + /// combine 2 convolution operators to one + friend SeparatedConvolution combine(const std::shared_ptr> left, + const std::shared_ptr> right) { + SeparatedConvolution result; + if (left and right) { + return combine(*left, *right); + } else if (left) { + return *left; + } else if (right) { + return *right; + } else { + MADNESS_EXCEPTION("can't combine empty SeparatedConvolutions",1); + } + return result; + } }; @@ -1604,7 +1769,7 @@ namespace madness { if (bc(0,0) == BC_PERIODIC) { fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true); } - return SeparatedConvolution(world, coeff, expnt, bc, k); + return SeparatedConvolution(world, coeff, expnt, lo, eps, bc, k); } @@ -1627,94 +1792,45 @@ namespace madness { if (bc(0,0) == BC_PERIODIC) { fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true); } - return new SeparatedConvolution(world, coeff, expnt, bc, k); + return new SeparatedConvolution(world, coeff, expnt, lo, eps, bc, k); } /// Factory function generating separated kernel for convolution with BSH kernel in general NDIM template - static - inline - SeparatedConvolution BSHOperator(World& world, - double mu, - double lo, - double eps, - const BoundaryConditions& bc=FunctionDefaults::get_bc(), - int k=FunctionDefaults::get_k()) - { + static inline + SeparatedConvolution + BSHOperator(World& world, double mu, double lo, double eps, + const BoundaryConditions& bc=FunctionDefaults::get_bc(), + int k=FunctionDefaults::get_k()) { if (eps>1.e-4) { if (world.rank()==0) print("the accuracy in BSHOperator is too small, tighten the threshold",eps); MADNESS_EXCEPTION("0",1); } - const Tensor& cell_width = FunctionDefaults::get_cell_width(); - double hi = cell_width.normf(); // Diagonal width of cell - if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation - - GFit fit=GFit::BSHFit(mu,lo,hi,eps,false); - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); - - if (bc(0,0) == BC_PERIODIC) { - fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); - } - - return SeparatedConvolution(world, coeff, expnt, bc, k); + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_BSH),bc,k); } /// Factory function generating separated kernel for convolution with BSH kernel in general NDIM template - static - inline - SeparatedConvolution* BSHOperatorPtr(World& world, - double mu, - double lo, - double eps, - const BoundaryConditions& bc=FunctionDefaults::get_bc(), - int k=FunctionDefaults::get_k()) - { - if (eps>1.e-4) { - if (world.rank()==0) print("the accuracy in BSHOperator is too small, tighten the threshold",eps); - MADNESS_EXCEPTION("0",1); - } - const Tensor& cell_width = FunctionDefaults::get_cell_width(); - double hi = cell_width.normf(); // Diagonal width of cell - if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation - - GFit fit=GFit::BSHFit(mu,lo,hi,eps,false); - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); - - if (bc(0,0) == BC_PERIODIC) { - fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); + static inline + SeparatedConvolution* + BSHOperatorPtr(World& world, double mu, double lo, double eps, + const BoundaryConditions& bc=FunctionDefaults::get_bc(), + int k=FunctionDefaults::get_k()) { + if (eps>1.e-4) { + if (world.rank()==0) print("the accuracy in BSHOperator is too small, tighten the threshold",eps); + MADNESS_EXCEPTION("0",1); } - - return new SeparatedConvolution(world, coeff, expnt, bc, k); + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_BSH),bc,k); } /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D - static - inline - SeparatedConvolution BSHOperator3D(World& world, - double mu, - double lo, - double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), - int k=FunctionDefaults<3>::get_k()) - - { - const Tensor& cell_width = FunctionDefaults<3>::get_cell_width(); - double hi = cell_width.normf(); // Diagonal width of cell - if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation - - GFit fit=GFit::BSHFit(mu,lo,hi,eps,false); - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); - - if (bc(0,0) == BC_PERIODIC) { - fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); - } - return SeparatedConvolution(world, coeff, expnt, bc, k); + static inline SeparatedConvolution + BSHOperator3D(World& world, double mu, double lo, double eps, + const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + int k=FunctionDefaults<3>::get_k()) { + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_BSH),bc,k); } /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D @@ -1744,8 +1860,7 @@ namespace madness { } /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D - static - inline + static inline SeparatedConvolution* PeriodicBSHOperatorPtr3D(World& world, Vector args, double mu, @@ -1769,49 +1884,59 @@ namespace madness { return new SeparatedConvolution(world, args, coeff, expnt, bc, k); } - /// Factory function generating separated kernel for convolution with (1 - exp(-mu*r))/(2 mu) in 3D - /// note that the 1/2mu factor is not included here, nor is the term 1/(2mu) - static inline SeparatedConvolution SlaterF12Operator(World& world, - double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), - int k=FunctionDefaults<3>::get_k()) { - - const Tensor& cell_width = FunctionDefaults<3>::get_cell_width(); - double hi = cell_width.normf(); // Diagonal width of cell - if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation - - GFit fit=GFit::SlaterFit(mu,lo,hi,eps,false); - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); + static inline SeparatedConvolution + SlaterF12Operator(World& world, double mu, double lo, double eps, + const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), int k=FunctionDefaults<3>::get_k()) { + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F12),bc,k); + } - if (bc(0,0) == BC_PERIODIC) { - fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); - } - return SeparatedConvolution(world, coeff, expnt, bc, k, false, mu); + static inline SeparatedConvolution SlaterF12sqOperator(World& world, + double mu, double lo, double eps, + const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + int k=FunctionDefaults<3>::get_k()) { + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F212),bc,k); } /// Factory function generating separated kernel for convolution with exp(-mu*r) in 3D - /// Note that the 1/(2mu) factor is not included, this is just the exponential function static inline SeparatedConvolution SlaterOperator(World& world, double mu, double lo, double eps, const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), int k=FunctionDefaults<3>::get_k()) { + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_SLATER),bc,k); + } - const Tensor& cell_width = FunctionDefaults<3>::get_cell_width(); - double hi = cell_width.normf(); // Diagonal width of cell - if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation + /// Factory function generating separated kernel for convolution with exp(-mu*r*r) - GFit fit=GFit::SlaterFit(mu,lo,hi,eps,false); - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); + /// lo and eps are not used here + template + static inline SeparatedConvolution GaussOperator(World& world, + double mu, double lo=0.0, double eps=0.0, + const BoundaryConditions& bc=FunctionDefaults::get_bc(), + int k=FunctionDefaults::get_k()) { + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_GAUSS),bc,k); + } - if (bc(0,0) == BC_PERIODIC) { - fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); - } - SeparatedConvolution tmp(world, coeff, expnt, bc, k, false, mu); - tmp.is_slaterf12 = false; - return tmp; + /// Factory function generating separated kernel for convolution with exp(-mu*r*r) in 3D + + /// lo and eps are not used here + template + static inline SeparatedConvolution* GaussOperatorPtr(World& world, + double mu, double lo=0.0, double eps=0.0, + const BoundaryConditions& bc = FunctionDefaults::get_bc(), + int k = FunctionDefaults::get_k()) { + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_GAUSS),bc,k); + } + + + /// Factory function generating separated kernel for convolution with exp(-mu*r) in 3D + /// Note that the 1/(2mu) factor of SlaterF12Operator is not included, this is just the exponential function + template + static inline SeparatedConvolution* SlaterOperatorPtr_ND(World& world, + double mu, double lo, double eps, + const BoundaryConditions& bc = FunctionDefaults::get_bc(), + int k = FunctionDefaults::get_k()) { + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_SLATER),bc,k); } /// Factory function generating separated kernel for convolution with exp(-mu*r) in 3D @@ -1820,41 +1945,64 @@ namespace madness { double mu, double lo, double eps, const BoundaryConditions<3>& bc = FunctionDefaults<3>::get_bc(), int k = FunctionDefaults<3>::get_k()) { - - const Tensor& cell_width = FunctionDefaults<3>::get_cell_width(); - double hi = cell_width.normf(); // Diagonal width of cell - if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation - - GFit fit=GFit::SlaterFit(mu,lo,hi,eps,false); - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); - - if (bc(0,0) == BC_PERIODIC) { - fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); - } - SeparatedConvolution* tmp=new SeparatedConvolution(world, coeff, expnt, bc, k, false, mu); - tmp->is_slaterf12 = false; - return tmp; + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_SLATER),bc,k); } /// Factory function generating separated kernel for convolution with (1 - exp(-mu*r))/(2 mu) in 3D + + /// includes the factor 1/(2 mu) static inline SeparatedConvolution* SlaterF12OperatorPtr(World& world, double mu, double lo, double eps, const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), int k=FunctionDefaults<3>::get_k()) { + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F12),bc,k); + } - const Tensor& cell_width = FunctionDefaults<3>::get_cell_width(); - double hi = cell_width.normf(); // Diagonal width of cell - if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation - GFit fit=GFit::SlaterFit(mu,lo,hi,eps,false); - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); + /// Factory function generating separated kernel for convolution with 1/(2 mu)*(1 - exp(-mu*r))/r in 3D - if (bc(0,0) == BC_PERIODIC) { - fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); - } - return new SeparatedConvolution(world, coeff, expnt, bc, k, false, mu); + /// fg = (1 - exp(-gamma r12)) / r12 = 1/r12 - exp(-gamma r12)/r12 = coulomb - bsh + /// includes the factor 1/(2 mu) + static inline SeparatedConvolution + FGOperator(World& world, double mu, double lo, double eps, + const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + int k=FunctionDefaults<3>::get_k()) { + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_FG12),bc,k); + } + + /// Factory function generating separated kernel for convolution with 1/(2 mu)*(1 - exp(-mu*r))/r in 3D + + /// fg = (1 - exp(-gamma r12)) / r12 = 1/r12 - exp(-gamma r12)/r12 = coulomb - bsh + /// includes the factor 1/(2 mu) + static inline SeparatedConvolution* + FGOperatorPtr(World& world, double mu, double lo, double eps, + const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + int k=FunctionDefaults<3>::get_k()) { + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_FG12),bc,k); + } + + /// Factory function generating separated kernel for convolution with (1/(2 mu)*(1 - exp(-mu*r)))^2/r in 3D + + /// f2g = (1/(2 gamma) (1 - exp(-gamma r12)))^2 / r12 + /// = 1/(4 gamma) * [ 1/r12 - 2 exp(-gamma r12)/r12 + exp(-2 gamma r12)/r12 ] + /// includes the factor 1/(2 mu)^2 + static inline SeparatedConvolution* + F2GOperatorPtr(World& world, double mu, double lo, double eps, + const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + int k=FunctionDefaults<3>::get_k()) { + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F2G12),bc,k); + } + + /// Factory function generating separated kernel for convolution with (1/(2 mu)*(1 - exp(-mu*r)))^2/r in 3D + + /// f2g = (1/(2 gamma) (1 - exp(-gamma r12)))^2 / r12 + /// = 1/(4 gamma) * [ 1/r12 - 2 exp(-gamma r12)/r12 + exp(-2 gamma r12)/r12 ] + /// includes the factor 1/(2 mu)^2 + static inline SeparatedConvolution + F2GOperator(World& world, double mu, double lo, double eps, + const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + int k=FunctionDefaults<3>::get_k()) { + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F2G12),bc,k); } @@ -1869,7 +2017,7 @@ namespace madness { Tensor coeffs(1), exponents(1); exponents(0L) = exponent; coeffs(0L)=pow(exponent/M_PI,0.5*3.0); // norm of the gaussian - return SeparatedConvolution(world, coeffs, exponents); + return SeparatedConvolution(world, coeffs, exponents, 1.e-8, eps); } @@ -1885,7 +2033,7 @@ namespace madness { Tensor coeffs(1), exponents(1); exponents(0L) = exponent; coeffs(0L)=pow(exponent/M_PI,0.5*NDIM); // norm of the gaussian - return SeparatedConvolution(world, coeffs, exponents); + return SeparatedConvolution(world, coeffs, exponents, 1.e-8, eps); } /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D @@ -1909,7 +2057,7 @@ namespace madness { if (bc(0,0) == BC_PERIODIC) { fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); } - return new SeparatedConvolution(world, coeff, expnt, bc, k); + return new SeparatedConvolution(world, coeff, expnt, lo, eps, bc, k); } @@ -2021,6 +2169,7 @@ namespace madness { } + namespace archive { template struct ArchiveLoadImpl*> { diff --git a/src/madness/mra/test6.cc b/src/madness/mra/test6.cc index ee05cd4f85e..6b3cdc5ce8e 100644 --- a/src/madness/mra/test6.cc +++ b/src/madness/mra/test6.cc @@ -36,6 +36,7 @@ /// \brief test various functionality for 6d functions #include +#include using namespace madness; @@ -158,9 +159,8 @@ static double V(const Vector& r) { /// test f(1,2) = g(1) h(2) int test_hartree_product(World& world, const long& k, const double thresh) { - - print("entering hartree_product"); - int nerror=0; + test_output output("testing hartree_product"); +// output.set_cout_to_terminal(); bool good; real_function_3d phi=real_factory_3d(world).f(gauss_3d); @@ -168,28 +168,31 @@ int test_hartree_product(World& world, const long& k, const double thresh) { { real_function_6d ij=hartree_product(phi,phi); + ij.print_size("ij before truncation"); ij.truncate(); + ij.print_size("ij after truncation"); double norm=ij.norm2(); - print("norm(ij)",norm); - double err=ij.err(gauss_6d); - good=is_small(err,thresh); - print(ok(good), "hartree_product(phi,phi) error:",err); - if (not good) nerror++; + print("norm(ij), error, thresh",norm, err, thresh); + + good=is_small(err,thresh*2.5); + output.checkpoint(good, "hartree_product(phi,phi)"); } { - real_function_6d ij=hartree_product(phisq,phi); - double err=ij.err(r2r); - good=is_small(err,thresh); - print(ok(good), "hartree_product(phi^2,phi) error:",err); - if (not good) nerror++; + real_function_6d iij=hartree_product(phisq,phi); + iij.print_size("iij before truncation"); + double norm=iij.norm2(); + double err=iij.err(r2r); + print("norm(iij), error, thresh",norm, err, thresh); + good=is_small(err,thresh*2.0); + output.checkpoint(good, "hartree_product(phi^2,phi)"); } print("all done\n"); - return nerror; + return (output.get_final_success()) ? 0 : 1; } /// test f(1,2)*g(1) @@ -643,8 +646,9 @@ int test(World& world, const long& k, const double thresh) { int main(int argc, char**argv) { - initialize(argc,argv); - World world(SafeMPI::COMM_WORLD); +// initialize(argc,argv); +// World world(SafeMPI::COMM_WORLD); + World& world= initialize(argc,argv); srand(time(nullptr)); startup(world,argc,argv); @@ -678,13 +682,14 @@ int main(int argc, char**argv) { } } - FunctionDefaults<3>::set_thresh(thresh); + FunctionDefaults<3>::set_thresh(thresh*0.1); FunctionDefaults<3>::set_k(k); FunctionDefaults<3>::set_cubic_cell(-L/2,L/2); FunctionDefaults<6>::set_thresh(thresh); FunctionDefaults<6>::set_k(k); FunctionDefaults<6>::set_cubic_cell(-L/2,L/2); FunctionDefaults<6>::set_tensor_type(tt); + FunctionDefaults<6>::set_truncate_mode(3); print("entering testsuite for 6-dimensional functions\n"); print("k ",k); @@ -696,26 +701,28 @@ int main(int argc, char**argv) { int error=0; - real_function_3d phi=real_factory_3d(world).f(gauss_3d); - double norm=phi.norm2(); - if (world.rank()==0) printf("phi.norm2() %12.8f\n",norm); - - real_function_3d phi2=2.0*phi*phi; - norm=phi2.norm2(); - if (world.rank()==0) printf("phi2.norm2() %12.8f\n",norm); - - test(world,k,thresh); - error+=test_hartree_product(world,k,thresh); - error+=test_convolution(world,k,thresh); - error+=test_multiply(world,k,thresh); - error+=test_add(world,k,thresh); - error+=test_exchange(world,k,thresh); - error+=test_inner(world,k,thresh); - error+=test_replicate(world,k,thresh); - - print(ok(error==0),error,"finished test suite\n"); + { + real_function_3d phi=real_factory_3d(world).f(gauss_3d); + double norm=phi.norm2(); + if (world.rank()==0) printf("phi.norm2() %12.8f\n",norm); + + real_function_3d phi2=2.0*phi*phi; + norm=phi2.norm2(); + if (world.rank()==0) printf("phi2.norm2() %12.8f\n",norm); + +// test(world,k,thresh); + error+=test_hartree_product(world,k,thresh); +// error+=test_convolution(world,k,thresh); +// error+=test_multiply(world,k,thresh); +// error+=test_add(world,k,thresh); +// error+=test_exchange(world,k,thresh); +// error+=test_inner(world,k,thresh); +// error+=test_replicate(world,k,thresh); + + print(ok(error==0),error,"finished test suite\n"); + world.gop.fence(); + } - world.gop.fence(); finalize(); return error; diff --git a/src/madness/mra/testinnerext.cc b/src/madness/mra/testinnerext.cc index c086d3c71b7..b2bd6e0f28a 100644 --- a/src/madness/mra/testinnerext.cc +++ b/src/madness/mra/testinnerext.cc @@ -82,12 +82,41 @@ bool test_loose1(std::string msg, double a, double b, double tol=thresh) { } +int test_tight_diffuse(World& world) { + FunctionDefaults<4>::set_thresh(1.e-5); + double a=1.e2; + double b=1.e-2; + for (int i=0; i<4; ++i) { + a=std::pow(10.0,double(i)); + b=std::pow(0.1,double(i)); + print("a,b",a,b); + + real_function_2d aa=real_factory_2d(world).functor([&](const coord_2d& r){return exp(-a*inner(r,r));}); + real_function_2d bb=real_factory_2d(world).functor([&](const coord_2d& r){return exp(-b*inner(r,r));}); + real_function_4d cc=real_factory_4d(world).functor([&](const coord_4d& r){return exp(-b*inner(r,r));}); + real_function_4d dd=CompositeFactory(world).particle1(aa).particle2(copy(aa)); + aa.print_size("exp(-1000 r^2"); + bb.print_size("exp(-0.001 r^2"); + double result=inner(cc,dd); + double refresult=std::pow(constants::pi/(a+b),2.0); + print("result,refresult,error",result,refresult,result-refresult); + MADNESS_CHECK(test(" inner(exp(-a r^2 , exp(-b r^2)) ", result,refresult)); + } + + + + + return 0; + + +} + int test_partial_inner(World& world) { - print("\ntesting partial inner\n"); bool do_low_rank=false; #if HAVE_GENTENSOR do_low_rank=true; #endif + print("\ntesting partial inner; low rank: ",do_low_rank,"\n"); real_function_1d one_1d=real_factory_1d(world).functor([](const coord_1d& r){return 1.0;}); real_function_2d one_2d=real_factory_2d(world).functor([](const coord_2d& r){return 1.0;}); @@ -268,6 +297,7 @@ int main(int argc, char** argv) { initialize<6>(world); test_partial_inner(world); + test_tight_diffuse(world); if (!smalltest) { real_function_3d alpha1 = real_factory_3d(world).f(alpha_func); diff --git a/src/madness/mra/testper.cc b/src/madness/mra/testper.cc index 61ba93b6c96..a044756375e 100644 --- a/src/madness/mra/testper.cc +++ b/src/madness/mra/testper.cc @@ -65,7 +65,8 @@ int test_per(World& world) { expnt[0] = 10000.0; coeff[0] = sqrt(expnt[0]/constants::pi); print(coeff,expnt); - SeparatedConvolution op(world, coeff, expnt); + double lo=1.e-6; thresh=1.e-4; // dummy values + SeparatedConvolution op(world, coeff, expnt, lo, thresh); Function f = FunctionFactory(world).f(constant).initial_level(3).norefine(); diff --git a/src/madness/mra/testsuite.cc b/src/madness/mra/testsuite.cc index faf62fb3e45..83251edd879 100644 --- a/src/madness/mra/testsuite.cc +++ b/src/madness/mra/testsuite.cc @@ -743,7 +743,8 @@ int test_op(World& world) { Tensor coeffs(1), exponents(1); exponents(0L) = 10.0; coeffs(0L) = pow(exponents(0L)/PI, 0.5*NDIM); - SeparatedConvolution op(world, coeffs, exponents); + double lo=1.e-6, thresh1=1.e-4; + SeparatedConvolution op(world, coeffs, exponents, lo, thresh1); START_TIMER; Function r = madness::apply(op,f); END_TIMER("apply"); diff --git a/src/madness/mra/vmra.h b/src/madness/mra/vmra.h index b8f291b72a9..b24819fd300 100644 --- a/src/madness/mra/vmra.h +++ b/src/madness/mra/vmra.h @@ -813,7 +813,9 @@ namespace madness { { world.gop.fence(); compress(world, f); - if ((void*)(&f) != (void*)(&g)) compress(world, g); +// if ((void*)(&f) != (void*)(&g)) compress(world, g); + compress(world, g); + std::vector*> left(f.size()); std::vector*> right(g.size()); @@ -1043,7 +1045,45 @@ namespace madness { } - /// Computes the square of a vector of functions --- q[i] = v[i]**2 + /// multiply a high-dimensional function with a low-dimensional function + + /// @param[in] f NDIM function of NDIM dimensions + /// @param[in] g LDIM function of LDIM + /// @param[in] v dimension indices of f to multiply + /// @return h[i](0,1,2,3) = f(0,1,2,3) * g[i](1,2,3) for v={1,2,3} + template + std::vector > partial_mul(const Function f, const std::vector > g, + const int particle) { + + World& world=f.world(); + std::vector > result(g.size()); + for (auto& r : result) r.set_impl(f, false); + + FunctionImpl* fimpl=f.get_impl().get(); + fimpl->make_redundant(false); + make_redundant(world,g,false); + world.gop.fence(); + + for (std::size_t i=0; i* gimpl=g[i].get_impl().get(); + result[i].get_impl()->multiply(fimpl,gimpl,particle); // stupid naming inconsistency + } + world.gop.fence(); + + fimpl->undo_redundant(false); + for (auto& ig : g) ig.get_impl()->undo_redundant(false); + world.gop.fence(); + return result; + } + + template + std::vector > multiply(const Function f, const std::vector > g, + const std::tuple v) { + return partial_mul(f,g,std::array({std::get<0>(v),std::get<1>(v),std::get<2>(v)})); + } + + +/// Computes the square of a vector of functions --- q[i] = v[i]**2 template std::vector< Function > square(World& world, @@ -1127,6 +1167,17 @@ namespace madness { return r; } + + /// Returns a deep copy of a vector of functions + template + std::vector< Function > + copy(const std::vector< Function >& v, bool fence=true) { + PROFILE_BLOCK(Vcopy); + std::vector< Function > r(v.size()); + if (v.size()>0) r=copy(v.front().world(),v,fence); + return r; + } + /// Returns a vector of `n` deep copies of a function template std::vector< Function > @@ -1344,7 +1395,6 @@ namespace madness { std::vector< Function > result(f.size()); for (unsigned int i=0; iis_slaterf12); result[i] = apply_only(*op[i], f[i], false); } @@ -1358,6 +1408,15 @@ namespace madness { } + /// Applies an operator to a vector of functions --- q[i] = apply(op,f[i]) + template + std::vector< Function > + apply(const SeparatedConvolution& op, + const std::vector< Function > f) { + return apply(op.get_world(),op,f); + } + + /// Applies an operator to a vector of functions --- q[i] = apply(op,f[i]) template std::vector< Function > @@ -1367,7 +1426,7 @@ namespace madness { PROFILE_BLOCK(Vapply); std::vector< Function >& ncf = *const_cast< std::vector< Function >* >(&f); - bool print_timings=(NDIM==6) and (world.rank()==0); + bool print_timings=(NDIM==6) and (world.rank()==0) and op.print_timings; double wall0=wall_time(); reconstruct(world, f); @@ -1397,15 +1456,6 @@ namespace madness { } reconstruct(world, result); - if (op.is_slaterf12) { - MADNESS_ASSERT(not op.destructive()); - if (typeid(T)!=typeid(R)) MADNESS_EXCEPTION("think again!",1); - for (unsigned int i=0; i full_tensor_copy() const {return copy(*this);} + Tensor full_tensor_copy() {return copy(*this);} bool is_assigned() const {return this->size()>0;}; bool has_data() const {return this->size()>0;}; diff --git a/src/madness/tensor/lowranktensor.h b/src/madness/tensor/lowranktensor.h index 1618038332d..d56be254275 100644 --- a/src/madness/tensor/lowranktensor.h +++ b/src/madness/tensor/lowranktensor.h @@ -695,9 +695,6 @@ class GenTensor { friend GenTensor transform_dir( const GenTensor& t, const Tensor& c, const int axis); - template - friend GenTensor outer( - const GenTensor& t1, const GenTensor& t2); std::string what_am_i() const { TensorType tt; @@ -815,7 +812,7 @@ void change_tensor_type(GenTensor& t, const TensorArgs& targs) { /// all other combinations are currently invalid. template GenTensor outer(const GenTensor& t1, - const GenTensor& t2, const TensorArgs final_tensor_args) { + const GenTensor& t2, const TensorArgs final_tensor_args=TensorArgs(-1.0,TT_2D)) { typedef TENSOR_RESULT_TYPE(T,Q) resultT; diff --git a/src/madness/tensor/srconf.h b/src/madness/tensor/srconf.h index 480b0e7f2c5..9ecb4858b48 100644 --- a/src/madness/tensor/srconf.h +++ b/src/madness/tensor/srconf.h @@ -693,10 +693,17 @@ namespace madness { Tensor make_vector_with_weights(const int dim) const { Tensor v=copy(vector_[dim].reshape(rank(),vector_[dim].size()/rank())); for (unsigned int r=0; r flat_vector_with_weights(const int dim) const { + return make_vector_with_weights(dim).reshape(rank(),vector_[dim].size()/rank()); + } + protected: /// return the number of coefficients unsigned int nCoeff() const { diff --git a/src/madness/tensor/tensor.h b/src/madness/tensor/tensor.h index c438f2dcc07..bdb0a77dfad 100644 --- a/src/madness/tensor/tensor.h +++ b/src/madness/tensor/tensor.h @@ -260,6 +260,9 @@ namespace madness { template T mynorm(T t) { return t*t; } + template double mynorm(int t) { + return double(t)*double(t); + } template T mynorm(std::complex t) { return std::norm(t); diff --git a/src/madness/tensor/tensor_lapack.h b/src/madness/tensor/tensor_lapack.h index f446236a5b0..1215aebbb9a 100644 --- a/src/madness/tensor/tensor_lapack.h +++ b/src/madness/tensor/tensor_lapack.h @@ -60,6 +60,20 @@ namespace madness { void svd_result(Tensor& a, Tensor& U, Tensor< typename Tensor::scalar_type >& s, Tensor& VT, Tensor& work); + /// SVD - MATLAB syntax + + /// call as + /// auto [U,s,VT] = svd(A); + /// with A=U*S*VT + template + std::tuple, Tensor< typename Tensor::scalar_type >, Tensor> + svd(const Tensor& A) { + Tensor U,VT; + Tensor< typename Tensor::scalar_type > s; + svd(A,U,s,VT); + return std::make_tuple(U,s,VT); + } + /// Solves linear equations /// \ingroup linalg diff --git a/src/madness/world/test_utilities.h b/src/madness/world/test_utilities.h index 345b02f262c..d87cfe9a1d4 100644 --- a/src/madness/world/test_utilities.h +++ b/src/madness/world/test_utilities.h @@ -16,6 +16,8 @@ struct test_output { test_output(std::string line) { std::cout << ltrim_to_length(line,70); logger << std::scientific << std::setprecision(8) ; + time_begin=cpu_time(); + time_last_checkpoint=time_begin; set_cout_to_logger(); } @@ -42,7 +44,9 @@ struct test_output { if (not have_checkpoints) print(""); // first checkpoint have_checkpoints=true; std::cout << " " << ltrim_to_length(message,66); - print_success_fail(std::cout,success,time); + double time1=cpu_time()-time_last_checkpoint; + time_last_checkpoint=cpu_time(); + print_success_fail(std::cout,success,time1); if (not success) { print_and_clear_log(); } @@ -65,7 +69,8 @@ struct test_output { set_cout_to_terminal(); if (have_checkpoints) std::cout << ltrim_to_length("--> final result -->",70); success = success and final_success; - print_success_fail(std::cout,success); + double time_end=cpu_time(); + print_success_fail(std::cout,success,time_end-time_begin); if (not success) print_and_clear_log(); return (success) ? 0 : 1; } @@ -94,6 +99,8 @@ struct test_output { bool cout_set_to_logger=false; // do not change this directly! bool have_checkpoints=false; std::streambuf* stream_buffer_cout; + double time_begin=0.0; + double time_last_checkpoint=0.0; }; diff --git a/src/madness/world/timing_utilities.h b/src/madness/world/timing_utilities.h index 8b2a80266fd..35b443b12ed 100644 --- a/src/madness/world/timing_utilities.h +++ b/src/madness/world/timing_utilities.h @@ -17,7 +17,7 @@ struct timer { sss = cpu_time(); } - void tag(const std::string msg) { + double tag(const std::string msg) { world.gop.fence(); double tt1 = wall_time() - ttt; double ss1 = cpu_time() - sss; @@ -29,10 +29,11 @@ struct timer { } ttt = wall_time(); sss = cpu_time(); + return ss1; } - void end(const std::string msg) { - tag(msg); + double end(const std::string msg) { + return tag(msg); // world.gop.fence(); // double tt1 = wall_time() - ttt; // double ss1 = cpu_time() - sss;