local orbitals in lrcc2

src/madness/chem/CC2.cc

@@ -474,7 +474,7 @@ double CC2::solve_mp2_coupled(Pairs<CCPair>& doubles, Info& info) {
         if (world.rank()==0) print_header3("Starting iteration " + std::to_string(int(iter)) + " of MP2");
 
         // compute the coupling between the pair functions
-        Pairs<real_function_6d> coupling=compute_local_coupling(pair_vec);
+        Pairs<real_function_6d> coupling=compute_local_coupling(pair_vec, info);
         auto coupling_vec=Pairs<real_function_6d>::pairs2vector(coupling,triangular_map);
         if (parameters.debug()) print_size(world, coupling_vec, "couplingvector");
 
@@ -556,15 +556,16 @@ double CC2::solve_mp2_coupled(Pairs<CCPair>& doubles, Info& info) {
 /// add the coupling terms for local MP2
 
 /// @return \sum_{k\neq i} f_ki |u_kj> + \sum_{l\neq j} f_lj |u_il>
-Pairs<real_function_6d> CC2::compute_local_coupling(const Pairs<real_function_6d>& pairs) const {
+Pairs<real_function_6d> CC2::compute_local_coupling(const Pairs<real_function_6d>& pairs, const Info& info) {
 
-    const int nmo = nemo->get_calc()->amo.size();
+    const int nmo = info.mo_ket.size();
+    World& world=pairs.allpairs.begin()->second.world();
 
     // temporarily make all N^2 pair functions
     typedef std::map<std::pair<int, int>, real_function_6d> pairsT;
     pairsT quadratic;
-    for (int k = parameters.freeze(); k < nmo; ++k) {
-        for (int l = parameters.freeze(); l < nmo; ++l) {
+    for (int k = info.parameters.freeze(); k < nmo; ++k) {
+        for (int l = info.parameters.freeze(); l < nmo; ++l) {
             if (l >= k) {
                 quadratic[std::make_pair(k, l)] = pairs(k, l);
             } else {
@@ -577,28 +578,29 @@ Pairs<real_function_6d> CC2::compute_local_coupling(const Pairs<real_function_6d
     world.gop.fence();
 
     // the coupling matrix is the Fock matrix, skipping diagonal elements
-    Tensor<double> fock1 = nemo->compute_fock_matrix(nemo->get_calc()->amo, nemo->get_calc()->aocc);
+    // Tensor<double> fock1 = nemo->compute_fock_matrix(nemo->get_calc()->amo, nemo->get_calc()->aocc);
+    Tensor<double> fock1 = copy(info.fock);
     for (int k = 0; k < nmo; ++k) {
         if (fock1(k, k) > 0.0) MADNESS_EXCEPTION("positive orbital energies", 1);
         fock1(k, k) = 0.0;
     }
 
     Pairs<real_function_6d> coupling;
-    for (int i = parameters.freeze(); i < nmo; ++i) {
+    for (int i = info.parameters.freeze(); i < nmo; ++i) {
         for (int j = i; j < nmo; ++j) {
             coupling.insert(i, j, real_factory_6d(world).compressed());
         }
     }
 
-    for (int i = parameters.freeze(); i < nmo; ++i) {
+    for (int i = info.parameters.freeze(); i < nmo; ++i) {
         for (int j = i; j < nmo; ++j) {
-            for (int k = parameters.freeze(); k < nmo; ++k) {
+            for (int k = info.parameters.freeze(); k < nmo; ++k) {
                 if (fock1(k, i) != 0.0) {
                     coupling(i, j).gaxpy(1.0, quadratic[std::make_pair(k, j)], fock1(k, i), false);
                 }
             }
 
-            for (int l = parameters.freeze(); l < nmo; ++l) {
+            for (int l = info.parameters.freeze(); l < nmo; ++l) {
                 if (fock1(l, j) != 0.0) {
                     coupling(i, j).gaxpy(1.0, quadratic[std::make_pair(i, l)], fock1(l, j), false);
                 }
@@ -695,14 +697,21 @@ CC2::iterate_lrcc2_pairs(World& world, const CC_vecfunction& cc2_s,
     }
     for (const auto& p : pair_vec) p.function().print_size("u before iter");
 
+            // compute the coupling between the pair functions
+            if (world.rank()==0) print("computing local coupling in the universe");
+            Pairs<real_function_6d> coupling=compute_local_coupling(pair_vec, info);
+            auto coupling_vec=Pairs<real_function_6d>::pairs2vector(coupling,triangular_map);
+        if (info.parameters.debug()) print_size(world, coupling_vec, "couplingvector");
+
+
     // iterate the pair
     MacroTaskIteratePair t1;
     MacroTask task1(world, t1);
     // temporary fix: create dummy functions to that the cloud is not confused
-    real_function_6d tmp=real_factory_6d(world).functor([](const coord_6d& r){return 0.0;});
-    std::vector<real_function_6d> vdummy_6d(pair_vec.size(),tmp);         // dummy vectors
+    // real_function_6d tmp=real_factory_6d(world).functor([](const coord_6d& r){return 0.0;});
+    // std::vector<real_function_6d> vdummy_6d(pair_vec.size(),tmp);         // dummy vectors
     const std::size_t maxiter=10;
-    auto unew = task1(pair_vec, vdummy_6d, cc2_s, lrcc2_s, info, maxiter);
+    auto unew = task1(pair_vec, coupling_vec, cc2_s, lrcc2_s, info, maxiter);
 
     for (const auto& u : unew) u.print_size("u after iter");
     // get some statistics
@@ -770,7 +779,7 @@ CC2::solve_cc2(CC_vecfunction& singles, Pairs<CCPair>& doubles, Info& info) cons
 
             // compute the coupling between the pair functions
             if (world.rank()==0) print("computing local coupling in the universe");
-            Pairs<real_function_6d> coupling=compute_local_coupling(pair_vec);
+            Pairs<real_function_6d> coupling=compute_local_coupling(pair_vec, info);
             auto coupling_vec=Pairs<real_function_6d>::pairs2vector(coupling,triangular_map);
             timer1.tag("computing local coupling");
 
@@ -862,7 +871,7 @@ CC2::solve_lrcc2(Pairs<CCPair>& gs_doubles, const CC_vecfunction& gs_singles, co
     const double omega_cis = ex_singles.omega;
 
     for (size_t iter = 0; iter < parameters.iter_max(); iter++) {
-        print_header2("Macroiteration " + std::to_string(int(iter)) + " of LRCC2 for excitation energy "+std::to_string(ex_singles.omega));
+        if (world.rank()==0) print_header2("Macroiteration " + std::to_string(int(iter)) + " of LRCC2 for excitation energy "+std::to_string(ex_singles.omega));
         update_reg_residues_ex(world, gs_singles, ex_singles, ex_doubles, info);
         bool dconv = iterate_lrcc2_pairs(world, gs_singles, ex_singles, ex_doubles, info);
         bool sconv = iterate_lrcc2_singles(world, gs_singles, gs_doubles, ex_singles, ex_doubles, info);

src/madness/chem/CC2.h

@@ -366,7 +366,8 @@ class CC2 : public OptimizationTargetInterface, public QCPropertyInterface {
             if (ctype==CT_CC2) update_reg_residues_gs(world, singles,gs_doubles, info);
             else if(ctype==CT_LRCC2) update_reg_residues_ex(world, singles2,singles,ex_doubles, info);
 
-            CCPotentials::print_convergence(singles.name(0),rmsresidual,rmsresidual,omega-old_omega,iter);
+            if (world.rank()==0) CCPotentials::print_convergence(singles.name(0),rmsresidual,
+                rmsresidual,omega-old_omega,iter);
             converged = (R2vector_error < info.parameters.dconv_3D());
 
             time.info();
@@ -485,19 +486,19 @@ class CC2 : public OptimizationTargetInterface, public QCPropertyInterface {
     }
 
     /// forward to the other function (converting CCPair to real_function)
-    Pairs<real_function_6d> compute_local_coupling(const std::vector<CCPair> &vpairs) const {
+    static Pairs<real_function_6d> compute_local_coupling(const std::vector<CCPair> &vpairs, const Info& info) {
         // create new pairs structure
         Pairs<CCPair> pairs;
         for (auto& tmp_pair : vpairs) pairs.insert(tmp_pair.i, tmp_pair.j, tmp_pair);
         auto ccpair2function = [](const CCPair& a) {return a.function();};
-        return compute_local_coupling(pairs.convert<real_function_6d>(pairs,ccpair2function));
+        return compute_local_coupling(pairs.convert<real_function_6d>(pairs,ccpair2function), info);
 
     };
 
     /// add the coupling terms for local MP2
 
     /// \sum_{k\neq i} f_ki |u_kj> + \sum_{l\neq j} f_lj |u_il>
-    Pairs<real_function_6d> compute_local_coupling(const Pairs<real_function_6d>& pairs) const;
+    static Pairs<real_function_6d> compute_local_coupling(const Pairs<real_function_6d>& pairs, const Info& info);
 
 
     double solve_mp2_coupled(Pairs<CCPair> &doubles, Info& info);

src/madness/chem/CCPotentials.cc

@@ -2939,7 +2939,7 @@ CCPotentials::get_CC2_singles_potential_ex(World& world, const CC_vecfunction& g
         const double s4a = inner(world, xbra, Vs4a).sum();
         const double s4b = inner(world, xbra, Vs4b).sum();
         const double s4c = inner(world, xbra, Vs4c).sum();
-        std::cout << std::fixed << std::setprecision(10) << "functional response energies:" << "\n<x|ccs>=" << ccs
+        if (world.rank()==0) std::cout << std::fixed << std::setprecision(10) << "functional response energies:" << "\n<x|ccs>=" << ccs
                   << "\n<x|S2b>=" << s2b << "\n<x|S2c>=" << s2c << "\n<x|s4a>=" << s4a << "\n<x|s4b>="
                   << s4b << "\n<x|s4c>=" << s4c << "\n";
         // debug end

src/madness/chem/CCPotentials.h

@@ -40,6 +40,7 @@ class CCPotentials {
         info.U2 = nemo->ncf->U2();
         info.intermediate_potentials = get_potentials;
         info.orbital_energies = orbital_energies_;
+        info.fock=nemo->compute_fock_matrix(nemo->get_calc()->amo, nemo->get_calc()->aocc);
         return info;
     }
 

src/madness/chem/CCStructures.cc

@@ -177,7 +177,7 @@ void CCParameters::set_derived_values() {
     set_derived_value("tight_thresh_3d",thresh_3D()*0.1);
     set_derived_value("thresh_ue",tight_thresh_6D());
     set_derived_value("dconv_6d",3.0*thresh_6D());
-    set_derived_value("dconv_3d",3.0*thresh_3D());
+    set_derived_value("dconv_3d",0.3*thresh_6D());
     set_derived_value("econv",0.1*dconv_6D());
     set_derived_value("econv_pairs",econv());
 

src/madness/chem/CCStructures.h

@@ -230,7 +230,7 @@ struct CCParameters : public QCCalculationParametersBase {
         initialize < double > ("thresh_Ue", thresh_operators, "ue threshold");
         initialize < double > ("econv", thresh, "overal convergence threshold ");
         initialize < double > ("econv_pairs", 0.1*thresh, "convergence threshold for pairs");
-        initialize < double > ("dconv_3d", 0.03*thresh, "convergence for cc singles");
+        initialize < double > ("dconv_3d", 0.3*thresh, "convergence for cc singles");
         initialize < double > ("dconv_6d", 3.0*thresh, "convergence for cc doubles");
         initialize < std::size_t > ("iter_max", 10, "max iterations");
         initialize < std::size_t > ("iter_max_3d", 10, "max iterations for singles");
@@ -1208,6 +1208,7 @@ struct Info {
     std::vector<madness::Vector<double,3>> molecular_coordinates;
     CCParameters parameters;
     std::vector<double> orbital_energies;
+    Tensor<double> fock;
     CCIntermediatePotentials intermediate_potentials;
     Function<double,3> R_square, U2, R;;
     std::vector<Function<double,3>> U1;
@@ -1233,6 +1234,7 @@ struct Info {
         records+=cloud.store(world,mo_ket);
         records+=cloud.store(world,parameters);
         records+=cloud.store(world,orbital_energies);
+        records+=cloud.store(world,fock);
         records+=cloud.store(world,intermediate_potentials);
         records+=cloud.store(world,R_square);
         records+=cloud.store(world,molecular_coordinates);
@@ -1251,6 +1253,7 @@ struct Info {
         mo_ket=cloud.forward_load<std::vector<Function<double,3>>>(world,recordlist);
         parameters=cloud.forward_load<CCParameters>(world,recordlist);
         orbital_energies=cloud.forward_load<std::vector<double>>(world,recordlist);
+        fock=cloud.forward_load<Tensor<double>>(world,recordlist);
         intermediate_potentials=cloud.forward_load<CCIntermediatePotentials>(world,recordlist);
         R_square=cloud.forward_load<Function<double,3>>(world,recordlist);
         molecular_coordinates=cloud.forward_load<std::vector<madness::Vector<double,3>>>(world,recordlist);