Merge branch 'master' into pr-pno-update

src/apps/chem/TDHF.cc

@@ -40,7 +40,13 @@ struct TDHF_allocator {
 
 /// \param world    the world
 /// \param parser   the parser
-TDHF::TDHF(World &world, const commandlineparser &parser)
+TDHF::TDHF(World &world, const commandlineparser &parser) : TDHF(world,parser,std::make_shared<Nemo>(world,parser)) {}
+
+///  ctor with command line parser, constructs SCF and Nemo objects on-the-fly, and delegates further
+
+/// \param world    the world
+/// \param parser   the parser
+TDHF::TDHF(World &world, const commandlineparser &parser, std::shared_ptr<Nemo> nemo)
         : world(world),
           reference_(std::make_shared<Nemo>(world, parser)),
           parameters(world, parser.value("input")),
@@ -50,7 +56,6 @@ TDHF::TDHF(World &world, const commandlineparser &parser)
           Q(),
           msg(world) {
 
-
     if (parameters.do_oep()) {
         std::shared_ptr<OEP> oep(new OEP(world, parser));
         set_reference(oep);

src/apps/chem/TDHF.h

@@ -208,6 +208,8 @@ class TDHF {
 
     TDHF(World &world, const commandlineparser &parser);
 
+    TDHF(World &world, const commandlineparser &parser, std::shared_ptr<Nemo> nemo);
+
     void initialize();
 
     ///  sets the reference wave function (nemo or oep)

src/apps/chem/mp2.cc

@@ -75,9 +75,8 @@ void load_balance(const real_function_6d& f, const bool leaf) {
 
 /// ctor
 MP2::MP2(World& world, const commandlineparser& parser) : world(world),
-                                                          param(world), corrfac(world), correlation_energy(0.0),
-                                                          coords_sum(-1.0), Q12(world), ttt(0.0), sss(0.0) {
-
+                   param(world, parser), corrfac(world), correlation_energy(0.0),
+                   coords_sum(-1.0), Q12(world), ttt(0.0), sss(0.0) {
     {
         std::shared_ptr<Nemo> nemo = std::shared_ptr<Nemo>(new Nemo(world, parser));
 

src/apps/chem/mp2.h

@@ -318,7 +318,7 @@ class MP2 : public OptimizationTargetInterface {
         bool do_oep1 = false;
 
         /// ctor reading out the input file
-        Parameters(World& world) {
+        Parameters(World& world, const commandlineparser& parser) {
 
             /// the map with initial values
             initialize < double > ("thresh", 1.e-3, "recommended values: 1.e-4 < econv < 1.e-8");
@@ -330,14 +330,15 @@ class MP2 : public OptimizationTargetInterface {
             initialize < bool > ("restart", true);
             initialize < int > ("maxiter", 5);
 
-            read_and_set_derived_values(world);
+            read_and_set_derived_values(world,parser);
 
             // print final parameters
             if (world.rank() == 0) print("mp2", "end");
         }
 
-        void read_and_set_derived_values(World& world) {
-            read(world, "input", "mp2");
+        void read_and_set_derived_values(World& world, const commandlineparser& parser) {
+            read(world, parser.value("input"), "mp2");
+
             set_derived_value("dconv", sqrt(get<double>("econv")) * 0.1);
         	set_derived_value("thresh",get<double>("econv"));
         }

src/apps/chem/nemo.cc

@@ -289,13 +289,13 @@ std::shared_ptr<Fock<double,3>> Nemo::make_fock_operator() const {
 /// compute the Fock matrix from scratch
 tensorT Nemo::compute_fock_matrix(const vecfuncT& nemo, const tensorT& occ) const {
 	// apply all potentials (J, K, Vnuc) on the nemos
-	vecfuncT psi, Jnemo, Knemo, pcmnemo, JKVpsi, Unemo;
+	vecfuncT Jnemo, Knemo, pcmnemo, JKVpsi, Unemo;
 
     vecfuncT R2nemo=mul(world,R_square,nemo);
     truncate(world,R2nemo);
 
     // compute potentials the Fock matrix: J - K + Vnuc
-	compute_nemo_potentials(nemo, psi, Jnemo, Knemo, pcmnemo, Unemo);
+	compute_nemo_potentials(nemo, Jnemo, Knemo, pcmnemo, Unemo);
 
 //    vecfuncT JKUpsi=add(world, sub(world, Jnemo, Knemo), Unemo);
     vecfuncT JKUpsi=Unemo+Jnemo-Knemo;
@@ -319,7 +319,7 @@ double Nemo::solve(const SCFProtocol& proto) {
 	// Therefore set all tolerance thresholds to zero, also in the mul_sparse
 
 	// apply all potentials (J, K, Vnuc) on the nemos
-	vecfuncT psi, Jnemo, Knemo, pcmnemo, Unemo;
+	vecfuncT Jnemo, Knemo, pcmnemo, Unemo;
 
 	std::vector<double> energies(1,0.0);	// contains the total energy and all its contributions
 	double energy=0.0;
@@ -344,7 +344,7 @@ double Nemo::solve(const SCFProtocol& proto) {
 	    truncate(world,R2nemo);
 
 		// compute potentials the Fock matrix: J - K + Vnuc
-		compute_nemo_potentials(nemo, psi, Jnemo, Knemo, pcmnemo, Unemo);
+		compute_nemo_potentials(nemo, Jnemo, Knemo, pcmnemo, Unemo);
 
 		// compute the energy
 		std::vector<double> oldenergies=energies;
@@ -371,7 +371,7 @@ double Nemo::solve(const SCFProtocol& proto) {
             if (world.rank()==0) print("F max off-diagonal  ",max_fock_offidag);
 
             // canonicalize the orbitals, rotate subspace and potentials
-        	tensorT overlap = matrix_inner(world, psi, psi, true);
+            tensorT overlap = matrix_inner(world, R2nemo, nemo, true);
             tensorT U=calc->get_fock_transformation(world,overlap,
                     fock,calc->aeps,calc->aocc,FunctionDefaults<3>::get_thresh());
 
@@ -566,7 +566,7 @@ std::vector<double> Nemo::compute_energy_regularized(const vecfuncT& nemo, const
 /// @param[out]	Knemo	exchange operator applied on the nemos
 /// @param[out]	Vnemo	nuclear potential applied on the nemos
 /// @param[out]	Unemo	regularized nuclear potential applied on the nemos
-void Nemo::compute_nemo_potentials(const vecfuncT& nemo, vecfuncT& psi,
+void Nemo::compute_nemo_potentials(const vecfuncT& nemo,
 		vecfuncT& Jnemo, vecfuncT& Knemo, vecfuncT& pcmnemo,
 		vecfuncT& Unemo) const {
 

src/apps/chem/nemo.h

@@ -629,12 +629,11 @@ class Nemo: public NemoBase {
 	/// to use these potentials in the fock matrix computation they must
 	/// be multiplied by the nuclear correlation factor
 	/// @param[in]	nemo	the nemo orbitals
-	/// @param[out]	psi		the reconstructed, full orbitals
 	/// @param[out]	Jnemo	Coulomb operator applied on the nemos
 	/// @param[out]	Knemo	exchange operator applied on the nemos
 	/// @param[out]	pcmnemo	PCM (solvent) potential applied on the nemos
 	/// @param[out]	Unemo	regularized nuclear potential applied on the nemos
-	void compute_nemo_potentials(const vecfuncT& nemo, vecfuncT& psi,
+	void compute_nemo_potentials(const vecfuncT& nemo,
 			vecfuncT& Jnemo, vecfuncT& Knemo, vecfuncT& pcmnemo,
 			vecfuncT& Unemo) const;
 

src/examples/cc2.cc

@@ -55,73 +55,72 @@ using namespace madness;
 
 #ifdef USE_GENTENSOR
 
-int main(int argc, char** argv) {
+int main(int argc, char **argv) {
 
-	initialize(argc, argv);
+    initialize(argc, argv);
 
-	World world(SafeMPI::COMM_WORLD);
+    World world(SafeMPI::COMM_WORLD);
 
-	if (world.rank() == 0) {
-		print("\n  CC2: Coupled Cluster approximate Doubles  \n");
-		printf("starting at time %.1f\n", wall_time());
-		print("\nmain() compiled at ",__TIME__," on ",__DATE__);
+    if (world.rank() == 0) {
+        print("\n  CC2: Coupled Cluster approximate Doubles  \n");
+        printf("starting at time %.1f\n", wall_time());
+        print("\nmain() compiled at ", __TIME__, " on ", __DATE__);
 
-	}
-	startup(world,argc,argv);
-	std::cout.precision(6);
-	FunctionDefaults<3>::set_truncate_mode(1);
-	if(world.rank()==0) print("Truncate mode set to ",FunctionDefaults<3>::get_truncate_mode());
+    }
+    startup(world, argc, argv);
+    std::cout.precision(6);
+    FunctionDefaults<3>::set_truncate_mode(1);
+    if (world.rank() == 0) print("Truncate mode set to ", FunctionDefaults<3>::get_truncate_mode());
 
 #ifdef GITREVISION
-	const  char* gitrev =  GITREVISION;
-	const std::string gitrevision(gitrev);
-	if (world.rank()==0) {
-		print("           git revision ...",gitrevision);
-	}
+    const  char* gitrev =  GITREVISION;
+    const std::string gitrevision(gitrev);
+    if (world.rank()==0) {
+        print("           git revision ...",gitrevision);
+    }
 #endif
 
-typedef std::vector<real_function_3d> vecfuncT;
-
 // set the tensor type
-TensorType tt=TT_2D;
-FunctionDefaults<6>::set_tensor_type(tt);
-FunctionDefaults<6>::set_apply_randomize(true);
+    TensorType tt = TT_2D;
+    FunctionDefaults<6>::set_tensor_type(tt);
+    FunctionDefaults<6>::set_apply_randomize(true);
 
 // Process 0 reads input information and broadcasts
-const char * inpname = "input";
-for (int i=1; i<argc; i++) {
-    if (argv[i][0] != '-') {
-        inpname = argv[i];
-        break;
+    const char *inpname = "input";
+    for (int i = 1; i < argc; i++) {
+        if (argv[i][0] != '-') {
+            inpname = argv[i];
+            break;
+        }
     }
-}
-if (world.rank() == 0) print("input filename: ", inpname);
+    if (world.rank() == 0) print("input filename: ", inpname);
 
 //SCF calc(world,input.c_str());
-commandlineparser parser(argc,argv);
-std::shared_ptr<Nemo> nemo(new Nemo(world,parser));
-std::shared_ptr<SCF> calc=nemo->get_calc();
-if (world.rank()==0) {
-    calc->molecule.print();
-    print("\n");
-    calc->param.print("cc2");
-}
-double hf_energy = nemo->value();
-if(world.rank()==0) std::cout << "\n\n\n\n\n\n Reference Calclation Ended\n SCF Energy is: " << hf_energy
-		<<"\n current wall-time: " << wall_time()
-		<<"\n current cpu-time: " << cpu_time()<< "\n\n\n";
+    commandlineparser parser(argc, argv);
+    std::shared_ptr<Nemo> nemo(new Nemo(world, parser));
+    std::shared_ptr<SCF> calc = nemo->get_calc();
+    if (world.rank() == 0) {
+        calc->molecule.print();
+        print("\n");
+        calc->param.print("cc2");
+    }
+    double hf_energy = nemo->value();
+    if (world.rank() == 0)
+        std::cout << "\n\n\n\n\n\n Reference Calclation Ended\n SCF Energy is: " << hf_energy
+                  << "\n current wall-time: " << wall_time()
+                  << "\n current cpu-time: " << cpu_time() << "\n\n\n";
 
 
 // Make CC2
-CC2 cc2(world,inpname,nemo);
+    CC2 cc2(world, parser, nemo);
 
-cc2.solve();
+    cc2.solve();
 
-if(world.rank() == 0) printf("\nfinished at time %.1fs\n\n", wall_time());
-world.gop.fence();
-finalize();
+    if (world.rank() == 0) printf("\nfinished at time %.1fs\n\n", wall_time());
+    world.gop.fence();
+    finalize();
 
-	return 0;
+    return 0;
 }// end main
 
 #else
@@ -131,9 +130,9 @@ int main(int argc, char** argv) {
     startup(world,argc,argv);
     if(world.rank() == 0) {
 
-    	print("\nYou can't run cc2 because you have configured MADNESS ");
-    	print("without the --enable-gentensor flag");
-    	print("You need to reconfigure and recompile\n");
+        print("\nYou can't run cc2 because you have configured MADNESS ");
+        print("without the --enable-gentensor flag");
+        print("You need to reconfigure and recompile\n");
 
     }
     world.gop.fence();

src/madness/mra/funcimpl.h

@@ -2027,7 +2027,7 @@ namespace madness {
             world.gop.fence(); // fence is necessary if trickle down is used afterwards
             // trickle down and undo redundand shouldnt change anything if only the diagonal elements are considered above -> check this
             f->trickle_down(true); // fence must be true otherwise undo_redundant will have trouble
-            f->undo_redundant(true);
+//            f->undo_redundant(true);
             f->verify_tree();
             //if (fence) world.gop.fence(); // unnecessary, fence is activated in undo_redundant
 

src/madness/tensor/RandomizedMatrixDecomposition.h

@@ -14,6 +14,8 @@ struct RMDFactory {
 	long maxrank_=LONG_MAX;
 	long oversampling_=10;
 
+	RMDFactory() {}
+
 	RMDFactory& maxrank(const long mr) {
 		maxrank_=mr;
 		return *this;