Merge pull request #299 from m-a-d-n-e-s-s/TDHF-Cleanup

TDHF Cleanup
m-a-d-n-e-s-s · Oct 12, 2018 · 397af84 · 397af84
2 parents 69b168c + 31fc4c6
commit 397af84
Show file tree

Hide file tree

Showing 19 changed files with 2,160 additions and 1,021 deletions.
diff --git a/src/apps/chem/CC2.cc b/src/apps/chem/CC2.cc
@@ -17,7 +17,7 @@ namespace madness {
     const CalcType ctype = parameters.calculation;
 
     if(ctype==CT_TDHF){
-      TDHF tdhf(world,parameters,nemo);
+      TDHF tdhf(world,nemo);
       std::vector<CC_vecfunction> ccs;
       for(size_t k=0;k<parameters.excitations_.size();k++){
 	CC_vecfunction tmp;
@@ -65,7 +65,7 @@ namespace madness {
       if(world.rank()==0) std::cout << std::fixed << std::setprecision(10) << " CC2 Correlation Energy =" << cc2_correlation_energy << "\n";
 
     }else if(ctype==CT_LRCCS){
-      TDHF tdhf(world,parameters,nemo);
+      TDHF tdhf(world,nemo);
       std::vector<CC_vecfunction> ccs;
       for(size_t k=0;k<parameters.excitations_.size();k++){
 	CC_vecfunction tmp;
@@ -302,7 +302,7 @@ namespace madness {
 	CCTimer time_ex(world,"LRCC2 Calculation for Excitation " + std::to_string(int(excitation)));
 	CC_vecfunction lrcc2_s = vccs[xxx];
 	// needed to assign an omega
-	const vecfuncT backup = copy(world,lrcc2_s.get_vecfunction());
+	const vector_real_function_3d backup = copy(world,lrcc2_s.get_vecfunction());
 	CC_vecfunction test(backup,RESPONSE,parameters.freeze);
 	test.excitation = lrcc2_s.excitation;
 	iterate_ccs_singles(test);
@@ -360,7 +360,7 @@ namespace madness {
   // Solve the CCS equations for the ground state (debug potential and check HF convergence)
   std::vector<CC_vecfunction> CC2::solve_ccs() {
     output.section("SOLVE CCS");
-    TDHF tdhf(world,parameters,nemo);
+    TDHF tdhf(world,nemo);
     std::vector<CC_vecfunction> excitations;
     for(size_t k=0;k<parameters.excitations_.size();k++){
 	CC_vecfunction tmp;

diff --git a/src/apps/chem/CC2.h b/src/apps/chem/CC2.h
@@ -25,6 +25,16 @@ namespace madness {
   class CC2{
   public:
 
+	  CC2(World &world_,const CCParameters& param,const Nemo &nemo_)
+	 	: world(world_),
+	 	parameters(param),
+	 	nemo(nemo_),
+	 	CCOPS(world,nemo,parameters),
+	 	output(CCOPS.output)
+	   {
+		  parameters.sanity_check(world);
+	   }
+
     CC2(World &world_,const std::string &inputFileName,const Nemo &nemo_)
 	: world(world_),
 	parameters(inputFileName,nemo_.get_calc()->param.lo),
@@ -67,8 +77,6 @@ namespace madness {
       else if(nemo.nuclear_correlation->type() == NuclearCorrelationFactor::Two) nuc="Two";
       if(world.rank() == 0) std::cout << "Nuclear Correlation Factor is " << nuc << std::endl;
 
-      //output.section("Testing Section in Constructor");
-      //CCOPS.test_fill_tree();
     }
 
     void
@@ -217,7 +225,7 @@ namespace madness {
 
 	// get potentials
 	CCTimer time_V(world,assign_name(ctype) + "-Singles Potential");
-	vecfuncT V;
+	vector_real_function_3d V;
 	if(ctype == CT_CC2) V=CCOPS.get_CC2_singles_potential_gs(singles,gs_doubles);
 	else if(ctype == CT_LRCC2) V=CCOPS.get_CC2_singles_potential_ex(singles2,gs_doubles,singles,ex_doubles);
 	else if(ctype == CT_LRCCS) V=CCOPS.get_CCS_potential_ex(singles);
@@ -248,7 +256,7 @@ namespace madness {
 
 	// apply bsh operators
 	CCTimer time_applyG(world,"Apply G-Operators");
-	vecfuncT GV=apply<SeparatedConvolution<double, 3>, double, 3>(world,G,V);
+	vector_real_function_3d GV=apply<SeparatedConvolution<double, 3>, double, 3>(world,G,V);
 	world.gop.fence();
 	time_applyG.info();
 
@@ -258,15 +266,15 @@ namespace madness {
 	// Normalize Singles if it is excited state
 	if(ctype==CT_LRCCS or ctype==CT_LRCC2 or ctype==CT_ADC2){
 	  output("Normalizing new singles");
-	  const vecfuncT x = GV;
-	  const vecfuncT xbra = mul(world,nemo.nuclear_correlation->square(),GV);
+	  const vector_real_function_3d x = GV;
+	  const vector_real_function_3d xbra = mul(world,nemo.nuclear_correlation->square(),GV);
 	  const double norm = sqrt(inner(world,xbra,x).sum());
 	  if(world.rank()==0) std::cout << " Norm was " <<std::fixed<< std::setprecision(parameters.output_prec) << norm << "\n";
 	  scale(world,GV,1.0/norm);
 	} else output("Singles not normalized");
 
 	// residual
-	const vecfuncT residual=sub(world,singles.get_vecfunction(),GV);
+	const vector_real_function_3d residual=sub(world,singles.get_vecfunction(),GV);
 
 	// information
 	const Tensor<double> R2xinnerx=inner(world,mul(world,nemo.nuclear_correlation->square(),singles.get_vecfunction()),singles.get_vecfunction());
@@ -292,8 +300,8 @@ namespace madness {
 	  output("\nMake 2nd order energy update:");
 	  // include nuclear factors
 	  {
-	    vecfuncT bra_res=mul(world,nemo.nuclear_correlation->square(),residual);
-	    vecfuncT bra_GV=mul(world,nemo.nuclear_correlation->square(),GV);
+	    vector_real_function_3d bra_res=mul(world,nemo.nuclear_correlation->square(),residual);
+	    vector_real_function_3d bra_GV=mul(world,nemo.nuclear_correlation->square(),GV);
 	    double Rtmp=inner(world,bra_res,V).sum();
 	    double Rtmp2=inner(world,bra_GV,GV).sum();
 	    const double Rdelta=(0.5 * Rtmp / Rtmp2);
@@ -307,7 +315,7 @@ namespace madness {
 
 	// update singles
 	singles.omega=omega;
-	vecfuncT new_singles=GV;
+	vector_real_function_3d new_singles=GV;
 	if(parameters.kain) new_singles=solver.update(singles.get_vecfunction(),residual).x;
 	print_size(world,new_singles,"new_singles");
 	truncate(world,new_singles);