Symmetry operations

src/apps/chem/CCStructures.cc

@@ -22,21 +22,21 @@ namespace madness{
   CCMessenger::section(const std::string& msg) const {
     if(world.rank() == 0){
       std::cout << "\n" << std::setw(msg.size() + 10) << std::setfill('*') << "\n";
-      std::setfill(' ');
+      std::cout << std::setfill(' ');
       output(msg);
       std::cout << std::setw(msg.size() + 10) << std::setfill('*') << "\n\n";
-      std::setfill(' ');
+      std::cout << std::setfill(' ');
     }
   }
 
   void
   CCMessenger::subsection(const std::string& msg) const {
     if(world.rank() == 0){
       std::cout << "\n" << std::setw(msg.size() + 5) << std::setfill('-') << "\n";
-      std::setfill(' ');
+      std::cout << std::setfill(' ');
       output(msg);
       std::cout << std::setw(msg.size() + 5) << std::setfill('-') << "\n";
-      std::setfill(' ');
+      std::cout << std::setfill(' ');
     }
   }
 
@@ -55,7 +55,8 @@ namespace madness{
       if(norm != 12345.6789) s_norm=", ||result||=" + std::to_string(norm);
 
       if(world.rank() == 0){
-	std::cout << std::setfill(' ') << std::scientific << std::setprecision(2) << "Timer: " << time_wall << " (Wall), " << time_cpu << " (CPU)" << s_norm << ", (" + operation + ")" << "\n";
+    	  std::cout << std::setfill(' ') << std::scientific << std::setprecision(2)
+    	  << "Timer: " << time_wall << " (Wall), " << time_cpu << " (CPU)" << s_norm << ", (" + operation + ")" << "\n";
       }
     }
   }
@@ -86,12 +87,14 @@ namespace madness{
 
   madness::CC_vecfunction
   CC_vecfunction::copy() const {
-    std::vector<CCFunction> vn;
-    for(auto x : functions){
-      const CCFunction fn(madness::copy(x.second.function),x.second.i,x.second.type);
-      vn.push_back(fn);
-    }
-    return CC_vecfunction(vn,type);
+	  std::vector<CCFunction> vn;
+	  for(auto x : functions){
+		  const CCFunction fn(madness::copy(x.second.function),x.second.i,x.second.type);
+		  vn.push_back(fn);
+	  }
+	  CC_vecfunction result(vn,type);
+	  result.irrep=irrep;
+	  return result;
   }
 
   std::string

src/apps/chem/CCStructures.h

@@ -391,12 +391,18 @@ namespace madness{
 	functions.insert(std::make_pair(freeze+i,tmp));
       }
     }
-    CC_vecfunction(const std::vector<CCFunction> &v,const FuncType type_): type(type_),omega(0.0),excitation(-1),current_error(99.9),delta(0.0){
-      for(auto x:v){
-	functions.insert(std::make_pair(x.i,x));
-      }
+
+    CC_vecfunction(const std::vector<CCFunction> &v, const FuncType type_)
+    	: type(type_),omega(0.0),excitation(-1),current_error(99.9),delta(0.0) {
+    	for(auto x:v) functions.insert(std::make_pair(x.i,x));
+    }
+
+    /// copy ctor (shallow)
+    CC_vecfunction(const CC_vecfunction &other)
+    	: functions(other.functions),type(other.type), omega(other.omega),
+		  excitation(other.excitation),current_error(other.current_error),
+		  delta(other.delta), irrep(other.irrep) {
     }
-    CC_vecfunction(const CC_vecfunction &other) : functions(other.functions),type(other.type), omega(other.omega),excitation(other.excitation),current_error(other.current_error),delta(other.delta) {}
 
     typedef std::map<std::size_t, CCFunction> CC_functionmap;
     CC_functionmap functions;
@@ -411,6 +417,7 @@ namespace madness{
     int excitation; /// the excitation number
     double current_error;
     double delta; // Last difference in Energy
+    std::string irrep="null";	/// excitation irrep (direct product of x function and corresponding orbital)
 
     std::string
     name() const;

src/apps/chem/CMakeLists.txt

@@ -10,12 +10,13 @@ set(MADCHEM_HEADERS
     molecular_optimizer.h projector.h
     SCFOperators.h CCStructures.h CalculationParameters.h
     electronic_correlation_factor.h cheminfo.h vibanal.h molopt.h TDHF.h CC2.h CCPotentials.h
-    pcm.h SCFProtocol.h AC.h)
+    pcm.h SCFProtocol.h AC.h pointgroupoperator.h pointgroupsymmetry.h)
 set(MADCHEM_SOURCES
     correlationfactor.cc molecule.cc molecularbasis.cc vibanal.cc
     corepotential.cc atomutil.cc lda.cc cheminfo.cc
     distpm.cc SCF.cc gth_pseudopotential.cc nemo.cc mp2.cc pcm.cc
-    SCFOperators.cc TDHF.cc GuessFactory.cc CCStructures.cc CC2.cc CCPotentials.cc AC.cc)
+    SCFOperators.cc TDHF.cc GuessFactory.cc CCStructures.cc CC2.cc CCPotentials.cc 
+    AC.cc pointgroupsymmetry.cc)
 if(LIBXC_FOUND)
   list(APPEND MADCHEM_SOURCES xcfunctional_libxc.cc)
 else()
@@ -51,15 +52,17 @@ install(FILES sto-3g sto-6g 6-31g coredata/mcp coredata/mcp2 coredata/mcp_guess
 # Add unit tests    
 if(ENABLE_UNITTESTS)
 
+  SET(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
   # The list of unit test source files
-  set(CHEM_TEST_SOURCES test_SCFOperators.cc test_dft.cc)
+  set(CHEM_TEST_SOURCES test_SCFOperators.cc test_dft.cc test_pointgroupsymmetry.cc)
 
   add_unittests(chem CHEM_TEST_SOURCES "MADchem;MADgtest")
 
   # Create other test executables not included in the unit tests
   set(CHEM_OTHER_TESTS testxc)
   foreach(_test ${CHEM_OTHER_TESTS})  
-    add_executable(${_test} EXCLUDE_FROM_ALL ${_test}.cc)
+#    add_executable(${_test} EXCLUDE_FROM_ALL ${_test}.cc)
+    add_executable(${_test} ${_test}.cc)
     target_link_libraries(${_test} MADchem)
   endforeach()
 

src/apps/chem/CalculationParameters.h

@@ -72,6 +72,7 @@ struct CalculationParameters {
     bool plotcoul;              ///< If true plot the total coulomb potential at convergence
     bool localize;              ///< If true solve for localized orbitals
     bool localize_pm;           ///< If true use PM for localization
+    std::string symmetry;		///< use point group symmetry for all orbitals: default/full/schoenflies
     bool restart;               ///< If true restart from orbitals on disk
     bool restartao;             ///< If true restart from orbitals projected into AO basis (STO3G) on disk
     bool no_compute;            ///< If true use orbitals on disk, set value to computed
@@ -140,7 +141,7 @@ struct CalculationParameters {
     void serialize(Archive& ar) {
         ar & charge & smear & econv & dconv & k & L & maxrotn & nvalpha & nvbeta
         & nopen & maxiter & nio & spin_restricted;
-        ar & plotlo & plothi & plotdens & plotcoul & localize & localize_pm
+        ar & plotlo & plothi & plotdens & plotcoul & localize & localize_pm & symmetry
         & restart & restartao & save & no_compute &no_orient & maxsub & orbitalshift & npt_plot & plot_cell & aobasis;
         ar & nalpha & nbeta & nmo_alpha & nmo_beta & lo;
         ar & core_type & derivatives & conv_only_dens & dipole;
@@ -170,6 +171,7 @@ struct CalculationParameters {
     , plotcoul(false)
     , localize(true)
     , localize_pm(true)
+    , symmetry("default")
     , restart(false)
     , restartao(false)
     , no_compute(false)
@@ -331,6 +333,14 @@ struct CalculationParameters {
             else if (s == "boys") {
                 localize_pm = false;
             }
+            else if (s == "symmetry") {
+            	symmetry="full";
+            	std::string buf,buf1;
+                std::getline(f,buf);
+            	std::stringstream ff(buf);
+            	ff >> buf1;
+                if (buf1.size()>0) symmetry=buf1;
+            }
             else if (s == "restart") {
                 restart = true;
             }
@@ -532,6 +542,26 @@ struct CalculationParameters {
         }
 
         lo = molecule.smallest_length_scale();
+
+        // set molecular and computational point groups
+        // use highest point group unless specified by user
+
+        // complain if symmetry has been set to anything other than c1
+        if ((symmetry!="default" and symmetry!="c1") and localize) {
+        	error("\n\nsymmetry and localization cannot be used at the same time\n"
+        			"switch from local to canonical orbitals (keyword canon)\n\n");
+        }
+
+        // no symmetry keyword specified
+    	if (symmetry=="default") {
+    		if (localize) symmetry="c1";
+    		else symmetry=molecule.pointgroup_;
+
+    	// symmetry keyword specified without pointgroup
+    	} else if (symmetry=="full") {
+    		symmetry=molecule.pointgroup_;
+    	}
+
     }
 
     void print(World& world) const {
@@ -592,6 +622,7 @@ struct CalculationParameters {
             madness::print("  localized orbitals ", loctype);
         else
             madness::print("  canonical orbitals ");
+        madness::print("   comp. point group ", symmetry);
         if (derivatives)
             madness::print("    calc derivatives ");
         if (dipole)

src/apps/chem/Makefile.am

@@ -22,8 +22,9 @@ thisinclude_HEADERS = correlationfactor.h molecule.h molecularbasis.h \
                       corepotential.h atomutil.h SCF.h xcfunctional.h \
                       mp2.h nemo.h potentialmanager.h gth_pseudopotential.h \
                       molecular_optimizer.h projector.h \
-                      SCFOperators.h CCStructures.h \
-                      electronic_correlation_factor.h cheminfo.h vibanal.h molopt.h TDHF.h CC2.h CCPotentials.h AC.h GuessFactory.h
+                      SCFOperators.h CCStructures.h pointgroupoperator.h pointgroupsymmetry.h \
+                      electronic_correlation_factor.h cheminfo.h vibanal.h molopt.h TDHF.h \
+                      CC2.h CCPotentials.h AC.h GuessFactory.h 
 
 testxc_SOURCES = testxc.cc xcfunctional.h
 testxc_LDADD = libMADchem.la $(MRALIBS)
@@ -39,8 +40,9 @@ plotxc_LDADD = libMADchem.la $(MRALIBS)
 
 libMADchem_la_SOURCES = correlationfactor.cc molecule.cc molecularbasis.cc vibanal.cc \
                        corepotential.cc atomutil.cc lda.cc cheminfo.cc \
-                       distpm.cc SCF.cc gth_pseudopotential.cc nemo.cc mp2.cc pcm.cc\
-                       SCFOperators.cc xcfunctional_ldaonly.cc TDHF.cc CCStructures.cc CC2.cc CCPotentials.cc AC.cc GuessFactory.cc\
+                       distpm.cc SCF.cc gth_pseudopotential.cc nemo.cc mp2.cc pcm.cc pointgroupsymmetry.cc \
+                       SCFOperators.cc xcfunctional_ldaonly.cc TDHF.cc CCStructures.cc CC2.cc CCPotentials.cc \
+                       AC.cc GuessFactory.cc\
                        $(thisinclude_HEADERS)
 libMADchem_la_LDFLAGS = -version-info 0:0:0
 

src/apps/chem/SCFOperators.cc

@@ -187,6 +187,11 @@ real_function_3d Coulomb::compute_potential(const madness::Nemo* nemo) const {
         print("number of electrons in Coulomb",nel);
     }
     density.truncate();
+//    if (do_R2) {
+//        density.print_size("density with R2");
+//    } else {
+//        density.print_size("density without R2");
+//    }
     return nemo->get_calc()->make_coulomb_potential(density);
 }
 

src/apps/chem/SCFProtocol.h

@@ -102,7 +102,10 @@ class SCFProtocol {
         }
     }
 
-    bool finished() const {return converged;}
+    bool finished() const {
+        if(world.rank()==0) printf("\nending protocol at time %8.1fs \n",wall_time());
+    	return converged;
+    }
 
     /// go to the next level
     SCFProtocol& operator++() {
@@ -113,7 +116,6 @@ class SCFProtocol {
             if(world.rank()==0) print("protocol: thresh",thresh,"econv ",econv,"dconv",dconv);
         } else {
             converged=true;
-            if(world.rank()==0) printf("\nending protocol at time %8.1fs \n",wall_time());
         }
 
         // update restart data on file