Merge pull request #498 from m-a-d-n-e-s-s/pr-cc-refactoring

major merge including functionality for MP3, ccpairfunctions, low-rank functions and various other stuff.

src/apps/cc2/cc2.cc

@@ -93,13 +93,10 @@ int main(int argc, char **argv) {
             calc->param.print("dft","end");
             print("\n");
             cc2.tdhf->get_parameters().print("response","end");
+            print("\n");
+            nemo->molecule().print();
         }
         double hf_energy = nemo->value();
-        if (world.rank() == 0)
-            std::cout << "\n\n\n\n\n\n Reference Calculation Ended\n SCF Energy is: " << hf_energy
-                      << "\n current wall-time: " << wall_time()
-                      << "\n current cpu-time: " << cpu_time() << "\n\n\n";
-
         cc2.solve();
 
         if (world.rank() == 0) printf("\nfinished at time %.1fs\n\n", wall_time());

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();
+
+            double mp3_energy=mp3_correction+mp2_energy;
+            if(world.rank() == 0) {
+                printf("final hf/mp2/mp3/total energy %12.8f %12.8f %12.8f %12.8f\n",
+                       hf_energy,mp2_energy,mp3_correction,hf_energy+mp3_energy);
             }
         } catch (std::exception& e) {
 

src/examples/heat.cc

@@ -130,7 +130,9 @@ int main(int argc, char** argv) {
     expnt[0] = 1.0/(4.0*c*tstep);
     coeff[0] = pow(4.0*constants::pi*c*tstep,-1.5);
 
-    operatorT G(world, coeff, expnt);
+    double lo_dummy=1.e-4;
+    double thresh_dummy=1.e-6;
+    operatorT G(world, coeff, expnt, lo_dummy, thresh_dummy);
 
     functionT ut = G(u0);
 

src/examples/heat2.cc

@@ -147,7 +147,7 @@ int main(int argc, char** argv) {
     real_tensor expnt(1), coeff(1);
     expnt[0] = 1.0/(4.0*c*tstep*0.5);
     coeff[0] = pow(4.0*constants::pi*c*tstep*0.5,-1.5);
-    real_convolution_3d G(world, coeff, expnt);
+    real_convolution_3d G(world, coeff, expnt, 1.e-4, 1.e-4);
 
     // Propagate forward 50 time steps
     real_function_3d u = u0;

src/examples/smooth.h

@@ -607,7 +607,7 @@ class smooth {
         Tensor<double> coeffs(1), exponents(1);
         exponents(0L) = 1.0 / (2.0 * 0.04);
         coeffs(0L) = pow((1.0 / (2.0 * 0.04)) / M_PI, 0.5);
-        SeparatedConvolution<double, 1> op(world, coeffs, exponents);
+        SeparatedConvolution<double, 1> op(world, coeffs, exponents,1.e-5,1.e-5);
         real_function_1d sf = apply(op, f);
         double diff = (f - sf).norm2();
         output("||f - smoothed_f||_1D=" + stringify(diff));

src/examples/tdse_example.cc

@@ -175,7 +175,7 @@ void converge(World& world, functionT& potn, functionT& psi, double& eps) {
 
     Tensor<double> coeff(1); coeff[0] = 1.0/pow(constants::pi*tmax,1.5);
     Tensor<double> expnt(1); expnt[0] = 1.0/tmax;
-    operatorT* op = new operatorT(world,coeff,expnt);
+    operatorT* op = new operatorT(world,coeff,expnt,0.0,0.0);
 
     for (int iter=0; iter<20; iter++) {
         if (world.rank() == 0) print("ITER",iter);
@@ -197,7 +197,7 @@ void converge(World& world, functionT& potn, functionT& psi, double& eps) {
             delete op;
             coeff[0] = 1.0/sqrt(constants::pi*tmax);
             expnt[0] = 1.0/tmax;
-            op = new operatorT(world,coeff,expnt);
+            op = new operatorT(world,coeff,expnt,0.0,0.0);
         }
     }
     delete op;

src/madness/chem/CC2.h

@@ -12,6 +12,7 @@
 #include<madness/chem/SCF.h>
 #include<madness/chem/nemo.h>
 #include<madness/chem/CCPotentials.h>
+#include<madness/chem/mp3.h>
 #include <madness/mra/operator.h>
 #include <madness/mra/mra.h>
 #include <madness/mra/vmra.h>
@@ -131,12 +132,15 @@ class CC2 : public OptimizationTargetInterface, public QCPropertyInterface {
     /// solve the CC2 ground state equations, returns the correlation energy
     void solve();
 
+
     std::vector<CC_vecfunction>
     solve_ccs();
 
-    /// solve the MP2 equations (uncoupled -> Canonical Orbitals)
-    double
-    solve_mp2(Pairs<CCPair>& doubles);
+    double compute_mp3(const Pairs<CCPair>& mp2pairs) const {
+        MP3 mp3(CCOPS);
+        double mp3_contribution=mp3.mp3_energy_contribution_macrotask_driver(mp2pairs);
+        return mp3_contribution;
+    }
 
     double
     solve_cc2(CC_vecfunction& tau, Pairs<CCPair>& u);
@@ -404,28 +408,6 @@ class CC2 : public OptimizationTargetInterface, public QCPropertyInterface {
     iterate_lrcc2_pairs(const CC_vecfunction& cc2_s, const Pairs<CCPair>& cc2_d, const CC_vecfunction lrcc2_s,
                         Pairs<CCPair>& lrcc2_d);
 
-    bool update_constant_part_mp2(CCPair& pair) {
-        MADNESS_ASSERT(pair.ctype == CT_MP2);
-        MADNESS_ASSERT(pair.type == GROUND_STATE);
-        if (parameters.no_compute_mp2_constantpart()) {
-            pair.constant_part=real_factory_6d(world);
-            load(pair.constant_part,pair.name()+"_const");
-        }
-        if (pair.constant_part.is_initialized()) return false;
-
-        // make screening Operator
-        real_convolution_6d Gscreen = BSHOperator<6>(world, sqrt(-2.0 * CCOPS.get_epsilon(pair.i, pair.j)),
-                                                     parameters.lo(), parameters.thresh_bsh_6D());
-        Gscreen.modified() = true;
-
-        const CCFunction& moi = CCOPS.mo_ket(pair.i);
-        const CCFunction& moj = CCOPS.mo_ket(pair.j);
-
-        pair.constant_part = CCOPS.make_constant_part_mp2(moi, moj, &Gscreen);
-        save(pair.constant_part, pair.name() + "_const");
-        return true;
-    }
-
     bool update_constant_part_cc2_gs(const CC_vecfunction& tau, CCPair& pair) {
         MADNESS_ASSERT(pair.ctype == CT_CC2);
         MADNESS_ASSERT(pair.type == GROUND_STATE);