Merge branch 'develop' of github-cDFT:jimlutsko/classicalDFT into dev…

…elop

CMakeLists.txt

@@ -1,6 +1,8 @@
 cmake_minimum_required(VERSION 3.10)
 
 set(USE_OMP true)
+set(USE_MPI true)
+set(USE_SLEPC true)
 
 set(USE_FMT_WEIGHTS_BEFORE_JUN_2021 false)
 set(USE_FMT_WEIGHTS_BEFORE_JAN_2023 false)
@@ -45,6 +47,8 @@ set(all_SRCS
   ${PROJECT_SOURCE_DIR}/src/DDFT.cpp
   ${PROJECT_SOURCE_DIR}/src/Arnoldi.cpp
   ${PROJECT_SOURCE_DIR}/src/Eigenvalues.cpp
+  ${PROJECT_SOURCE_DIR}/src/DFT_Petsc.cpp
+  ${PROJECT_SOURCE_DIR}/src/DFT_Slepc.cpp
   ${PROJECT_SOURCE_DIR}/src/Log_Det.cpp
   ${PROJECT_SOURCE_DIR}/src/Interaction.cpp
   ${PROJECT_SOURCE_DIR}/src/GitSHA1.cpp
@@ -66,6 +70,34 @@ link_directories("/usr/local/lib")
 
 target_link_libraries(${PROJECT_NAME} PUBLIC Boost::serialization)
 
+message("USE_SLEPC = " ${USE_SLEPC})
+
+if(USE_SLEPC)
+  find_package(PkgConfig REQUIRED) 
+  set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig:$ENV{PETSC_DIR}/lib/pkgconfig:$ENV{PKG_CONFIG_PATH}")
+  pkg_check_modules(PETSC_PKG REQUIRED IMPORTED_TARGET $ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig/PETSc.pc)
+  pkg_check_modules(SLEPC_PKG REQUIRED IMPORTED_TARGET $ENV{SLEPC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig/slepc.pc)
+  target_compile_definitions(${PROJECT_NAME} PUBLIC USE_SLEPC)
+  target_link_libraries(${PROJECT_NAME} PUBLIC PkgConfig::PETSC_PKG PkgConfig::SLEPC_PKG)
+endif()
+
+message("USE_MPI = " ${USE_MPI})
+
+if(USE_MPI)
+  find_package(MPI REQUIRED)
+  include_directories(${MPI_INCLUDE_PATH})
+
+  if(MPI_COMPILE_FLAGS)
+    set_target_properties(${PROJECT_NAME} PROPERTIES COMPILE_FLAGS "${MPI_COMPILE_FLAGS}")
+  endif()
+  if(MPI_LINK_FLAGS)
+    set_target_properties(${PROJECT_NAME} PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif()
+
+  #target_compile_definitions(${PROJECT_NAME} PUBLIC USE_MPI)
+  target_link_libraries(${PROJECT_NAME} PUBLIC MPI::MPI_CXX)
+endif()
+
 # enable omp if it exists
 
 message("USE_OMP = " ${USE_OMP})

additionalCode/DFT_Armadillo.cpp

@@ -16,6 +16,9 @@ using namespace std;
 #define v1_DATA (*(static_cast<arma::vec*>(v1.data_)))
 #define v2_DATA (*(static_cast<arma::vec*>(v2.data_)))
 
+
+const char *DFT_Vec::get_library_name() { return "Armadillo";}
+
 DFT_Vec::DFT_Vec(unsigned N) : data_(new arma::vec(N)){}
 DFT_Vec::DFT_Vec(const DFT_Vec& v) : data_(new arma::vec(v.size())) {DATA = v_DATA;}
 DFT_Vec::DFT_Vec(): data_(new arma::vec(1)) {}

additionalCode/DFT_Eigen.cpp

@@ -22,6 +22,7 @@ using namespace Eigen;
 #define v1_DATA (*(static_cast<Eigen::VectorXd*>(v1.data_)))
 #define v2_DATA (*(static_cast<Eigen::VectorXd*>(v2.data_)))
 
+const char *DFT_Vec::get_library_name() { return "Eigen";}
 
 DFT_Vec::DFT_Vec(unsigned N) : data_(new Eigen::VectorXd(N)){}
 DFT_Vec::DFT_Vec(const DFT_Vec& v) : data_(new Eigen::VectorXd(v_DATA)){} //(v.size())) { DATA = v_DATA;}
@@ -71,6 +72,7 @@ double DFT_Vec::accu() const { return DATA.sum();}
 void DFT_Vec::save(ofstream &of) const
 {
   of << DATA.size() << endl;
+  of << setprecision(20);
   for(long i=0;i<DATA.size();i++)
     of << DATA[i] << " ";
   of << endl;

include/DFT.h

@@ -56,6 +56,7 @@ class DFT : public Dynamical_Matrix
   // Accessors
   const Lattice& get_lattice() const {return allSpecies_.front()->getLattice();}
   const Density& getDensity(int species) const {return allSpecies_[species]->getDensity();}
+  const Density& get_density(int species) const {return allSpecies_[species]->getDensity();}
   Interaction_Base* getInteraction(int which) { return Interactions_[which];}
 
   Species *getSpecies(int s) const { return allSpecies_[s];}
@@ -79,6 +80,17 @@ class DFT : public Dynamical_Matrix
   void convert_dF_to_alias_derivs(vector<DFT_Vec> &x_);
   void convert_dF_to_alias_derivs();
 
+  void dF_to_H_dot_dF()
+  {
+    for(int s = 0; s<allSpecies_.size(); s++)
+      {
+	DFT_Vec &dF = getDF(s);
+	DFT_Vec ff(dF.size());
+	matrix_dot_v1(dF,ff);
+	dF.set(ff);
+      }
+  }
+
   void set_offset(bool val) { offset_ = val;}
 
   // A few actions  
@@ -87,8 +99,8 @@ class DFT : public Dynamical_Matrix
   void write_density_vtk(int i, string &of) const {allSpecies_[i]->getDensity().write_VTK_File
       (of);}
   void perturb_density(DFT_Vec& perturbation, int species = 0) { allSpecies_[species]->perturb_density(perturbation);} 
-  double calculateFreeEnergyAndDerivatives(bool onlyFex = false);
-  double evaluate(bool onlyFex = false) { return calculateFreeEnergyAndDerivatives(onlyFex);}  
+  double calculateFreeEnergyAndDerivatives(bool onlyFex = false, bool H_dot_Force = false);
+  double evaluate(bool onlyFex = false, bool H_dot_Force = false) { return calculateFreeEnergyAndDerivatives(onlyFex,H_dot_Force);}  
 
   // Bulk Thermodynamics
 
@@ -143,6 +155,8 @@ class DFT : public Dynamical_Matrix
 
   void set_full_hessian(bool full) { full_hessian_ = full;}
 
+  double get_rms_force() const { return rms_force_;}
+
   friend class boost::serialization::access;
   template<class Archive> void serialize(Archive & ar, const unsigned int version)
   {
@@ -162,14 +176,14 @@ class DFT : public Dynamical_Matrix
   FMT *fmt_;
   vector<External_Field*> external_fields_;
 
-
   double F_id_  = 0.0; ///< Ideal part of free energy
   double F_ext_ = 0.0; ///< External field contribution to free energy (including chemical potential)
   double F_hs_  = 0.0; ///< Hard-sphere contribution to free energy
   double F_mf_  = 0.0; ///< Mean-field contribution to free energy
 
   mutable bool full_hessian_ = true; // used in matrix_holder to distinguish full hessian from excess hessian.
-
+  mutable double rms_force_ = 0.0;   // updated everytime the forces are calculated.
+
   bool offset_ = false;
 };
 

include/DFT_Factory.h

@@ -105,12 +105,12 @@ class DFT_Factory
       options_.read(SourceInput_.c_str(), verbose_);
 
     if(include_log_ && theLog_ == NULL) theLog_ = new Log(log_file_name_.c_str(),-1,-1,NULL, -1, true, verbose_);
-    if(verbose_ && theLog_ != NULL) *theLog_ << myColor::GREEN << "=================================" << myColor::RESET << endl << "#" << endl;
+    if(verbose_ && theLog_ != NULL) *theLog_ << myColor::GREEN << "#=================================" << myColor::RESET << endl << "#" << endl;
 
     if(verbose_ && theLog_ != NULL) *theLog_ << myColor::RED << myColor::BOLD << "Input parameters:" << myColor::RESET << endl <<  "#" << endl;
 
     if(verbose_ && theLog_ != NULL) options_.write(*theLog_);
-    if(verbose_ && theLog_ != NULL) *theLog_ <<  myColor::GREEN << "=================================" <<  myColor::RESET << endl;
+    if(verbose_ && theLog_ != NULL) *theLog_ <<  myColor::GREEN << "#=================================" <<  myColor::RESET << endl;
 
     if (dy_<0) dy_ = dx_;
     if (dz_<0) dz_ = dx_;
@@ -189,8 +189,8 @@ class DFT_Factory
       }
     /////////////////////////////////////////////////////
     // Report
-    if(verbose_ && theLog_ != NULL) *theLog_ <<  myColor::GREEN << "=================================" <<  myColor::RESET << endl;
-    if(verbose_ && theLog_ != NULL) *theLog_ <<  myColor::GREEN << "Summary:" <<  myColor::RESET << endl << endl;  
+    if(verbose_ && theLog_ != NULL) *theLog_ <<  myColor::GREEN << "#=================================" <<  myColor::RESET << endl;
+    if(verbose_ && theLog_ != NULL) *theLog_ <<  myColor::GREEN << "#Summary:" <<  myColor::RESET << endl << endl;  
 #ifdef USE_OMP
     if(verbose_ && theLog_ != NULL) *theLog_ << "\tomp max threads  : " << omp_get_max_threads() << endl;
 #endif      
@@ -253,7 +253,7 @@ class DFT_Factory
 
     /////////////////////////////////////////////////////
     // Report
-    if(verbose_ && theLog_ != NULL) *theLog_ <<  myColor::GREEN << "=================================" <<  myColor::RESET << endl;
+    if(verbose_ && theLog_ != NULL) *theLog_ <<  myColor::GREEN << "#=================================" <<  myColor::RESET << endl;
     if(verbose_ && theLog_ != NULL) *theLog_ <<  myColor::GREEN << "Temperature set to " << kT_ << " giving:" <<  myColor::RESET << endl << endl;  
     if(verbose_ && theLog_ != NULL) *theLog_ << "\tHSD1                        = " << hsd1_ << endl;
     if(potential1_ && interaction1_)
@@ -274,7 +274,7 @@ class DFT_Factory
     check();
     /////////////////////////////////////////////////////
     // Thermodynamics
-    if(verbose_ && theLog_ != NULL) *theLog_ <<  myColor::GREEN << "=================================" <<  myColor::RESET << endl;
+    if(verbose_ && theLog_ != NULL) *theLog_ <<  myColor::GREEN << "#=================================" <<  myColor::RESET << endl;
     if(verbose_ && theLog_ != NULL) *theLog_ <<  myColor::GREEN << "Thermodynamics:" <<  myColor::RESET << endl << endl;
 
     xv_ = xl_ = xs1_ = xs2_ = -1;

include/DFT_LinAlg.h

@@ -22,6 +22,8 @@ class DFT_Vec
   DFT_Vec();
   ~DFT_Vec();
 
+  static const char *get_library_name();
+
   DFT_Vec& operator= (const DFT_Vec& c){set(c); return *this;}
 
   void   set(unsigned pos, double val);
@@ -73,7 +75,6 @@ class DFT_Vec
 
   void Schur(const DFT_Vec &v1, const DFT_Vec &v2);
 
-
   // There may be value in having library-specific implementations  of these functions
   // but for now, I do not see it. 
   friend ostream &operator<<(ostream &of, const DFT_Vec &v)

include/DFT_Slepc.h

@@ -0,0 +1,102 @@
+#ifdef USE_SLEPC
+
+#include <slepceps.h>
+#include <slepcsys.h>
+#include <petscerror.h>
+
+#include "Dynamical_Matrix.h"
+#include "Log.h"
+
+class DFT_Petsc
+{
+ public:
+  DFT_Petsc(Dynamical_Matrix &dm, Log &theLog, int argc, char**argv);
+  ~DFT_Petsc();
+
+  PetscErrorCode Solve_g_x_equals_b(DFT_Vec & dft_x, DFT_Vec &dft_b, stringstream &ret);
+  virtual PetscErrorCode write_version_info(ostream &os);
+
+  void set_abs_tol(double x) { abs_tol_ = x;}
+  void set_rel_tol(double x) { rel_tol_ = x;}
+  void set_max_its(int x)    { max_its_ = x;}
+
+  void update_message(string s);
+
+  PetscErrorCode MatMultA1(Vec &x,Vec &y);
+  PetscErrorCode MatGetDiagonalA1(Vec &diag);
+  PetscErrorCode MatMultTransposeA1(Vec &x,Vec &y);
+
+  PetscErrorCode MatMultG1(Vec &x,Vec &y);
+  PetscErrorCode MatGetDiagonalG1(Vec &diag);
+
+  void DFT_to_SLEPC(const DFT_Vec &vin,  Vec& v) const;
+  void SLEPC_to_DFT(const  Vec &v, DFT_Vec &vout) const;
+
+ protected:
+  Dynamical_Matrix &dm_;
+  Log &theLog_;
+  long Ntot_      = 0;
+  long Ndynamic_  = 0;
+
+  double abs_tol_ = 1e-5;
+  double rel_tol_ = 1e-7;
+  int    max_its_ = 100000;
+
+  Mat A_;
+};
+
+
+class DFT_Slepc : public DFT_Petsc
+{
+ public:
+  DFT_Slepc(Dynamical_Matrix &dm, Log &theLog, int argc, char**argv);
+  ~DFT_Slepc();
+
+  virtual PetscErrorCode write_version_info(ostream &os);
+
+  int run_eigenproblem(int argc, char**argv);
+
+  void set_input_vec(string  filename);
+  void set_input_vec(const DFT_Vec &v);
+
+  void set_eps_tol(double v) { eps_tol_ = v;}
+  void set_num_eig(int n) { num_eigenvalues_ = n;}
+  void set_two_sided(bool two_sided) { two_sided_ = two_sided;}
+  void set_smallest(bool s) { smallest_ = s;}
+
+  PetscErrorCode MatMultA1(Vec &x,Vec &y);
+  PetscErrorCode MatGetDiagonalA1(Vec &diag);
+  PetscErrorCode MatMultTransposeA1(Vec &x,Vec &y);
+
+  PetscErrorCode check();
+  PetscErrorCode check(DFT_Vec &eigen, double kr);
+  PetscErrorCode DisplaySolution(EPS &eps);
+
+  PetscErrorCode init_from_dft_vector();
+  PetscErrorCode init_with_random_vector();
+
+  PetscErrorCode write_output_vectors_dft();
+
+  PetscErrorCode get_eigenvalue(double &eigenvalue, double &eigenvalue_imag, int which = 0) const;
+  PetscErrorCode get_eigenvector(DFT_Vec &eigenvec, DFT_Vec &eigenvector_imag, int which = 0) const;
+  PetscErrorCode get_eigenvector_left(DFT_Vec &eigenvec, DFT_Vec &eigenvector_imag, int which = 0) const;
+  PetscErrorCode get_number_converged(int &nconv) const;
+
+  PetscErrorCode compute_error(double &err, int which = 0) const;
+
+  double  get_eigenvalue (int which = 0) const {double  kr,ki; get_eigenvalue(kr,ki,which);  return kr;}
+  DFT_Vec get_eigenvector(int which = 0) const {DFT_Vec vr(Ntot_),vi(Ntot_); get_eigenvector(vr,vi,which); return vr;}
+  double  get_res_error(int which = 0) const {double err; compute_error(err, which); return err;}
+
+ protected:
+  DFT_Vec input_vector_dft_;
+
+  double eps_tol_      = 1e-4; // relative tolerance
+  int num_eigenvalues_ = 1;
+  bool two_sided_      = false;
+  bool smallest_       = true;
+
+  EPS eps_;
+};
+
+#endif //USE_SLEPC

include/Dynamical_Matrix.h

@@ -35,15 +35,18 @@ class Dynamical_Matrix
 
   long get_Ntot() const { return get_dimension(0)*get_dimension(1)*get_dimension(2);}
 
-  void set_verbose(bool v)            { dm_verbose_ = v;}
-  void set_use_squared_matrix(bool v) { use_squared_matrix_  = v;}
+  void set_verbose(bool v)                 { dm_verbose_ = v;}
+  void set_use_squared_matrix(bool v)      { use_squared_matrix_  = v;}
+  void set_dm_in_alias_coordinates(bool v) { dm_in_alias_coordinates_ = v;}
 
   bool get_use_squared_matrix() const { return use_squared_matrix_;}
+  bool is_matrix_in_alias_coordinates() const { return dm_in_alias_coordinates_;}
 
 protected:
 
-  bool dm_verbose_         = false;
-  bool use_squared_matrix_ = false;
+  bool dm_verbose_              = false;
+  bool use_squared_matrix_      = false;
+  bool dm_in_alias_coordinates_ = false; // Not compatible with DDFT dynamics (not implemented)
 };
 
 #endif

include/Log.h

@@ -38,7 +38,7 @@
 class LogStreamBuf: public std::stringbuf
 {
  public:
- LogStreamBuf(ofstream& log) : log_(log) {}
+  LogStreamBuf(ofstream& log) : log_(log) {}
  ~LogStreamBuf(){ if(pbase() != pptr()) putOutput();}
 
   // When we sync the stream with the output. 
@@ -49,7 +49,7 @@ class LogStreamBuf: public std::stringbuf
 
   void putOutput()
   {
-    cout << std::boolalpha;
+    //    cout << std::boolalpha;
     // Called by destructor.
     // destructor can not call virtual methods.
     const string &s = str();
@@ -126,10 +126,11 @@ class Log: public std::ostream
       if(verbose) *this << "*****************************************************************" << endl;
       if(verbose) if(prog != NULL) *this << prog << " version " << Major << "." << Minor << endl;
       if(verbose) *this << std::ctime(&now_time) << " " << endl;
-      if(verbose) *this << "Using:\tLib " << PROJECT_NAME << endl
+      if(verbose) *this << "#Using:\tLib " << PROJECT_NAME << endl
 	    << "\tversion: " << PROJECT_VER << endl
 	    << "\tgit revision: " << g_GIT_SHA1 << endl;
       //      *this << "Library built " << _TIMEZ_  << endl;
+      if(verbose) *this << "\tLinear Algebra library: " << DFT_Vec::get_library_name() << endl;
       if(verbose) if(numtasks > 0) *this << " MPI: numtasks = " << numtasks << endl;
       if(verbose) *this << "*****************************************************************" << endl  << endl;      
     }
@@ -148,7 +149,7 @@ class Log: public std::ostream
    */
   void write(Options &options)
   {
-    *this << "Input parameters:" << endl <<  "#" << endl;
+    *this << "#Input parameters:" << endl <<  "#" << endl;
     options.write(*this);
     *this << "=================================" << endl;
     this->flush();

include/Minimizer.h

@@ -31,6 +31,7 @@ class Minimizer
   void setForceTerminationCriterion(double v) {forceLimit_ = v;}
   void set_terminimation_criterion(double v) {forceLimit_ = v;}  
   void setVerbose(bool verbose) { verbose_ = verbose;}
+  void set_use_squared_forces(bool b) { use_squared_forces_ = b;}
 
   const DFT_Vec&  get_fixed_direction() {return fixed_direction_;}
   void set_fixed_direction(const DFT_Vec& fixed, bool already_using_density_alias);
@@ -42,7 +43,7 @@ class Minimizer
   virtual void   cleanup() {} // called when minimization is finished to allow decendents to clean up user output.
 
 protected:
-  double get_energy_and_forces() {calls_++;  return dft_->calculateFreeEnergyAndDerivatives(false);}
+  double get_energy_and_forces() {calls_++;  return dft_->calculateFreeEnergyAndDerivatives(false, use_squared_forces_);}
   virtual double getDF_DX();
   virtual double step() = 0;
   virtual bool   should_stop() const;
@@ -52,6 +53,7 @@ class Minimizer
   DFT *dft_ = NULL; // calculates energy and forces
 
   bool verbose_ = false;
+  bool use_squared_forces_ = false;
 
   int    calls_        = 0;
   int    step_counter_ = 0;