diff --git a/cmake_modules/FindCLANG_FORMAT.cmake b/cmake_modules/FindCLANG_FORMAT.cmake index 509b9171..1eb3b1d3 100644 --- a/cmake_modules/FindCLANG_FORMAT.cmake +++ b/cmake_modules/FindCLANG_FORMAT.cmake @@ -5,7 +5,7 @@ # CLANG_FORMAT_VERSION - The version of clang-format found find_program(CLANG_FORMAT_EXECUTABLE - NAMES clang-format clang-format-6.0 + NAMES clang-format clang-format-15.0 DOC "clang-format executable") mark_as_advanced(CLANG_FORMAT_EXECUTABLE) @@ -14,14 +14,14 @@ if (CLANG_FORMAT_EXECUTABLE) execute_process(COMMAND ${CLANG_FORMAT_EXECUTABLE} -version OUTPUT_VARIABLE clang_format_version ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE) - if (clang_format_version MATCHES "^clang-format version .*") - # clang_format_version sample: "clang-format version 6.0.0-1ubuntu2 (tags/RELEASE_600/final)" - string(REGEX REPLACE "clang-format version ([.0-9]+).*" "\\1" - CLANG_FORMAT_VERSION "${clang_format_version}") + message(STATUS ${clang_format_version}) + + if (clang_format_version MATCHES "clang-format version ([0-9]+)") + set(CLANG_FORMAT_FOUND_VERSION_MAJOR ${CMAKE_MATCH_1}) + message(STATUS "Found clang-format version major: ${CLANG_FORMAT_FOUND_VERSION_MAJOR}") endif() - if (NOT CLANG_FORMAT_VERSION STREQUAL "6.0.0" AND - NOT CLANG_FORMAT_VERSION STREQUAL "6.0.1") - message(SEND_ERROR "Wrong clang-format version. Please use clang-format 6.0.") + if (NOT CLANG_FORMAT_FOUND_VERSION_MAJOR STREQUAL "15") + message(SEND_ERROR "Wrong clang-format version. Please use clang-format 15") endif() endif() diff --git a/src/AndersonMix.cc b/src/AndersonMix.cc index 2892c9f9..02be2622 100644 --- a/src/AndersonMix.cc +++ b/src/AndersonMix.cc @@ -26,7 +26,7 @@ static const double min_det_mat = 0.01; static const double max_theta = 0.5; static const double min_theta = -3.; -//#define DEBUG 0 +// #define DEBUG 0 template AndersonMix::AndersonMix(const int m, const double beta, T& x) : m_(m), x_(x) @@ -236,7 +236,7 @@ void AndersonMix::update(T& f, T& work, ostream& os, const bool verbose) } } - //#ifdef DEBUG + // #ifdef DEBUG if (os.good() && mm_ > 0 && verbose) { os << "Anderson extrapolation:"; @@ -244,7 +244,7 @@ void AndersonMix::update(T& f, T& work, ostream& os, const bool verbose) os << " theta[" << j << "]=" << theta_[j]; os << endl; } - //#endif + // #endif } // update x_ diff --git a/src/AndersonMix.h b/src/AndersonMix.h index 92a1569c..6f390767 100644 --- a/src/AndersonMix.h +++ b/src/AndersonMix.h @@ -32,7 +32,7 @@ class AndersonMix : public Mixing static Timer update_tm_; - virtual void postprocessUpdate(){}; + virtual void postprocessUpdate() {}; T& x_; // current trial solution diff --git a/src/ChebyshevApproximationFunction.h b/src/ChebyshevApproximationFunction.h index 0fa71d80..aac5c6c3 100644 --- a/src/ChebyshevApproximationFunction.h +++ b/src/ChebyshevApproximationFunction.h @@ -17,8 +17,8 @@ class ChebyshevApproximationFunction private: public: - ChebyshevApproximationFunction(){}; - virtual ~ChebyshevApproximationFunction(){}; // destructor + ChebyshevApproximationFunction() {}; + virtual ~ChebyshevApproximationFunction() {}; // destructor virtual std::vector eval(const std::vector& x) = 0; }; diff --git a/src/ChebyshevApproximationInterface.h b/src/ChebyshevApproximationInterface.h index b0c3606f..a7fc9db2 100644 --- a/src/ChebyshevApproximationInterface.h +++ b/src/ChebyshevApproximationInterface.h @@ -49,7 +49,7 @@ class ChebyshevApproximationInterface public: // constructor - ChebyshevApproximationInterface(){}; + ChebyshevApproximationInterface() {}; // build the Chebyshev coefficients for the interval [a, b] void computeChebyshevCoeffs(); @@ -66,7 +66,7 @@ class ChebyshevApproximationInterface max_order_ = order; order_ = order; } - virtual ~ChebyshevApproximationInterface(){}; + virtual ~ChebyshevApproximationInterface() {}; static void printTimers(std::ostream& os) { diff --git a/src/Constraint.h b/src/Constraint.h index 7a7c51a0..d8cf779e 100644 --- a/src/Constraint.h +++ b/src/Constraint.h @@ -23,7 +23,7 @@ class Constraint names_; // names of atoms involved in the constraint public: - virtual ~Constraint(){}; + virtual ~Constraint() {}; virtual bool enforce(void) = 0; virtual bool project_out_forces(void) = 0; diff --git a/src/ConstraintSet.cc b/src/ConstraintSet.cc index 70981518..f4869caa 100644 --- a/src/ConstraintSet.cc +++ b/src/ConstraintSet.cc @@ -147,7 +147,7 @@ bool ConstraintSet::addConstraint(Ions& ions, const vector& argv) // check if a constraint (name1,name2) or (name2,name1) is defined if (pc->type() == "distance" && ((pc->names(0) == name1 && pc->names(1) == name2) - || (pc->names(1) == name1 && pc->names(0) == name2))) + || (pc->names(1) == name1 && pc->names(0) == name2))) found = true; } diff --git a/src/ConstraintSet.h b/src/ConstraintSet.h index a76dddbc..b9953344 100644 --- a/src/ConstraintSet.h +++ b/src/ConstraintSet.h @@ -35,7 +35,7 @@ class ConstraintSet bool addConstraint(Ions&, const std::vector& argv); public: - ConstraintSet(){}; + ConstraintSet() {}; void clear(); diff --git a/src/Control.cc b/src/Control.cc index 01843fa9..f2b0c330 100644 --- a/src/Control.cc +++ b/src/Control.cc @@ -271,7 +271,8 @@ void Control::print(std::ostream& os) os << " Multigrid preconditioning for wave functions:" << std::endl; os << " # of Multigrid levels : " << mg_levels_ << std::endl; os << " # of pre-smoothing steps : " << mg_npresmoothing_ << std::endl; - os << " # of post-smoothing steps : " << mg_npostsmoothing_ << std::endl; + os << " # of post-smoothing steps : " << mg_npostsmoothing_ + << std::endl; } else { @@ -503,7 +504,8 @@ void Control::sync(void) memset(&float_buffer[0], 0, size_float_buffer * sizeof(float)); } - auto bcast_check = [](int mpirc) { + auto bcast_check = [](int mpirc) + { if (mpirc != MPI_SUCCESS) { (*MPIdata::sout) << "MPI Bcast of Control failed!!!" << std::endl; @@ -1497,9 +1499,9 @@ void Control::setOptions(const boost::program_options::variables_map& vm) mg_npresmoothing_ = vm["Preconditioner.npresmoothing"].as(); mg_npostsmoothing_ = vm["Preconditioner.npostsmoothing"].as(); precond_precision_ = vm["Preconditioner.precision"].as(); - assert(precond_precision_==32 ||precond_precision_==64); + assert(precond_precision_ == 32 || precond_precision_ == 64); - precond_factor = vm["Quench.step_length"].as(); + precond_factor = vm["Quench.step_length"].as(); if (precond_factor < 0.) { switch (lap_type) @@ -1994,10 +1996,14 @@ void Control::setROMOptions(const boost::program_options::variables_map& vm) else if (str.compare("none") == 0) rom_pri_option.rom_stage = ROMStage::UNSUPPORTED; - rom_pri_option.restart_file_fmt = vm["ROM.offline.restart_filefmt"].as(); - rom_pri_option.restart_file_minidx = vm["ROM.offline.restart_min_idx"].as(); - rom_pri_option.restart_file_maxidx = vm["ROM.offline.restart_max_idx"].as(); - rom_pri_option.basis_file = vm["ROM.offline.basis_file"].as(); + rom_pri_option.restart_file_fmt + = vm["ROM.offline.restart_filefmt"].as(); + rom_pri_option.restart_file_minidx + = vm["ROM.offline.restart_min_idx"].as(); + rom_pri_option.restart_file_maxidx + = vm["ROM.offline.restart_max_idx"].as(); + rom_pri_option.basis_file + = vm["ROM.offline.basis_file"].as(); str = vm["ROM.offline.variable"].as(); if (str.compare("orbitals") == 0) @@ -2007,14 +2013,19 @@ void Control::setROMOptions(const boost::program_options::variables_map& vm) else rom_pri_option.variable = ROMVariable::NONE; - rom_pri_option.save_librom_snapshot = vm["ROM.offline.save_librom_snapshot"].as(); - rom_pri_option.librom_snapshot_freq = vm["ROM.offline.librom_snapshot_freq"].as(); + rom_pri_option.save_librom_snapshot + = vm["ROM.offline.save_librom_snapshot"].as(); + rom_pri_option.librom_snapshot_freq + = vm["ROM.offline.librom_snapshot_freq"].as(); rom_pri_option.compare_md = vm["ROM.basis.compare_md"].as(); - rom_pri_option.num_orbbasis = vm["ROM.basis.number_of_orbital_basis"].as(); - rom_pri_option.num_potbasis = vm["ROM.basis.number_of_potential_basis"].as(); - rom_pri_option.pot_rom_file = vm["ROM.potential_rom_file"].as(); - } // onpe0 + rom_pri_option.num_orbbasis + = vm["ROM.basis.number_of_orbital_basis"].as(); + rom_pri_option.num_potbasis + = vm["ROM.basis.number_of_potential_basis"].as(); + rom_pri_option.pot_rom_file + = vm["ROM.potential_rom_file"].as(); + } // onpe0 // synchronize all processors syncROMOptions(); @@ -2031,7 +2042,8 @@ void Control::syncROMOptions() mmpi.bcast(rom_pri_option.basis_file, comm_global_); mmpi.bcast(rom_pri_option.pot_rom_file, comm_global_); - auto bcast_check = [](int mpirc) { + auto bcast_check = [](int mpirc) + { if (mpirc != MPI_SUCCESS) { (*MPIdata::sout) << "MPI Bcast of Control failed!!!" << std::endl; @@ -2044,31 +2056,38 @@ void Control::syncROMOptions() mpirc = MPI_Bcast(&rom_stage, 1, MPI_SHORT, 0, comm_global_); bcast_check(mpirc); - mpirc = MPI_Bcast(&rom_pri_option.restart_file_minidx, 1, MPI_INT, 0, comm_global_); + mpirc = MPI_Bcast( + &rom_pri_option.restart_file_minidx, 1, MPI_INT, 0, comm_global_); bcast_check(mpirc); - mpirc = MPI_Bcast(&rom_pri_option.restart_file_maxidx, 1, MPI_INT, 0, comm_global_); + mpirc = MPI_Bcast( + &rom_pri_option.restart_file_maxidx, 1, MPI_INT, 0, comm_global_); bcast_check(mpirc); - mpirc = MPI_Bcast(&rom_pri_option.save_librom_snapshot, 1, MPI_C_BOOL, 0, comm_global_); + mpirc = MPI_Bcast( + &rom_pri_option.save_librom_snapshot, 1, MPI_C_BOOL, 0, comm_global_); bcast_check(mpirc); - mpirc = MPI_Bcast(&rom_pri_option.librom_snapshot_freq, 1, MPI_INT, 0, comm_global_); + mpirc = MPI_Bcast( + &rom_pri_option.librom_snapshot_freq, 1, MPI_INT, 0, comm_global_); bcast_check(mpirc); short rom_var = (short)static_cast(rom_pri_option.variable); - mpirc = MPI_Bcast(&rom_var, 1, MPI_SHORT, 0, comm_global_); + mpirc = MPI_Bcast(&rom_var, 1, MPI_SHORT, 0, comm_global_); bcast_check(mpirc); rom_pri_option.rom_stage = static_cast(rom_stage); - rom_pri_option.variable = static_cast(rom_var); + rom_pri_option.variable = static_cast(rom_var); - mpirc = MPI_Bcast(&rom_pri_option.compare_md, 1, MPI_C_BOOL, 0, comm_global_); + mpirc + = MPI_Bcast(&rom_pri_option.compare_md, 1, MPI_C_BOOL, 0, comm_global_); bcast_check(mpirc); - mpirc = MPI_Bcast(&rom_pri_option.num_orbbasis, 1, MPI_INT, 0, comm_global_); + mpirc + = MPI_Bcast(&rom_pri_option.num_orbbasis, 1, MPI_INT, 0, comm_global_); bcast_check(mpirc); - mpirc = MPI_Bcast(&rom_pri_option.num_potbasis, 1, MPI_INT, 0, comm_global_); + mpirc + = MPI_Bcast(&rom_pri_option.num_potbasis, 1, MPI_INT, 0, comm_global_); bcast_check(mpirc); } diff --git a/src/Control.h b/src/Control.h index 6ae93692..1f97472d 100644 --- a/src/Control.h +++ b/src/Control.h @@ -216,7 +216,7 @@ class Control Control(); - ~Control(){}; + ~Control() {}; Control(const Control& ct) { (void)ct; }; void printRestartLink(); @@ -406,7 +406,6 @@ class Control short mg_npresmoothing_; short mg_npostsmoothing_; - // dielectric model for solvation short diel; // Parameters for MG solver/ preconditioner for Poisson problem @@ -483,7 +482,7 @@ class Control // Number of v-cycles for hartree solution short vh_its; - + // convergence tolerance for solving Poisson problem using PCG. float poisson_conv_tol; diff --git a/src/DMStrategy.h b/src/DMStrategy.h index 69322347..7839f92f 100644 --- a/src/DMStrategy.h +++ b/src/DMStrategy.h @@ -17,7 +17,7 @@ class DMStrategy virtual void initialize(OrbitalsType& orbitals) = 0; virtual int update(OrbitalsType& orbitals) = 0; - virtual ~DMStrategy(){}; + virtual ~DMStrategy() {}; // tells if strategy needs an updated H matrix // to update DM diff --git a/src/DavidsonSolver.cc b/src/DavidsonSolver.cc index 636de977..2b097be9 100644 --- a/src/DavidsonSolver.cc +++ b/src/DavidsonSolver.cc @@ -131,9 +131,9 @@ int DavidsonSolver::checkConvergence( // { // n++; // } -//#ifdef USE_MPI +// #ifdef USE_MPI // MPI_Bcast(&n, 1, MPI_INT, 0, comm_); -//#endif +// #endif // if (n <= numst) return; // // // build H matrix in basis of evect @@ -159,9 +159,9 @@ int DavidsonSolver::checkConvergence( // index = i; // } // } -//#ifdef USE_MPI +// #ifdef USE_MPI // MPI_Bcast(&index, 1, MPI_INT, 0, comm_); -//#endif +// #endif // if (index != n - 1) // { // if (onpe0 && ct.verbose > 2) diff --git a/src/DistMatrix/DistMatrix.h b/src/DistMatrix/DistMatrix.h index e07c8d56..317e7248 100644 --- a/src/DistMatrix/DistMatrix.h +++ b/src/DistMatrix/DistMatrix.h @@ -258,7 +258,7 @@ class DistMatrix return *this; } - ~DistMatrix() {} + ~DistMatrix() { } void identity(void); diff --git a/src/DistMatrix/DistVector.h b/src/DistMatrix/DistVector.h index 1f345191..d800f89c 100644 --- a/src/DistMatrix/DistVector.h +++ b/src/DistMatrix/DistVector.h @@ -31,7 +31,7 @@ class DistVector : public DistMatrix { } - DistVector(const int m) : DistMatrix("noname", m, 1) {} + DistVector(const int m) : DistMatrix("noname", m, 1) { } DistVector(const std::vector& v) : DistMatrix("noname", v.size(), 1) { diff --git a/src/DistMatrix/SubMatrices.cc b/src/DistMatrix/SubMatrices.cc index 93b7b08a..1d6c2782 100644 --- a/src/DistMatrix/SubMatrices.cc +++ b/src/DistMatrix/SubMatrices.cc @@ -12,7 +12,7 @@ #include -//#define DEBUG +// #define DEBUG namespace dist_matrix { @@ -161,12 +161,12 @@ void SubMatrices::gather(const DistMatrix& mat) // compute displacements std::vector remote_displ(npes_, 0); for (int pe = 1; pe < npes_; pe++) - remote_displ[pe] = (type_displ)( - submat_indexing_.getRemoteSize(pe - 1) + remote_displ[pe - 1]); + remote_displ[pe] = (type_displ)(submat_indexing_.getRemoteSize(pe - 1) + + remote_displ[pe - 1]); std::vector my_displ(npes_distmat_, 0); for (int pe = 1; pe < npes_distmat_; pe++) - my_displ[pe] = (type_displ)( - submat_indexing_.getMySize(pe - 1) + my_displ[pe - 1]); + my_displ[pe] = (type_displ)(submat_indexing_.getMySize(pe - 1) + + my_displ[pe - 1]); // build buffer for data array to send std::vector buf_remote_val(tot_remote_size, 0.); diff --git a/src/DistMatrix2SquareLocalMatrices.h b/src/DistMatrix2SquareLocalMatrices.h index 354c27d4..ed9e1513 100644 --- a/src/DistMatrix2SquareLocalMatrices.h +++ b/src/DistMatrix2SquareLocalMatrices.h @@ -61,7 +61,7 @@ class DistMatrix2SquareLocalMatrices "Work", gids, comm, empty_mat, *submat_indexing_)); } - ~DistMatrix2SquareLocalMatrices() {} + ~DistMatrix2SquareLocalMatrices() { } void convert(const dist_matrix::DistMatrix& dmat, SquareLocalMatrices& lmat); diff --git a/src/DistributedIonicData.h b/src/DistributedIonicData.h index 0132ccab..f99beaf1 100644 --- a/src/DistributedIonicData.h +++ b/src/DistributedIonicData.h @@ -39,7 +39,7 @@ class DistributedIonicData public: DistributedIonicData( const std::vector&, const std::vector&); - DistributedIonicData(){}; + DistributedIonicData() {}; int size() const { return (int)ion_names_.size(); }; diff --git a/src/DotProductDiagonal.h b/src/DotProductDiagonal.h index a5d5dec2..81e250dc 100644 --- a/src/DotProductDiagonal.h +++ b/src/DotProductDiagonal.h @@ -15,7 +15,7 @@ template class DotProductDiagonal : public DotProductManager { public: - DotProductDiagonal(){}; + DotProductDiagonal() {}; double dotProduct(T& phi0, const T& phi1) override; }; diff --git a/src/DotProductManager.h b/src/DotProductManager.h index 6491d696..ed3a68da 100644 --- a/src/DotProductManager.h +++ b/src/DotProductManager.h @@ -13,9 +13,9 @@ template class DotProductManager { public: - DotProductManager(){}; + DotProductManager() {}; - virtual ~DotProductManager(){}; + virtual ~DotProductManager() {}; virtual double dotProduct(T& a, const T& b) = 0; }; diff --git a/src/EigenDMStrategy.h b/src/EigenDMStrategy.h index 619456e3..318965b9 100644 --- a/src/EigenDMStrategy.h +++ b/src/EigenDMStrategy.h @@ -27,9 +27,9 @@ class EigenDMStrategy : public DMStrategy bool needH() const override { return true; } - void stripDM() override {} - void dressDM() override {} - void reset() override {} + void stripDM() override { } + void dressDM() override { } + void reset() override { } }; #endif diff --git a/src/ExtendedGridOrbitals.cc b/src/ExtendedGridOrbitals.cc index 48eafce6..ddd18207 100644 --- a/src/ExtendedGridOrbitals.cc +++ b/src/ExtendedGridOrbitals.cc @@ -723,7 +723,7 @@ int ExtendedGridOrbitals::write( hid_t dtype_id = outHdfDataType(ct.out_restart_info); dset_id = H5Dcreate2(file_id, datasetname.c_str(), dtype_id, - filespace, H5P_DEFAULT, plist_id, H5P_DEFAULT); + filespace, H5P_DEFAULT, plist_id, H5P_DEFAULT); if (dset_id < 0) { (*MPIdata::serr) << "ExtendedGridOrbitals::write_func_hdf5(), " @@ -1823,9 +1823,10 @@ void ExtendedGridOrbitals::set(std::string file_path, int rdim) CAROM::BasisReader reader(file_path); CAROM::Matrix* orbital_basis = reader.getSpatialBasis(rdim); - Control& ct = *(Control::instance()); - Mesh* mymesh = Mesh::instance(); - pb::GridFunc gf_psi(mymesh->grid(), ct.bcWF[0], ct.bcWF[1], ct.bcWF[2]); + Control& ct = *(Control::instance()); + Mesh* mymesh = Mesh::instance(); + pb::GridFunc gf_psi( + mymesh->grid(), ct.bcWF[0], ct.bcWF[1], ct.bcWF[2]); CAROM::Vector psi; for (int i = 0; i < rdim; ++i) { diff --git a/src/ExtendedGridOrbitals.h b/src/ExtendedGridOrbitals.h index e1fb667c..85f532ba 100644 --- a/src/ExtendedGridOrbitals.h +++ b/src/ExtendedGridOrbitals.h @@ -356,8 +356,8 @@ class ExtendedGridOrbitals : public Orbitals void app_mask(const int, pb::GridFunc&, const short) const {}; - void applyMask(const bool = false){}; - void applyCorrMask(const bool = false){}; + void applyMask(const bool = false) {}; + void applyCorrMask(const bool = false) {}; #ifdef HAVE_MAGMA void multiplyByMatrix( @@ -406,7 +406,6 @@ class ExtendedGridOrbitals : public Orbitals #ifdef MGMOL_HAS_LIBROM void set(std::string file_path, int rdim); #endif - }; #endif diff --git a/src/FIRE.cc b/src/FIRE.cc index 22628ab9..13b4c249 100644 --- a/src/FIRE.cc +++ b/src/FIRE.cc @@ -21,7 +21,7 @@ FIRE::FIRE(OrbitalsType** orbitals, Ions& ions, std::shared_ptr lrs, MasksSet& masks, Electrostatic& electrostat, const double dt, MGmol& strategy) : IonicAlgorithm( - orbitals, ions, rho, constraints, lrs, masks, strategy), + orbitals, ions, rho, constraints, lrs, masks, strategy), orbitals_(orbitals), ions_(ions), rho_(rho), diff --git a/src/FIRE.h b/src/FIRE.h index f25e44ce..66fee8bf 100644 --- a/src/FIRE.h +++ b/src/FIRE.h @@ -41,7 +41,7 @@ class FIRE : public IonicAlgorithm MasksSet& masks, Electrostatic& electrostat, const double dt, MGmol&); - ~FIRE() override{}; + ~FIRE() override {}; }; #endif diff --git a/src/FullyOccupiedNonOrthoDMStrategy.h b/src/FullyOccupiedNonOrthoDMStrategy.h index 01a66830..d882a047 100644 --- a/src/FullyOccupiedNonOrthoDMStrategy.h +++ b/src/FullyOccupiedNonOrthoDMStrategy.h @@ -26,9 +26,9 @@ class FullyOccupiedNonOrthoDMStrategy : public DMStrategy bool needH() const override { return false; } - void stripDM() override {} - void dressDM() override {} - void reset() override {} + void stripDM() override { } + void dressDM() override { } + void reset() override { } }; #endif diff --git a/src/GrassmanCG.h b/src/GrassmanCG.h index c0bd4c76..01c5b394 100644 --- a/src/GrassmanCG.h +++ b/src/GrassmanCG.h @@ -30,7 +30,7 @@ class GrassmanCG : public GrassmanLineMinimization ProjectedMatricesInterface* proj_matrices, MGmol* mgmol_strategy, Ions& ions, std::ostream& os) : GrassmanLineMinimization( - hamiltonian, proj_matrices, mgmol_strategy, ions, os) + hamiltonian, proj_matrices, mgmol_strategy, ions, os) { } diff --git a/src/GrassmanCGSparse.h b/src/GrassmanCGSparse.h index 510871ac..f3e4a16a 100644 --- a/src/GrassmanCGSparse.h +++ b/src/GrassmanCGSparse.h @@ -29,7 +29,7 @@ class GrassmanCGSparse : public GrassmanLineMinimization ProjectedMatricesInterface* proj_matrices, MGmol* mgmol_strategy, Ions& ions, std::ostream& os) : GrassmanLineMinimization( - hamiltonian, proj_matrices, mgmol_strategy, ions, os) + hamiltonian, proj_matrices, mgmol_strategy, ions, os) { } diff --git a/src/GrassmanLineMinimization.h b/src/GrassmanLineMinimization.h index 0d1c6842..21c96c30 100644 --- a/src/GrassmanLineMinimization.h +++ b/src/GrassmanLineMinimization.h @@ -72,7 +72,7 @@ class GrassmanLineMinimization : public OrbitalsStepper delete sdir_; } - void setup(T&) override{}; + void setup(T&) override {}; int updateWF(T& orbitals, Ions& ions, const double precond_factor, const bool orthof, T& work_orbitals, const bool accelerate, diff --git a/src/GridMask.cc b/src/GridMask.cc index f180f418..8a5aac1a 100644 --- a/src/GridMask.cc +++ b/src/GridMask.cc @@ -15,7 +15,7 @@ #include "GridMask.h" #include "SubCell.h" -//#define DEBUG 1 +// #define DEBUG 1 Timer GridMask::init_tm_("GridMask::init"); @@ -373,7 +373,7 @@ bool GridMask::overlap_on_pe( assert(gm.mask_not_zero_[level][iloc] >= -1); overlap_12 = (overlap_12 || ((mask_not_zero_[level][iloc] + 1) - && (gm.mask_not_zero_[level][iloc] + 1))); + && (gm.mask_not_zero_[level][iloc] + 1))); } return overlap_12; diff --git a/src/GridMaskMax.h b/src/GridMaskMax.h index c231eb2a..add1075b 100644 --- a/src/GridMaskMax.h +++ b/src/GridMaskMax.h @@ -25,7 +25,7 @@ class GridMaskMax : public GridMask public: GridMaskMax(const unsigned short nclevels, const unsigned short subdivx, const pb::Grid& mygrid) - : GridMask(nclevels, subdivx, mygrid){}; + : GridMask(nclevels, subdivx, mygrid) {}; void apply(float* u, const unsigned short level, const unsigned short iloc, const bool /*first_application*/ = false) override diff --git a/src/GridMaskMult.h b/src/GridMaskMult.h index 0b5d2603..8e8a1ed5 100644 --- a/src/GridMaskMult.h +++ b/src/GridMaskMult.h @@ -25,7 +25,7 @@ class GridMaskMult : public GridMask public: GridMaskMult(const unsigned short nclevels, const unsigned short subdivx, const pb::Grid& mygrid) - : GridMask(nclevels, subdivx, mygrid){}; + : GridMask(nclevels, subdivx, mygrid) {}; void apply(float* u, const unsigned short level, const unsigned short iloc, const bool /*first_application*/ = false) override diff --git a/src/Hamiltonian.cc b/src/Hamiltonian.cc index 601fdc24..28aff992 100644 --- a/src/Hamiltonian.cc +++ b/src/Hamiltonian.cc @@ -137,7 +137,7 @@ void Hamiltonian::applyLocal(const int ncolors, T& phi, T& hphi) { // This loop is not thread safe as GridFunc ghost values filling // MPI calls may conflicts (all use the same tag) - //#pragma omp parallel for + // #pragma omp parallel for for (int i = 0; i < ncolors; i++) { using memory_space_type = typename T::memory_space_type; diff --git a/src/Hartree_CG.h b/src/Hartree_CG.h index e1a59657..6d048970 100644 --- a/src/Hartree_CG.h +++ b/src/Hartree_CG.h @@ -33,7 +33,7 @@ class Hartree_CG : public Poisson }; // Destructor - ~Hartree_CG() override {} + ~Hartree_CG() override { } void setup(const short nu1, const short nu2, const short max_sweeps, const double tol, const short max_nlevels, diff --git a/src/Ion.h b/src/Ion.h index cf52db97..8b9ee046 100644 --- a/src/Ion.h +++ b/src/Ion.h @@ -85,7 +85,7 @@ class Ion Ion(const Species& species, IonData data); Ion(const Ion& ion); - ~Ion(){}; + ~Ion() {}; void init(const double crds[3], const double velocity[3], const bool lock); void setup(); diff --git a/src/IonicAlgorithm.h b/src/IonicAlgorithm.h index 6fd35c78..16a9ab98 100644 --- a/src/IonicAlgorithm.h +++ b/src/IonicAlgorithm.h @@ -51,7 +51,7 @@ class IonicAlgorithm ConstraintSet& constraints, std::shared_ptr lrs, MasksSet& masks, MGmol&); - virtual ~IonicAlgorithm(){}; + virtual ~IonicAlgorithm() {}; virtual void init(HDFrestart* h5f_file); virtual int run1step(); diff --git a/src/IonicStepper.h b/src/IonicStepper.h index c1880a60..2fc2ec65 100644 --- a/src/IonicStepper.h +++ b/src/IonicStepper.h @@ -49,6 +49,6 @@ class IonicStepper virtual int write_hdf5(HDFrestart&) = 0; virtual int init(HDFrestart&) = 0; - virtual ~IonicStepper(){}; + virtual ~IonicStepper() {}; }; #endif diff --git a/src/Ions.cc b/src/Ions.cc index 248323ae..147aabb1 100644 --- a/src/Ions.cc +++ b/src/Ions.cc @@ -204,9 +204,9 @@ void Ions::setup() updateListIons(); - //#ifndef NDEBUG - // checkUnicityLocalIons(); - //#endif + // #ifndef NDEBUG + // checkUnicityLocalIons(); + // #endif setupInteractingIons(); @@ -1023,9 +1023,9 @@ void Ions::initFromRestartFile(HDFrestart& h5_file) // update list ions updateListIons(); - //#ifndef NDEBUG - // checkUnicityLocalIons(); - //#endif + // #ifndef NDEBUG + // checkUnicityLocalIons(); + // #endif } void Ions::readRestartPositions(HDFrestart& h5_file) @@ -2942,19 +2942,19 @@ bool Ions::inListIons(const double x, const double y, const double z) // check to see if ion is in list if (((t[0] >= list_boundary_left_[0] && t[0] <= list_boundary_right_[0]) || ((t[0] - lattice_[0]) >= list_boundary_left_[0] - && (t[0] - lattice_[0]) <= list_boundary_right_[0]) + && (t[0] - lattice_[0]) <= list_boundary_right_[0]) || ((t[0] + lattice_[0]) >= list_boundary_left_[0] - && (t[0] + lattice_[0]) <= list_boundary_right_[0])) + && (t[0] + lattice_[0]) <= list_boundary_right_[0])) && ((t[1] >= list_boundary_left_[1] && t[1] <= list_boundary_right_[1]) - || ((t[1] - lattice_[1]) >= list_boundary_left_[1] - && (t[1] - lattice_[1]) <= list_boundary_right_[1]) - || ((t[1] + lattice_[1]) >= list_boundary_left_[1] - && (t[1] + lattice_[1]) <= list_boundary_right_[1])) + || ((t[1] - lattice_[1]) >= list_boundary_left_[1] + && (t[1] - lattice_[1]) <= list_boundary_right_[1]) + || ((t[1] + lattice_[1]) >= list_boundary_left_[1] + && (t[1] + lattice_[1]) <= list_boundary_right_[1])) && ((t[2] >= list_boundary_left_[2] && t[2] <= list_boundary_right_[2]) - || ((t[2] - lattice_[2]) >= list_boundary_left_[2] - && (t[2] - lattice_[2]) <= list_boundary_right_[2]) - || ((t[2] + lattice_[2]) >= list_boundary_left_[2] - && (t[2] + lattice_[2]) <= list_boundary_right_[2]))) + || ((t[2] - lattice_[2]) >= list_boundary_left_[2] + && (t[2] - lattice_[2]) <= list_boundary_right_[2]) + || ((t[2] + lattice_[2]) >= list_boundary_left_[2] + && (t[2] + lattice_[2]) <= list_boundary_right_[2]))) inList = true; return inList; diff --git a/src/KBPsiMatrixInterface.h b/src/KBPsiMatrixInterface.h index deba0ced..6094c074 100644 --- a/src/KBPsiMatrixInterface.h +++ b/src/KBPsiMatrixInterface.h @@ -28,9 +28,9 @@ class KBPsiMatrixInterface static Timer computeLocalElement_tm_; public: - KBPsiMatrixInterface() : iterative_index_(-1){}; + KBPsiMatrixInterface() : iterative_index_(-1) {}; - virtual ~KBPsiMatrixInterface(){}; + virtual ~KBPsiMatrixInterface() {}; int getIterativeIndex() const { return iterative_index_; } void setOutdated() { iterative_index_ = -1; } diff --git a/src/KBprojector.h b/src/KBprojector.h index 20fcc208..af3ceaf7 100644 --- a/src/KBprojector.h +++ b/src/KBprojector.h @@ -84,7 +84,7 @@ class KBprojector assert(species_.dim_nl() < 1000); } - ~KBprojector() {} + ~KBprojector() { } virtual void clear() = 0; @@ -101,14 +101,18 @@ class KBprojector // axpySket for templated destination type virtual void axpySKet( - const short iloc, const double alpha, double* const) const = 0; + const short iloc, const double alpha, double* const) const + = 0; virtual void axpySKet( - const short iloc, const double alpha, float* const) const = 0; + const short iloc, const double alpha, float* const) const + = 0; virtual void axpyKet(const short iloc, const std::vector& alpha, - double* const dst) const = 0; + double* const dst) const + = 0; virtual void axpyKet(const short iloc, const std::vector& alpha, - float* const dst) const = 0; + float* const dst) const + = 0; bool onlyOneProjector() const { diff --git a/src/KBprojectorSparse.cc b/src/KBprojectorSparse.cc index 38b4c9bb..cb0f38a7 100644 --- a/src/KBprojectorSparse.cc +++ b/src/KBprojectorSparse.cc @@ -504,8 +504,9 @@ void KBprojectorSparse::setDProjector(const short iloc, const int icount) r3[jcount] = (KBPROJDTYPE)(t1 * z * x); r4[jcount] = (KBPROJDTYPE)(0.5 * t1 * (x * x - y * y)); - r5[jcount] = (KBPROJDTYPE)( - 0.5 * t1 * (sqrt3 * z * z - rr * inv_sqrt3)); + r5[jcount] = (KBPROJDTYPE)(0.5 * t1 + * (sqrt3 * z * z + - rr * inv_sqrt3)); } else { @@ -612,20 +613,22 @@ void KBprojectorSparse::setFProjector(const short iloc, const int icount) const double y2 = y * y; const double z2 = z * z; - r1[jcount] = (KBPROJDTYPE)( - sqrt5 * t1 * (3. * x2 - y2) * y); - r2[jcount] = (KBPROJDTYPE)( - 2. * sqrt2 * sqrt3 * sqrt5 * t1 * x * y * z); - r3[jcount] = (KBPROJDTYPE)( - sqrt3 * t1 * y * (4. * z2 - x2 - y2)); - r4[jcount] = (KBPROJDTYPE)( - sqrt2 * t1 * z * (2. * z2 - 3. * x2 - 3 * y2)); - r5[jcount] = (KBPROJDTYPE)( - sqrt3 * t1 * x * (4. * z2 - x2 - y2)); - r6[jcount] = (KBPROJDTYPE)( - sqrt2 * sqrt3 * sqrt5 * t1 * (x2 - y2) * z); - r7[jcount] = (KBPROJDTYPE)( - sqrt5 * t1 * (x2 - 3. * y2) * x); + r1[jcount] = (KBPROJDTYPE)(sqrt5 * t1 + * (3. * x2 - y2) * y); + r2[jcount] + = (KBPROJDTYPE)(2. * sqrt2 * sqrt3 * sqrt5 * t1 + * x * y * z); + r3[jcount] = (KBPROJDTYPE)(sqrt3 * t1 * y + * (4. * z2 - x2 - y2)); + r4[jcount] + = (KBPROJDTYPE)(sqrt2 * t1 * z + * (2. * z2 - 3. * x2 - 3 * y2)); + r5[jcount] = (KBPROJDTYPE)(sqrt3 * t1 * x + * (4. * z2 - x2 - y2)); + r6[jcount] = (KBPROJDTYPE)(sqrt2 * sqrt3 * sqrt5 + * t1 * (x2 - y2) * z); + r7[jcount] = (KBPROJDTYPE)(sqrt5 * t1 + * (x2 - 3. * y2) * x); } else { diff --git a/src/LBFGS.cc b/src/LBFGS.cc index e8064202..03dc6287 100644 --- a/src/LBFGS.cc +++ b/src/LBFGS.cc @@ -22,7 +22,7 @@ LBFGS::LBFGS(OrbitalsType** orbitals, Ions& ions, MasksSet& masks, MasksSet& corrmasks, Electrostatic& electrostat, const double dt, MGmol& strategy) : IonicAlgorithm( - orbitals, ions, rho, constraints, lrs, masks, strategy), + orbitals, ions, rho, constraints, lrs, masks, strategy), orbitals_(orbitals), ions_(ions), rho_(rho), @@ -70,8 +70,8 @@ void LBFGS::setup(const double dt) MasksSet* ref_corrmasks = ref_corrmasks_ ? ref_corrmasks_.get() : nullptr; ref_orbitals_ = std::shared_ptr( new OrbitalsType("LBFGS_ref", mygrid, mymesh->subdivx(), ct.numst, - ct.bcWF, (*orbitals_)->getProjMatrices(), ref_lrs_, ref_masks, - ref_corrmasks, local_cluster_)); + ct.bcWF, (*orbitals_)->getProjMatrices(), ref_lrs_, ref_masks, + ref_corrmasks, local_cluster_)); ref_orbitals_->assign(**orbitals_); } diff --git a/src/LDAonGridLibXC.cc b/src/LDAonGridLibXC.cc index b32d97ad..f5921f67 100644 --- a/src/LDAonGridLibXC.cc +++ b/src/LDAonGridLibXC.cc @@ -47,7 +47,7 @@ double LDAonGridLibXC::getExc() const // int ione=1; double exc = mygrid.vel() * LinearAlgebraUtils::MPdot( - np, &rho_.rho_[0][0], &exc_[0]); + np, &rho_.rho_[0][0], &exc_[0]); double sum = 0.; MGmol_MPI& mmpi = *(MGmol_MPI::instance()); diff --git a/src/LocGridOrbitals.cc b/src/LocGridOrbitals.cc index 92793a46..8b22ada8 100644 --- a/src/LocGridOrbitals.cc +++ b/src/LocGridOrbitals.cc @@ -1110,7 +1110,7 @@ int LocGridOrbitals::write( hid_t dtype_id = outHdfDataType(ct.out_restart_info); dset_id = H5Dcreate2(file_id, datasetname.c_str(), dtype_id, - filespace, H5P_DEFAULT, plist_id, H5P_DEFAULT); + filespace, H5P_DEFAULT, plist_id, H5P_DEFAULT); if (dset_id < 0) { (*MPIdata::serr) << "LocGridOrbitals::write_func_hdf5(), " @@ -2002,7 +2002,7 @@ void LocGridOrbitals::orthonormalize2states( { tmp[1] += vel * block_vector_.dot( - color_ic, color_jc, iloc); + color_ic, color_jc, iloc); } } } @@ -2062,7 +2062,7 @@ void LocGridOrbitals::orthonormalize2states( { tmp[1] += vel * block_vector_.dot( - color_ic, color_jc, iloc); + color_ic, color_jc, iloc); } } } @@ -2237,7 +2237,7 @@ void LocGridOrbitals::normalize() { diagS[gid] += vel * static_cast( - block_vector_.dot(color, color, iloc)); + block_vector_.dot(color, color, iloc)); } } } diff --git a/src/LocalMatrices2ReplicatedMatrix.h b/src/LocalMatrices2ReplicatedMatrix.h index cf073ed1..555a5137 100644 --- a/src/LocalMatrices2ReplicatedMatrix.h +++ b/src/LocalMatrices2ReplicatedMatrix.h @@ -42,7 +42,7 @@ class LocalMatrices2ReplicatedMatrix return pinstance_; } - LocalMatrices2ReplicatedMatrix() {} + LocalMatrices2ReplicatedMatrix() { } static void setup(const std::vector>& gids) { diff --git a/src/LocalizationRegions.cc b/src/LocalizationRegions.cc index 4f83fe01..47aaf209 100644 --- a/src/LocalizationRegions.cc +++ b/src/LocalizationRegions.cc @@ -1052,7 +1052,7 @@ void LocalizationRegions::computeSubdivOverlapGids() // subdiv center double cc[3] = { 0.5 * (2. * subdom_lower_left_[0] - + (2 * iloc + 1) * loc_length), + + (2 * iloc + 1) * loc_length), 0.5 * (2. * subdom_lower_left_[1] + subdom_div_lattice_[1]), 0.5 * (2. * subdom_lower_left_[2] + subdom_div_lattice_[2]) }; double t[3] = { center[0], center[1], center[2] }; diff --git a/src/LocalizationRegions.h b/src/LocalizationRegions.h index c2aad575..f835be49 100644 --- a/src/LocalizationRegions.h +++ b/src/LocalizationRegions.h @@ -176,7 +176,7 @@ class LocalizationRegions int globalNumLRs() const { return nglobal_; } - virtual ~LocalizationRegions() {} + virtual ~LocalizationRegions() { } virtual void setup() { diff --git a/src/MD_IonicStepper.cc b/src/MD_IonicStepper.cc index a1b4acb1..bb42e26e 100644 --- a/src/MD_IonicStepper.cc +++ b/src/MD_IonicStepper.cc @@ -265,8 +265,8 @@ int MD_IonicStepper::run() taup_[3 * ia + j] = factor * (2. * tau0_[3 * ia + j] - taum_[3 * ia + j] - + dt_ * dt_ * invmass * fion_[3 * ia + j] - + 0.5 * dt_ * gamma_ * taum_[3 * ia + j]); + + dt_ * dt_ * invmass * fion_[3 * ia + j] + + 0.5 * dt_ * gamma_ * taum_[3 * ia + j]); } } else @@ -402,7 +402,7 @@ void MD_IonicStepper::updateWithVelocityScaling(const double lambda) taup_[3 * ia + j] += lambda * (tau0_[3 * ia + j] - taum_[3 * ia + j] - + dt_ * dt_ * invmass * fion_[3 * ia + j]); + + dt_ * dt_ * invmass * fion_[3 * ia + j]); } } } diff --git a/src/MGOrbitalsPreconditioning.cc b/src/MGOrbitalsPreconditioning.cc index cdd6c44c..7ddbb244 100644 --- a/src/MGOrbitalsPreconditioning.cc +++ b/src/MGOrbitalsPreconditioning.cc @@ -28,8 +28,9 @@ MGOrbitalsPreconditioning::MGOrbitalsPreconditioning( Mesh* mymesh = Mesh::instance(); const pb::Grid& mygrid(mymesh->grid()); - precond_ = std::make_shared>( - lap_type_, mg_levels_, ct.mg_npresmoothing_, ct.mg_npostsmoothing_, mygrid, ct.bcWF); + precond_ + = std::make_shared>(lap_type_, mg_levels_, + ct.mg_npresmoothing_, ct.mg_npostsmoothing_, mygrid, ct.bcWF); } template diff --git a/src/MGmol.cc b/src/MGmol.cc index 0dcd2092..cec599a5 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -1135,13 +1135,14 @@ void MGmol::dumpRestart() #ifdef MGMOL_HAS_LIBROM // Save orbital snapshots - if (ct.getROMOptions().save_librom_snapshot > 0 && ct.AtomsDynamic() == AtomsDynamicType::Quench) + if (ct.getROMOptions().save_librom_snapshot > 0 + && ct.AtomsDynamic() == AtomsDynamicType::Quench) { - ierr = save_orbital_snapshot( - filename, *current_orbitals_); + ierr = save_orbital_snapshot(filename, *current_orbitals_); if (ierr < 0) - os_ << "WARNING: writing ROM snapshot data failed!!!" << std::endl; + os_ << "WARNING: writing ROM snapshot data failed!!!" + << std::endl; } #endif } @@ -1457,7 +1458,8 @@ void MGmol::getAtomicNumbers(std::vector& an) } template -void MGmol::updateDMandEnergy(OrbitalsType& orbitals, Ions& ions, double& eks) +void MGmol::updateDMandEnergy( + OrbitalsType& orbitals, Ions& ions, double& eks) { // initialize electronic density rho_->update(orbitals); @@ -1473,14 +1475,15 @@ void MGmol::updateDMandEnergy(OrbitalsType& orbitals, Ions& ions, std::shared_ptr> dm_strategy( DMStrategyFactory>::create(comm_, os_, ions, - rho_.get(), energy_.get(), electrostat_.get(), - hamiltonian_.get(), this, proj_matrices_.get(), &orbitals)); + rho_.get(), energy_.get(), electrostat_.get(), hamiltonian_.get(), + this, proj_matrices_.get(), &orbitals)); dm_strategy->update(orbitals); // evaluate energy and forces double ts = 0.; - eks = energy_->evaluateTotal(ts, proj_matrices_.get(), ions, orbitals, 2, os_); + eks = energy_->evaluateTotal( + ts, proj_matrices_.get(), ions, orbitals, 2, os_); } template diff --git a/src/MGmol.h b/src/MGmol.h index 6738c58f..db29f4bc 100644 --- a/src/MGmol.h +++ b/src/MGmol.h @@ -360,7 +360,8 @@ class MGmol : public MGmolInterface #ifdef MGMOL_HAS_LIBROM int save_orbital_snapshot(std::string file_path, OrbitalsType& orbitals); - void project_orbital(std::string file_path, int rdim, OrbitalsType& orbitals); + void project_orbital( + std::string file_path, int rdim, OrbitalsType& orbitals); #endif void updateDMandEnergy(OrbitalsType& orbitals, Ions& ions, double& eks); }; diff --git a/src/MGmolInterface.h b/src/MGmolInterface.h index 3c2384fd..a0138846 100644 --- a/src/MGmolInterface.h +++ b/src/MGmolInterface.h @@ -16,9 +16,9 @@ class MGmolInterface { public: - MGmolInterface() {} + MGmolInterface() { } - virtual ~MGmolInterface() {} + virtual ~MGmolInterface() { } virtual int setupFromInput(const std::string input_file) = 0; virtual int setupLRs(const std::string input_file) = 0; diff --git a/src/MGmol_prototypes.h b/src/MGmol_prototypes.h index 02e37890..e0f39301 100644 --- a/src/MGmol_prototypes.h +++ b/src/MGmol_prototypes.h @@ -10,8 +10,8 @@ #ifndef MGMOL_PROTOTYPES_H #define MGMOL_PROTOTYPES_H -#include "mgmol_config.h" #include "global.h" +#include "mgmol_config.h" #include namespace po = boost::program_options; @@ -27,7 +27,7 @@ int read_config(int argc, char** argv, std::string& lrs_filename, std::string& constraints_filename, float& total_spin, bool& with_spin); #ifdef MGMOL_HAS_LIBROM -void setupROMConfigOption(po::options_description &rom_cfg); +void setupROMConfigOption(po::options_description& rom_cfg); #endif #endif diff --git a/src/MLWFTransform.cc b/src/MLWFTransform.cc index 6bfccd29..191defa1 100644 --- a/src/MLWFTransform.cc +++ b/src/MLWFTransform.cc @@ -7,7 +7,7 @@ // This file is part of MGmol. For details, see https://github.com/llnl/mgmol. // Please also read this link https://github.com/llnl/mgmol/LICENSE -//#define DEBUG 1 +// #define DEBUG 1 #include "MGmol_blas1.h" #include @@ -138,4 +138,4 @@ void MLWFTransform::printTransform() fb.close(); } } -//#endif +// #endif diff --git a/src/MLWFTransform.h b/src/MLWFTransform.h index c35e9594..ceee27a1 100644 --- a/src/MLWFTransform.h +++ b/src/MLWFTransform.h @@ -29,7 +29,7 @@ class MLWFTransform : public OrbitalsTransform MLWFTransform(const int nst, const Vector3D& origin, const Vector3D& ll); - ~MLWFTransform() override {} + ~MLWFTransform() override { } void printTransform(); }; diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc index 97bcc825..04e4954b 100644 --- a/src/MVPSolver.cc +++ b/src/MVPSolver.cc @@ -235,7 +235,7 @@ int MVPSolver::solve(OrbitalsType& orbitals) ProjectedMatrices* current_proj_mat = (inner_it == 0) ? dynamic_cast*>( - orbitals.getProjMatrices()) + orbitals.getProjMatrices()) : proj_mat_work_; const int printE = (ct.verbose > 1) ? 1 : 0; diff --git a/src/MVP_DMStrategy.h b/src/MVP_DMStrategy.h index 04dbcad8..325e35bb 100644 --- a/src/MVP_DMStrategy.h +++ b/src/MVP_DMStrategy.h @@ -50,7 +50,7 @@ class MVP_DMStrategy : public DMStrategy const std::vector>& overlappingGids, ProjectedMatricesInterface* proj_matrices, const bool use_old_dm); - void initialize(OrbitalsType&) override{}; + void initialize(OrbitalsType&) override {}; int update(OrbitalsType& orbitals) override; // H is updated with MVP loop, so no need to compute it outside @@ -60,7 +60,7 @@ class MVP_DMStrategy : public DMStrategy void dressDM() override; - void reset() override {} + void reset() override { } }; #endif diff --git a/src/Masks4Orbitals.cc b/src/Masks4Orbitals.cc index fe055649..7eb6ce37 100644 --- a/src/Masks4Orbitals.cc +++ b/src/Masks4Orbitals.cc @@ -19,7 +19,7 @@ Masks4Orbitals::Masks4Orbitals( associateGids2Masks(overlap_gids); } -Masks4Orbitals::~Masks4Orbitals() {} +Masks4Orbitals::~Masks4Orbitals() { } void Masks4Orbitals::associateGids2Masks(const vector& overlap_gids) { diff --git a/src/Mesh.h b/src/Mesh.h index 123a40fe..803ab5a1 100644 --- a/src/Mesh.h +++ b/src/Mesh.h @@ -52,7 +52,7 @@ class Mesh delete myGrid_; delete myPEenv_; }; - Mesh(const Mesh&){}; + Mesh(const Mesh&) {}; public: static Mesh* instance() diff --git a/src/Mixing.h b/src/Mixing.h index 442c793c..5a3e8776 100644 --- a/src/Mixing.h +++ b/src/Mixing.h @@ -16,9 +16,9 @@ template class Mixing { public: - Mixing(){}; + Mixing() {}; - virtual ~Mixing(){}; + virtual ~Mixing() {}; virtual void update(T& res, T& work, std::ostream& os, const bool verbose) = 0; diff --git a/src/MultiDistanceConstraint.cc b/src/MultiDistanceConstraint.cc index acc622e5..19bfaa18 100644 --- a/src/MultiDistanceConstraint.cc +++ b/src/MultiDistanceConstraint.cc @@ -277,7 +277,7 @@ bool MultiDistanceConstraint::project_out_forces(void) // projection proj += dd * (dx * u_fion_[p][0] + dy * u_fion_[p][1] - + dz * u_fion_[p][2]); + + dz * u_fion_[p][2]); // gradient a_[p][0] += dd * dx; a_[p][1] += dd * dy; @@ -354,7 +354,7 @@ double MultiDistanceConstraint::projected_force(void) const // projection de -= dd * (dx * u_fion_[p][0] + dy * u_fion_[p][1] - + dz * u_fion_[p][2]); + + dz * u_fion_[p][2]); // gradient a_[p][0] += dd * dx; a_[p][1] += dd * dy; diff --git a/src/NOLMOTransform.cc b/src/NOLMOTransform.cc index 3f271bee..b53d1a10 100644 --- a/src/NOLMOTransform.cc +++ b/src/NOLMOTransform.cc @@ -17,7 +17,7 @@ #include "NOLMOTransform.h" #include "mputils.h" -//#define DEBUG 1 +// #define DEBUG 1 //////////////////////////////////////////////////////////////////////////////// diff --git a/src/Orbitals.h b/src/Orbitals.h index bc524aff..c744a7b2 100644 --- a/src/Orbitals.h +++ b/src/Orbitals.h @@ -26,7 +26,7 @@ class Orbitals Orbitals() { iterative_index_ = -10; } - virtual ~Orbitals(){}; + virtual ~Orbitals() {}; Orbitals(const Orbitals& A, const bool copy_data) { diff --git a/src/OrbitalsExtrapolation.h b/src/OrbitalsExtrapolation.h index e2d0ff0c..bd134bfe 100644 --- a/src/OrbitalsExtrapolation.h +++ b/src/OrbitalsExtrapolation.h @@ -22,7 +22,7 @@ template class OrbitalsExtrapolation { public: - OrbitalsExtrapolation() : orbitals_minus1_(nullptr) {} + OrbitalsExtrapolation() : orbitals_minus1_(nullptr) { } virtual ~OrbitalsExtrapolation(); diff --git a/src/OrbitalsExtrapolationOrder2.h b/src/OrbitalsExtrapolationOrder2.h index 82c71c67..b26228d1 100644 --- a/src/OrbitalsExtrapolationOrder2.h +++ b/src/OrbitalsExtrapolationOrder2.h @@ -16,7 +16,7 @@ template class OrbitalsExtrapolationOrder2 : public OrbitalsExtrapolation { public: - OrbitalsExtrapolationOrder2() {} + OrbitalsExtrapolationOrder2() { } void extrapolate_orbitals( OrbitalsType** orbitals, OrbitalsType* new_orbitals) override; diff --git a/src/OrbitalsExtrapolationOrder3.h b/src/OrbitalsExtrapolationOrder3.h index d6b7b4b4..d78554aa 100644 --- a/src/OrbitalsExtrapolationOrder3.h +++ b/src/OrbitalsExtrapolationOrder3.h @@ -24,7 +24,7 @@ class OrbitalsExtrapolationOrder3 : public OrbitalsExtrapolation OrbitalsExtrapolationOrder3() : initial_orbitals_minus2_(nullptr), orbitals_minus1_(nullptr), - orbitals_minus2_(nullptr){}; + orbitals_minus2_(nullptr) {}; ~OrbitalsExtrapolationOrder3() override { diff --git a/src/OrbitalsPreconditioning.h b/src/OrbitalsPreconditioning.h index 0501388c..80d22ec9 100644 --- a/src/OrbitalsPreconditioning.h +++ b/src/OrbitalsPreconditioning.h @@ -17,9 +17,9 @@ template class OrbitalsPreconditioning { public: - OrbitalsPreconditioning(){}; + OrbitalsPreconditioning() {}; - virtual ~OrbitalsPreconditioning(){}; + virtual ~OrbitalsPreconditioning() {}; virtual void setup(OrbitalsType& orbitals, MasksSet*, const std::shared_ptr&) diff --git a/src/OrbitalsStepper.h b/src/OrbitalsStepper.h index ffeedaf6..715912a6 100644 --- a/src/OrbitalsStepper.h +++ b/src/OrbitalsStepper.h @@ -16,9 +16,9 @@ template class OrbitalsStepper { public: - OrbitalsStepper() {} + OrbitalsStepper() { } - virtual ~OrbitalsStepper() {} + virtual ~OrbitalsStepper() { } virtual void setup(T&) = 0; @@ -27,7 +27,7 @@ class OrbitalsStepper const bool print_res, const double atol) = 0; - virtual void restartMixing(){}; + virtual void restartMixing() {}; }; #endif diff --git a/src/OrthoAndersonMix.h b/src/OrthoAndersonMix.h index 4df1e12b..0c14cabe 100644 --- a/src/OrthoAndersonMix.h +++ b/src/OrthoAndersonMix.h @@ -28,7 +28,7 @@ class OrthoAndersonMix : public AndersonMix { } - ~OrthoAndersonMix() override{}; + ~OrthoAndersonMix() override {}; }; #endif diff --git a/src/PBEonGridLibXC.cc b/src/PBEonGridLibXC.cc index 57231070..30f81b9e 100644 --- a/src/PBEonGridLibXC.cc +++ b/src/PBEonGridLibXC.cc @@ -150,7 +150,7 @@ double PBEonGridLibXC::getExc() const double exc = mygrid.vel() * LinearAlgebraUtils::MPdot( - np_, &rho_.rho_[0][0], &exc_[0]); + np_, &rho_.rho_[0][0], &exc_[0]); double sum = 0.; MGmol_MPI& mmpi = *(MGmol_MPI::instance()); diff --git a/src/Poisson.h b/src/Poisson.h index f7af6468..7ceabc80 100644 --- a/src/Poisson.h +++ b/src/Poisson.h @@ -66,7 +66,7 @@ class Poisson : public PoissonInterface double IntVhRhoc(void) const { return Int_vhrhoc_; } double IntVhRho_old(void) const { return Int_vhrho_old_; } - virtual void set_rhod(pb::GridFunc* /*rhod*/){}; + virtual void set_rhod(pb::GridFunc* /*rhod*/) {}; void set_vh(const pb::GridFunc& vh) { (*vh_) = vh; }; void set_vh(const std::vector& vh) { @@ -89,9 +89,11 @@ class Poisson : public PoissonInterface Int_vhrhoc_ = vel * vh_->gdot(rhoc); } - virtual void applyOperator(pb::GridFunc &vh, pb::GridFunc &lhs) + virtual void applyOperator( + pb::GridFunc& vh, pb::GridFunc& lhs) { - std::cerr << "ERROR: Abstract method Poisson::applyOperator()" << std::endl; + std::cerr << "ERROR: Abstract method Poisson::applyOperator()" + << std::endl; MPI_Abort(MPI_COMM_WORLD, 0); } }; diff --git a/src/PoissonInterface.h b/src/PoissonInterface.h index 37787395..659f34bb 100644 --- a/src/PoissonInterface.h +++ b/src/PoissonInterface.h @@ -21,7 +21,7 @@ class PoissonInterface static Timer poisson_tm_; public: - virtual ~PoissonInterface() {} + virtual ~PoissonInterface() { } static void printTimers(std::ostream& os) { poisson_tm_.print(os); } }; diff --git a/src/PolakRibiereSolver.cc b/src/PolakRibiereSolver.cc index ca25b1cc..b76d3bf0 100644 --- a/src/PolakRibiereSolver.cc +++ b/src/PolakRibiereSolver.cc @@ -160,7 +160,7 @@ bool PolakRibiereSolver::checkWolfeConditions( // const double dk =-1.*r_k_->dotProduct(*p_k_); const double dkm1 = -1. * r_km1_->dotProduct(*p_k_, - 2); // inverse(S) already included in r_km1_ + 2); // inverse(S) already included in r_km1_ const bool wolfe0 = (trial_step_energy <= eks_history_[0] + sigma_a_ * alpha_k * dkm1); diff --git a/src/Potentials.cc b/src/Potentials.cc index 960fd74d..768878c6 100644 --- a/src/Potentials.cc +++ b/src/Potentials.cc @@ -144,7 +144,7 @@ double Potentials::updateVtot(const std::vector>& rho) vtot_[idx] = (POTDTYPE)(ha2ry * ((double)v_nuc_[idx] + (double)v_ext_[idx] - + (double)vh_rho_[idx] + (double)vxc_rho_[idx])); + + (double)vh_rho_[idx] + (double)vxc_rho_[idx])); } double two = ha2ry; if (diel_) @@ -211,7 +211,7 @@ double Potentials::computeDeltaV(const std::vector>& rho) dv_[idx] = (POTDTYPE)(ha2ry * ((double)v_nuc_[idx] + (double)v_ext_[idx] - + (double)vh_rho_[idx] + (double)vxc_rho_[idx])); + + (double)vh_rho_[idx] + (double)vxc_rho_[idx])); } double two = ha2ry; if (diel_) @@ -895,11 +895,11 @@ void Potentials::initBackground() if (fabs(background_charge_) < 1.e-10) background_charge_ = 0.; } -void Potentials::evalIonDensityOnSamplePts( - Ions& ions, const std::vector &local_idx, std::vector &sampled_rhoc) +void Potentials::evalIonDensityOnSamplePts(Ions& ions, + const std::vector& local_idx, std::vector& sampled_rhoc) { Mesh* mymesh = Mesh::instance(); - MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); const pb::Grid& mygrid = mymesh->grid(); const char flag_filter = pot_type(0); @@ -907,8 +907,9 @@ void Potentials::evalIonDensityOnSamplePts( { if (onpe0) { - std::cout << "Potentials::evalIonDensityOnSamplePts - flag_filter s is not supported" - << std::endl; + std::cout << "Potentials::evalIonDensityOnSamplePts - flag_filter " + "s is not supported" + << std::endl; } mmpi.abort(); } @@ -931,8 +932,9 @@ void Potentials::evalIonDensityOnSamplePts( return; } -void Potentials::initializeRadialDataOnSampledPts( - const Vector3D& position, const Species& sp, const std::vector &local_idx, std::vector &sampled_rhoc) +void Potentials::initializeRadialDataOnSampledPts(const Vector3D& position, + const Species& sp, const std::vector& local_idx, + std::vector& sampled_rhoc) { assert(local_idx.size() == sampled_rhoc.size()); @@ -960,7 +962,7 @@ void Potentials::initializeRadialDataOnSampledPts( const RadialInter& lpot(sp.local_pot()); const Vector3D lattice(mygrid.ll(0), mygrid.ll(1), mygrid.ll(2)); - for(int k = 0; k < local_idx.size(); k++) + for (int k = 0; k < local_idx.size(); k++) { /* local_idx provides offset. @@ -979,8 +981,7 @@ void Potentials::initializeRadialDataOnSampledPts( /* accumulate ion species density */ const double r = position.minimage(point, lattice, ct.bcPoisson); - if (r < lrad) - sampled_rhoc[k] += sp.getRhoComp(r); + if (r < lrad) sampled_rhoc[k] += sp.getRhoComp(r); } return; diff --git a/src/Potentials.h b/src/Potentials.h index 3eb41fd3..9f3e7561 100644 --- a/src/Potentials.h +++ b/src/Potentials.h @@ -100,8 +100,9 @@ class Potentials void initializeSupersampledRadialDataOnMesh( const Vector3D& position, const Species& sp); - void initializeRadialDataOnSampledPts( - const Vector3D& position, const Species& sp, const std::vector &local_idx, std::vector &sampled_rhoc); + void initializeRadialDataOnSampledPts(const Vector3D& position, + const Species& sp, const std::vector& local_idx, + std::vector& sampled_rhoc); void rescaleRhoComp(); @@ -208,7 +209,8 @@ class Potentials void resetVhRho2Backup() { vh_rho_ = vh_rho_backup_; } - void evalIonDensityOnSamplePts(Ions& ions, const std::vector &local_idx, std::vector &sampled_rhoc); + void evalIonDensityOnSamplePts(Ions& ions, + const std::vector& local_idx, std::vector& sampled_rhoc); #ifdef HAVE_TRICUBIC void readExternalPot(const string filename, const char type); diff --git a/src/Preconditioning.cc b/src/Preconditioning.cc index 491b9974..8e320dae 100644 --- a/src/Preconditioning.cc +++ b/src/Preconditioning.cc @@ -12,16 +12,14 @@ template Preconditioning::Preconditioning(const short lap_type, const short maxlevels, - const short npresmooth, const short npostsmooth, - const pb::Grid& grid, const short bcWF[3]): - max_levels_(maxlevels), - npresmooth_(npresmooth), - npostsmooth_(npostsmooth) + const short npresmooth, const short npostsmooth, const pb::Grid& grid, + const short bcWF[3]) + : max_levels_(maxlevels), npresmooth_(npresmooth), npostsmooth_(npostsmooth) { - assert(npresmooth_>=0); - assert(npostsmooth_>=0); - assert(npresmooth_<100); - assert(npostsmooth_<100); + assert(npresmooth_ >= 0); + assert(npostsmooth_ >= 0); + assert(npresmooth_ < 100); + assert(npostsmooth_ < 100); for (short i = 0; i < 3; i++) bc_[i] = bcWF[i]; @@ -160,7 +158,7 @@ void Preconditioning::mg(pb::GridFuncVector& gfv_v, assert(gfv_work_[level] != nullptr); short ncycl = npresmooth_; - if (level == max_levels_) ncycl = npresmooth_+npostsmooth_; + if (level == max_levels_) ncycl = npresmooth_ + npostsmooth_; const double jacobi_factor = jacobi_factor_[level]; diff --git a/src/Preconditioning.h b/src/Preconditioning.h index 33aad1fb..6de10f86 100644 --- a/src/Preconditioning.h +++ b/src/Preconditioning.h @@ -30,7 +30,6 @@ class Preconditioning #endif private: - // V-cycle parameters const short max_levels_; const short npresmooth_; @@ -52,8 +51,8 @@ class Preconditioning public: Preconditioning(const short lap_type, const short maxlevels, - const short npresmooth, const short npostsmooth, - const pb::Grid& grid, const short bc[3]); + const short npresmooth, const short npostsmooth, const pb::Grid& grid, + const short bc[3]); Preconditioning(const Preconditioning&) = delete; diff --git a/src/ProjectedMatrices.h b/src/ProjectedMatrices.h index fa5add85..704bb4d0 100644 --- a/src/ProjectedMatrices.h +++ b/src/ProjectedMatrices.h @@ -151,7 +151,7 @@ class ProjectedMatrices : public ProjectedMatricesInterface { localHl_->copy(slH); } - void clearSparseH() override {} + void clearSparseH() override { } void consolidateH() override; diff --git a/src/ProjectedMatricesInterface.h b/src/ProjectedMatricesInterface.h index 53508a2e..c888785f 100644 --- a/src/ProjectedMatricesInterface.h +++ b/src/ProjectedMatricesInterface.h @@ -138,7 +138,7 @@ class ProjectedMatricesInterface : public ChebyshevApproximationFunction return (this->*(funcptr_))(nodes); } - virtual ~ProjectedMatricesInterface(){}; + virtual ~ProjectedMatricesInterface() {}; virtual void setup(const std::vector>& global_indexes) = 0; @@ -176,12 +176,12 @@ class ProjectedMatricesInterface : public ChebyshevApproximationFunction virtual void consolidateH() = 0; // return density matrix (inverse Gram if no unoccupied states) - virtual SquareLocalMatrices& - getLocalX() const = 0; + virtual SquareLocalMatrices& getLocalX() const + = 0; // return S**(-1)*H (or B**(-1)*H with Mehrstellen) - virtual SquareLocalMatrices& - getLocalT() const = 0; + virtual SquareLocalMatrices& getLocalT() const + = 0; // returns local Gram Matrix // virtual SquareLocalmatrices& getLocalS()const=0; @@ -223,8 +223,8 @@ class ProjectedMatricesInterface : public ChebyshevApproximationFunction virtual void applyInvS( SquareLocalMatrices& mat) = 0; - virtual double getLinDependent2states( - int& st1, int& st2, const bool) const = 0; + virtual double getLinDependent2states(int& st1, int& st2, const bool) const + = 0; virtual double checkCond(const double tol, const bool flag = true) = 0; virtual double getNel() const = 0; virtual void updateThetaAndHB() = 0; diff --git a/src/ProjectedMatricesSparse.h b/src/ProjectedMatricesSparse.h index fdbb4f88..245bca94 100644 --- a/src/ProjectedMatricesSparse.h +++ b/src/ProjectedMatricesSparse.h @@ -358,9 +358,9 @@ class ProjectedMatricesSparse : public ProjectedMatricesInterface const VariableSizeMatrix& getH() { return *matHB_; } void assignH(const VariableSizeMatrix& matH) { *matHB_ = matH; } - void printEigenvalues(std::ostream& /*os*/) const {} - void printOccupations(std::ostream& /*os*/) const override {} - void setHB2H() override {} + void printEigenvalues(std::ostream& /*os*/) const { } + void printOccupations(std::ostream& /*os*/) const override { } + void setHB2H() override { } DensityMatrixSparse& getDM() { return *dm_; } }; diff --git a/src/RadialSpline.cc b/src/RadialSpline.cc index 92610f75..9bb57d00 100644 --- a/src/RadialSpline.cc +++ b/src/RadialSpline.cc @@ -40,7 +40,7 @@ void RadialSpline::spline(const double yp1, const double ypn) y2_[i] = (sig - 1.0) * invp; u[i] = (6.0 * ((y_[i + 1] - y_[i]) / (x_[i + 1] - x_[i]) - - (y_[i] - y_[i - 1]) / (x_[i] - x_[i - 1])) + - (y_[i] - y_[i - 1]) / (x_[i] - x_[i - 1])) / (x_[i + 1] - x_[i - 1]) - sig * u[i - 1]) * invp; diff --git a/src/ReplicatedMatrix.cc b/src/ReplicatedMatrix.cc index 2f266ab2..885372d1 100644 --- a/src/ReplicatedMatrix.cc +++ b/src/ReplicatedMatrix.cc @@ -121,7 +121,7 @@ ReplicatedMatrix& ReplicatedMatrix::operator=(const ReplicatedMatrix& rhs) return *this; } -ReplicatedMatrix::~ReplicatedMatrix() {} +ReplicatedMatrix::~ReplicatedMatrix() { } void ReplicatedMatrix::getsub( const ReplicatedMatrix& src, int m, int n, int ia, int ja) diff --git a/src/ReplicatedMatrix2SquareLocalMatrices.h b/src/ReplicatedMatrix2SquareLocalMatrices.h index b3e05005..8218a282 100644 --- a/src/ReplicatedMatrix2SquareLocalMatrices.h +++ b/src/ReplicatedMatrix2SquareLocalMatrices.h @@ -35,14 +35,14 @@ class ReplicatedMatrix2SquareLocalMatrices return pinstance_; } - ReplicatedMatrix2SquareLocalMatrices() {} + ReplicatedMatrix2SquareLocalMatrices() { } static void setup(const std::vector>& gids) { global_indexes_ = gids; } - ~ReplicatedMatrix2SquareLocalMatrices() {} + ~ReplicatedMatrix2SquareLocalMatrices() { } void convert(const ReplicatedMatrix& dmat, SquareLocalMatrices& lmat); diff --git a/src/Rho.h b/src/Rho.h index 11063914..0c7642c1 100644 --- a/src/Rho.h +++ b/src/Rho.h @@ -69,7 +69,7 @@ class Rho std::vector> rho_; Rho(); - ~Rho(){}; + ~Rho() {}; const OrthoType getOrthoType() { return orbitals_type_; } diff --git a/src/SinCosOps.cc b/src/SinCosOps.cc index cd3f7a8c..01387106 100644 --- a/src/SinCosOps.cc +++ b/src/SinCosOps.cc @@ -696,7 +696,7 @@ void SinCosOps::compute( const double alpha = (double)orbitals1.psi(color)[index] * (double)orbitals2.psi( - jstate)[index]; + jstate)[index]; atmp[0] += alpha * cosx[ix]; atmp[1] += alpha * sinx[ix]; atmp[2] += alpha * cosy[iy]; diff --git a/src/Species.h b/src/Species.h index 87289da1..be051879 100644 --- a/src/Species.h +++ b/src/Species.h @@ -171,7 +171,7 @@ class Species #endif } - ~Species() {} + ~Species() { } unsigned short getAtomicNumber() const { return atomic_number_; } double getMass() const { return mass_ * SCMASS; } // SCMASS = 1822.89 diff --git a/src/SpreadPenaltyInterface.h b/src/SpreadPenaltyInterface.h index 19d2ba3a..538a438a 100644 --- a/src/SpreadPenaltyInterface.h +++ b/src/SpreadPenaltyInterface.h @@ -14,9 +14,9 @@ template class SpreadPenaltyInterface { public: - SpreadPenaltyInterface() {} + SpreadPenaltyInterface() { } - virtual ~SpreadPenaltyInterface() {} + virtual ~SpreadPenaltyInterface() { } // add penalty functional contribution to residual virtual void addResidual(T& phi, T& res) = 0; diff --git a/src/SpreadsAndCenters.h b/src/SpreadsAndCenters.h index f340ddfd..b83aee1c 100644 --- a/src/SpreadsAndCenters.h +++ b/src/SpreadsAndCenters.h @@ -53,7 +53,7 @@ class SpreadsAndCenters std::vector>& matr); public: - virtual ~SpreadsAndCenters() {} + virtual ~SpreadsAndCenters() { } void computeCenters(std::vector& centers) const { diff --git a/src/SquareSubMatrix2DistMatrix.h b/src/SquareSubMatrix2DistMatrix.h index fd1c3113..46a74568 100644 --- a/src/SquareSubMatrix2DistMatrix.h +++ b/src/SquareSubMatrix2DistMatrix.h @@ -39,7 +39,7 @@ class SquareSubMatrix2DistMatrix return pinstance_; } - SquareSubMatrix2DistMatrix() {} + SquareSubMatrix2DistMatrix() { } template void convert(const SquareSubMatrix& lmat, diff --git a/src/SubCell.cc b/src/SubCell.cc index 545fb923..3d63ddd5 100644 --- a/src/SubCell.cc +++ b/src/SubCell.cc @@ -46,8 +46,8 @@ SubCell::SubCell( outer_radius_ = sqrt(0.25 * (subcell_dimensions_[0] * subcell_dimensions_[0] - + subcell_dimensions_[1] * subcell_dimensions_[1] - + subcell_dimensions_[2] * subcell_dimensions_[2])); + + subcell_dimensions_[1] * subcell_dimensions_[1] + + subcell_dimensions_[2] * subcell_dimensions_[2])); inner_radius_ = std::min(subcell_dimensions_[0], subcell_dimensions_[1]); inner_radius_ = std::min(inner_radius_, subcell_dimensions_[2]); inner_radius_ *= 0.5; diff --git a/src/SubspaceProjector.h b/src/SubspaceProjector.h index ee2199cf..ebf293f3 100644 --- a/src/SubspaceProjector.h +++ b/src/SubspaceProjector.h @@ -29,7 +29,7 @@ class SubspaceProjector public: SubspaceProjector(T& subspace); - ~SubspaceProjector() {} + ~SubspaceProjector() { } void projectOut( T&, SquareLocalMatrices* mask = nullptr); diff --git a/src/XCFunctional.h b/src/XCFunctional.h index 0092ff64..fa57f16c 100644 --- a/src/XCFunctional.h +++ b/src/XCFunctional.h @@ -103,7 +103,7 @@ class XCFunctional } // virtual destructor needed to ensure proper deallocation - virtual ~XCFunctional() {} + virtual ~XCFunctional() { } virtual void computeXC(void) = 0; }; diff --git a/src/XConGrid.h b/src/XConGrid.h index 50693940..16d31281 100644 --- a/src/XConGrid.h +++ b/src/XConGrid.h @@ -18,9 +18,9 @@ class XConGrid public: static Timer get_xc_tm_; - XConGrid(){}; + XConGrid() {}; - virtual ~XConGrid(){}; + virtual ~XConGrid() {}; virtual void update() = 0; diff --git a/src/global.h b/src/global.h index a9176b60..7e98c91f 100644 --- a/src/global.h +++ b/src/global.h @@ -13,7 +13,7 @@ // this file should be included in every MGmol file // to enable global definitions, macros, ... -//#include "mgmol_memory.h" +// #include "mgmol_memory.h" #ifdef MGMOL_USE_MIXEDP typedef float ORBDTYPE; diff --git a/src/linear_algebra/mputils.cc b/src/linear_algebra/mputils.cc index b49ba366..794de6ca 100644 --- a/src/linear_algebra/mputils.cc +++ b/src/linear_algebra/mputils.cc @@ -632,7 +632,8 @@ void LAU_D::MPgemm(const char transa, const char transb, const int m, dgemm_tm.start(); // Transform char to magma_trans_t - auto convert_to_magma_trans = [](const char trans) { + auto convert_to_magma_trans = [](const char trans) + { if ((trans == 'N') || trans == 'n') return MagmaNoTrans; else if ((trans == 'T') || trans == 't') diff --git a/src/local_matrices/LocalMatrices.h b/src/local_matrices/LocalMatrices.h index 615cfab1..2c2fc314 100644 --- a/src/local_matrices/LocalMatrices.h +++ b/src/local_matrices/LocalMatrices.h @@ -63,7 +63,7 @@ class LocalMatrices void copy(const bml_matrix_t* A); #endif - virtual ~LocalMatrices(){}; + virtual ~LocalMatrices() {}; short nmat() const { return nmat_; } @@ -171,7 +171,7 @@ class LocalMatrices print(os, iloc); } - static void printTimers(std::ostream& /*os*/) {} + static void printTimers(std::ostream& /*os*/) { } void applyMask(const LocalMatrices& mask); diff --git a/src/local_matrices/LocalVector.h b/src/local_matrices/LocalVector.h index 2cf8992d..16ca3463 100644 --- a/src/local_matrices/LocalVector.h +++ b/src/local_matrices/LocalVector.h @@ -12,9 +12,9 @@ class LocalVector std::vector data_; public: - LocalVector(const int n) : data_(n) {} + LocalVector(const int n) : data_(n) { } - LocalVector(const std::vector& v) : data_(v) {} + LocalVector(const std::vector& v) : data_(v) { } DataType* data() { return data_.data(); } const DataType* data() const { return data_.data(); } diff --git a/src/md.cc b/src/md.cc index cc5d1d01..70e76c43 100644 --- a/src/md.cc +++ b/src/md.cc @@ -405,14 +405,16 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) h5f_file_.reset(); } - bool ROM_MVP = (ct.getROMOptions().rom_stage == ROMStage::ONLINE_PINNED_H2O_3DOF); + bool ROM_MVP + = (ct.getROMOptions().rom_stage == ROMStage::ONLINE_PINNED_H2O_3DOF); #ifdef MGMOL_HAS_LIBROM // ROM - initialize orbitals and density matrix if (ROM_MVP) { if (onpe0) os_ << "Setup ROM MVP solver..." << std::endl; - ExtendedGridOrbitals** extended_orbitals = reinterpret_cast**>(orbitals); - (*extended_orbitals)->set(ct.getROMOptions().basis_file, ct.numst); + ExtendedGridOrbitals** extended_orbitals + = reinterpret_cast**>(orbitals); + (*extended_orbitals)->set(ct.getROMOptions().basis_file, ct.numst); (*extended_orbitals)->orthonormalizeLoewdin(); (*extended_orbitals)->setDataWithGhosts(true); (*extended_orbitals)->setIterativeIndex(10); @@ -436,7 +438,7 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) bool extrapolated_flag = true; if (ct.dt <= 0.) extrapolated_flag = false; - int librom_snapshot_freq = ct.getROMOptions().librom_snapshot_freq; + int librom_snapshot_freq = ct.getROMOptions().librom_snapshot_freq; if (librom_snapshot_freq == -1) librom_snapshot_freq = ct.md_print_freq; MDfiles md_files; @@ -474,10 +476,13 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) PinnedH2O H2O_molecule; if (ct.getROMOptions().rom_stage == ROMStage::ONLINE_PINNED_H2O_3DOF) { - if (onpe0) os_ << "Rotate Pinned H2O molecule in timestep " << mdstep << std::endl; + if (onpe0) + os_ << "Rotate Pinned H2O molecule in timestep " << mdstep + << std::endl; H2O_molecule.rotate(positions, anumbers); if (onpe0) H2O_molecule.print(os_); - ROM_ions = new Ions(positions, anumbers, lattice, ions_->getSpecies()); + ROM_ions + = new Ions(positions, anumbers, lattice, ions_->getSpecies()); setupPotentials(*ROM_ions); force_on_ions = false; } @@ -555,15 +560,15 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) << std::endl; // Compute forces - if (force_on_ions) - force(**orbitals, ions); + if (force_on_ions) force(**orbitals, ions); #ifdef MGMOL_HAS_LIBROM if (ct.getROMOptions().rom_stage == ROMStage::ONLINE_PINNED_H2O_3DOF) { force(**orbitals, *ROM_ions); // Pinned H2O 3 DOF - if (onpe0) os_ << "Transpose rotate the PinnedH2O molecule" << std::endl; + if (onpe0) + os_ << "Transpose rotate the PinnedH2O molecule" << std::endl; std::vector forces; ROM_ions->getForces(forces); H2O_molecule.transpose_rotate(positions, anumbers, forces); @@ -575,14 +580,19 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) { if (onpe0) { - os_ << "Projecting orbitals onto ROM subspaces to compare " - << ((ct.getROMOptions().compare_md) ? "MD dynamics" : "force") << std::endl; - os_ << "Loading ROM basis " << ct.getROMOptions().basis_file << std::endl; - os_ << "ROM basis dimension = " << ct.getROMOptions().num_orbbasis << std::endl; + os_ << "Projecting orbitals onto ROM subspaces to compare " + << ((ct.getROMOptions().compare_md) ? "MD dynamics" + : "force") + << std::endl; + os_ << "Loading ROM basis " << ct.getROMOptions().basis_file + << std::endl; + os_ << "ROM basis dimension = " + << ct.getROMOptions().num_orbbasis << std::endl; } // Project orbitals to ROM subspace - project_orbital(ct.getROMOptions().basis_file, ct.getROMOptions().num_orbbasis, **orbitals); + project_orbital(ct.getROMOptions().basis_file, + ct.getROMOptions().num_orbbasis, **orbitals); if (ct.getROMOptions().compare_md) { // overwrite ions force for MD time stepping @@ -591,13 +601,13 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) else { // write ROM_ions force for one-step comparison - ROM_ions = new Ions(positions, anumbers, lattice, ions_->getSpecies()); + ROM_ions = new Ions( + positions, anumbers, lattice, ions_->getSpecies()); force(**orbitals, *ROM_ions); std::string zero = "0"; if (ions_->getNumIons() < 256 || ct.verbose > 2) { if (ct.verbose > 0) ROM_ions->printForcesGlobal(os_); - } else if (zero.compare(ct.md_print_filename) == 0) { @@ -707,7 +717,7 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) // note: extrapolation is going to modify it! if ((ct.out_restart_info > 2) && (((md_iteration_ % ct.checkpoint) == 0) - || (mdstep == ct.num_MD_steps))) + || (mdstep == ct.num_MD_steps))) proj_matrices_->saveDM(); preWFextrapolation(); @@ -756,14 +766,16 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) #ifdef MGMOL_HAS_LIBROM // Save orbital snapshots - if (ct.getROMOptions().save_librom_snapshot > 0 && - md_iteration_ % librom_snapshot_freq == 0) + if (ct.getROMOptions().save_librom_snapshot > 0 + && md_iteration_ % librom_snapshot_freq == 0) { int ierr = save_orbital_snapshot( - ct.md_print_filename + "_mdstep" + std::to_string(mdstep), **orbitals); + ct.md_print_filename + "_mdstep" + std::to_string(mdstep), + **orbitals); if (ierr < 0) - os_ << "WARNING md(): writing ROM snapshot data failed!!!" << std::endl; + os_ << "WARNING md(): writing ROM snapshot data failed!!!" + << std::endl; } #endif diff --git a/src/numerical_kernels/rho.cc b/src/numerical_kernels/rho.cc index 77dba143..eb493ebb 100644 --- a/src/numerical_kernels/rho.cc +++ b/src/numerical_kernels/rho.cc @@ -11,7 +11,7 @@ #include "Timer.h" #include "numerical_kernels.h" -//#define WTIMERS +// #define WTIMERS #ifdef WTIMERS Timer nonOrthoRhoKernel_tm("nonOrthoRhoKernel"); diff --git a/src/pb/Delh4.h b/src/pb/Delh4.h index 4a02a7b5..f6fecc93 100644 --- a/src/pb/Delh4.h +++ b/src/pb/Delh4.h @@ -19,7 +19,7 @@ template class Delxh4 : public FDoper { public: - Delxh4(const Grid& mygrid) : FDoper(mygrid) {} + Delxh4(const Grid& mygrid) : FDoper(mygrid) { } // A->B void apply(GridFunc& A, GridFunc& B) override @@ -27,7 +27,7 @@ class Delxh4 : public FDoper this->del1_4th(A, B, 0); } - ~Delxh4() override{}; + ~Delxh4() override {}; static short minNumberGhosts() { return 2; } }; @@ -35,7 +35,7 @@ template class Delyh4 : public FDoper { public: - Delyh4(const Grid& mygrid) : FDoper(mygrid) {} + Delyh4(const Grid& mygrid) : FDoper(mygrid) { } // A->B void apply(GridFunc& A, GridFunc& B) override @@ -43,7 +43,7 @@ class Delyh4 : public FDoper this->del1_4th(A, B, 1); } - ~Delyh4() override{}; + ~Delyh4() override {}; static short minNumberGhosts() { return 2; } }; @@ -51,7 +51,7 @@ template class Delzh4 : public FDoper { public: - Delzh4(const Grid& mygrid) : FDoper(mygrid) {} + Delzh4(const Grid& mygrid) : FDoper(mygrid) { } // A->B void apply(GridFunc& A, GridFunc& B) override @@ -59,7 +59,7 @@ class Delzh4 : public FDoper this->del1_4th(A, B, 2); } - ~Delzh4() override{}; + ~Delzh4() override {}; static short minNumberGhosts() { return 2; } }; diff --git a/src/pb/DielFunc.h b/src/pb/DielFunc.h index f91b0c22..05d991b0 100644 --- a/src/pb/DielFunc.h +++ b/src/pb/DielFunc.h @@ -90,7 +90,7 @@ class DielFunc : public GridFunc drho0_ = A.drho0_; } - ~DielFunc() override{}; + ~DielFunc() override {}; DielFunc& operator=(const T); diff --git a/src/pb/FDkernels.cc b/src/pb/FDkernels.cc index 76c75753..6d6cca98 100644 --- a/src/pb/FDkernels.cc +++ b/src/pb/FDkernels.cc @@ -309,28 +309,37 @@ void FDkernelDel2_6th(const Grid& grid, ScalarType* v, ScalarType* u, for (int iz = 0; iz < dim2; iz++) { - u[iiz] = (ScalarType)( - c0 * (double)v[iiz] - - + c1x * ((double)v[iiz - incx] + (double)v[iiz + incx]) - + c1y * ((double)v[iiz - incy] + (double)v[iiz + incy]) - + c1z * ((double)v[iiz - 1] + (double)v[iiz + 1]) - - + c2x - * ((double)v[iiz - incx2] - + (double)v[iiz + incx2]) - + c2y - * ((double)v[iiz - incy2] - + (double)v[iiz + incy2]) - + c2z * ((double)v[iiz - 2] + (double)v[iiz + 2]) - - + c3x - * ((double)v[iiz - incx3] - + (double)v[iiz + incx3]) - + c3y - * ((double)v[iiz - incy3] - + (double)v[iiz + incy3]) - + c3z * ((double)v[iiz - 3] + (double)v[iiz + 3])); + u[iiz] = (ScalarType)(c0 * (double)v[iiz] + + + c1x + * ((double)v[iiz - incx] + + (double)v[iiz + incx]) + + c1y + * ((double)v[iiz - incy] + + (double)v[iiz + incy]) + + c1z + * ((double)v[iiz - 1] + + (double)v[iiz + 1]) + + + c2x + * ((double)v[iiz - incx2] + + (double)v[iiz + incx2]) + + c2y + * ((double)v[iiz - incy2] + + (double)v[iiz + incy2]) + + c2z + * ((double)v[iiz - 2] + + (double)v[iiz + 2]) + + + c3x + * ((double)v[iiz - incx3] + + (double)v[iiz + incx3]) + + c3y + * ((double)v[iiz - incy3] + + (double)v[iiz + incy3]) + + c3z + * ((double)v[iiz - 3] + + (double)v[iiz + 3])); iiz++; } @@ -371,7 +380,7 @@ void FDkernelDel2_8th(const Grid& grid, ScalarType* v, ScalarType* u, const double c0 = -2. * (c1x + c2x + c3x + c4x + c1y + c2y + c3y + c4y + c1z - + c2z + c3z + c4z); + + c2z + c3z + c4z); const int incx = grid.inc(0); const int incy = grid.inc(1); @@ -402,36 +411,47 @@ void FDkernelDel2_8th(const Grid& grid, ScalarType* v, ScalarType* u, for (int iz = 0; iz < dim2; iz++) { - u[iiz] = (ScalarType)( - c0 * (double)v[iiz] - - + c1x * ((double)v[iiz - incx] + (double)v[iiz + incx]) - + c1y * ((double)v[iiz - incy] + (double)v[iiz + incy]) - + c1z * ((double)v[iiz - 1] + (double)v[iiz + 1]) - - + c2x - * ((double)v[iiz - incx2] - + (double)v[iiz + incx2]) - + c2y - * ((double)v[iiz - incy2] - + (double)v[iiz + incy2]) - + c2z * ((double)v[iiz - 2] + (double)v[iiz + 2]) - - + c3x - * ((double)v[iiz - incx3] - + (double)v[iiz + incx3]) - + c3y - * ((double)v[iiz - incy3] - + (double)v[iiz + incy3]) - + c3z * ((double)v[iiz - 3] + (double)v[iiz + 3]) - - + c4x - * ((double)v[iiz - incx4] - + (double)v[iiz + incx4]) - + c4y - * ((double)v[iiz - incy4] - + (double)v[iiz + incy4]) - + c4z * ((double)v[iiz - 4] + (double)v[iiz + 4])); + u[iiz] = (ScalarType)(c0 * (double)v[iiz] + + + c1x + * ((double)v[iiz - incx] + + (double)v[iiz + incx]) + + c1y + * ((double)v[iiz - incy] + + (double)v[iiz + incy]) + + c1z + * ((double)v[iiz - 1] + + (double)v[iiz + 1]) + + + c2x + * ((double)v[iiz - incx2] + + (double)v[iiz + incx2]) + + c2y + * ((double)v[iiz - incy2] + + (double)v[iiz + incy2]) + + c2z + * ((double)v[iiz - 2] + + (double)v[iiz + 2]) + + + c3x + * ((double)v[iiz - incx3] + + (double)v[iiz + incx3]) + + c3y + * ((double)v[iiz - incy3] + + (double)v[iiz + incy3]) + + c3z + * ((double)v[iiz - 3] + + (double)v[iiz + 3]) + + + c4x + * ((double)v[iiz - incx4] + + (double)v[iiz + incx4]) + + c4y + * ((double)v[iiz - incy4] + + (double)v[iiz + incy4]) + + c4z + * ((double)v[iiz - 4] + + (double)v[iiz + 4])); iiz++; } @@ -496,20 +516,26 @@ void FDkernelDel2_4th_Mehr(const Grid& grid, ScalarType* v, ScalarType* u, for (int iz = 0; iz < dim2; iz++) { - u0[iz] = (ScalarType)( - c0mehr4 * (double)v0[iz] - + czmehr4 * (double)(v0[iz - 1] + v0[iz + 1]) - + cymehr4 * (double)(vmy[iz] + vpy[iz]) - + cxmehr4 * (double)(vmx[iz] + vpx[iz]) - + cxzmehr4 - * (double)(vmx[iz - 1] + vmx[iz + 1] + vpx[iz - 1] - + vpx[iz + 1]) - + cyzmehr4 - * (double)(vmy[iz - 1] + vmy[iz + 1] + vpy[iz - 1] - + vpy[iz + 1]) - + cxymehr4 - * (double)(vmxmy[iz] + vpxmy[iz] + vmxpy[iz] - + vpxpy[iz])); + u0[iz] + = (ScalarType)(c0mehr4 * (double)v0[iz] + + czmehr4 + * (double)(v0[iz - 1] + v0[iz + 1]) + + cymehr4 * (double)(vmy[iz] + vpy[iz]) + + cxmehr4 * (double)(vmx[iz] + vpx[iz]) + + cxzmehr4 + * (double)(vmx[iz - 1] + + vmx[iz + 1] + + vpx[iz - 1] + + vpx[iz + 1]) + + cyzmehr4 + * (double)(vmy[iz - 1] + + vmy[iz + 1] + + vpy[iz - 1] + + vpy[iz + 1]) + + cxymehr4 + * (double)(vmxmy[iz] + vpxmy[iz] + + vmxpy[iz] + + vpxpy[iz])); } iiy += incy; @@ -568,11 +594,12 @@ void FDkernelRHS_4th_Mehr1(const Grid& grid, ScalarType* v, ScalarType* rhs, for (int iz = 0; iz < dim2; iz++) { - u0[iz] = (ScalarType)( - c0 * (double)v0[iz] - + c1 - * (double)(vmx[iz] + vpx[iz] + vmy[iz] + vpy[iz] - + v0[iz - 1] + v0[iz + 1])); + u0[iz] = (ScalarType)(c0 * (double)v0[iz] + + c1 + * (double)(vmx[iz] + vpx[iz] + + vmy[iz] + vpy[iz] + + v0[iz - 1] + + v0[iz + 1])); } iiy += incy; diff --git a/src/pb/FDoper.cc b/src/pb/FDoper.cc index 41530b1b..fbe1a1c9 100644 --- a/src/pb/FDoper.cc +++ b/src/pb/FDoper.cc @@ -216,10 +216,14 @@ void FDoper::del1_6th( for (int iz = 0; iz < dim2; iz++) { - u[iiz] = (T)( - e1 * ((double)v[iiz + incc] - (double)v[iiz - incc]) - + e2 * ((double)v[iiz - incc2] - (double)v[iiz + incc2]) - + e3 * ((double)v[iiz + incc3] - (double)v[iiz - incc3])); + u[iiz] + = (T)(e1 * ((double)v[iiz + incc] - (double)v[iiz - incc]) + + e2 + * ((double)v[iiz - incc2] + - (double)v[iiz + incc2]) + + e3 + * ((double)v[iiz + incc3] + - (double)v[iiz - incc3])); iiz++; } @@ -278,11 +282,17 @@ void FDoper::del1_8th( for (int iz = 0; iz < dim2; iz++) { - u[iiz] = (T)( - e1 * ((double)v[iiz + incc] - (double)v[iiz - incc]) - + e2 * ((double)v[iiz - incc2] - (double)v[iiz + incc2]) - + e3 * ((double)v[iiz + incc3] - (double)v[iiz - incc3]) - + e4 * ((double)v[iiz - incc4] - (double)v[iiz + incc4])); + u[iiz] + = (T)(e1 * ((double)v[iiz + incc] - (double)v[iiz - incc]) + + e2 + * ((double)v[iiz - incc2] + - (double)v[iiz + incc2]) + + e3 + * ((double)v[iiz + incc3] + - (double)v[iiz - incc3]) + + e4 + * ((double)v[iiz - incc4] + - (double)v[iiz + incc4])); iiz++; } @@ -609,24 +619,26 @@ void FDoper::rhs_4th_Mehr2(GridFunc& A, T* const u) const for (int iz = 0; iz < dim2; iz++) { - u0[iz] = (T)( - c0 * (double)v[iiz] - - + c1 - * (double)(v[iiz - incx_] + v[iiz + incx_] - + v[iiz - incy_] + v[iiz + incy_] - + v[iiz - 1] + v[iiz + 1]) - - + c2 - * (double)(v[iiz - incx_ - incy_] - + v[iiz + incx_ - incy_] - + v[iiz - incx_ + incy_] - + v[iiz + incx_ + incy_] - + v[iiz - incy_ - 1] + v[iiz - incy_ + 1] - + v[iiz + incy_ - 1] + v[iiz + incy_ + 1] - + v[iiz - incx_ - 1] + v[iiz - incx_ + 1] - + v[iiz + incx_ - 1] - + v[iiz + incx_ + 1])); + u0[iz] = (T)(c0 * (double)v[iiz] + + + c1 + * (double)(v[iiz - incx_] + v[iiz + incx_] + + v[iiz - incy_] + v[iiz + incy_] + + v[iiz - 1] + v[iiz + 1]) + + + c2 + * (double)(v[iiz - incx_ - incy_] + + v[iiz + incx_ - incy_] + + v[iiz - incx_ + incy_] + + v[iiz + incx_ + incy_] + + v[iiz - incy_ - 1] + + v[iiz - incy_ + 1] + + v[iiz + incy_ - 1] + + v[iiz + incy_ + 1] + + v[iiz - incx_ - 1] + + v[iiz - incx_ + 1] + + v[iiz + incx_ - 1] + + v[iiz + incx_ + 1])); iiz++; } diff --git a/src/pb/FDoperInterface.h b/src/pb/FDoperInterface.h index 07bad10d..5752a4bb 100644 --- a/src/pb/FDoperInterface.h +++ b/src/pb/FDoperInterface.h @@ -37,7 +37,7 @@ class FDoperInterface using memory_space_type = MemorySpace::Host; #endif - virtual ~FDoperInterface() {} + virtual ~FDoperInterface() { } static void printTimers(std::ostream& os) { diff --git a/src/pb/Grid.h b/src/pb/Grid.h index dbfe13bd..a4b28180 100644 --- a/src/pb/Grid.h +++ b/src/pb/Grid.h @@ -112,7 +112,7 @@ class Grid Vector3D closestGridPt(Vector3D coords) const; - ~Grid() {} + ~Grid() { } template void getSinCosFunctions(std::vector& sinx, std::vector& siny, diff --git a/src/pb/GridFunc.cc b/src/pb/GridFunc.cc index d80bad01..fe5ac4ed 100644 --- a/src/pb/GridFunc.cc +++ b/src/pb/GridFunc.cc @@ -359,7 +359,7 @@ GridFunc& GridFunc::operator=(const GridFunc& func) template void GridFunc::setZero() { - memset(uu_, 0, grid_.sizeg()*sizeof(T)); + memset(uu_, 0, grid_.sizeg() * sizeof(T)); updated_boundaries_ = true; } @@ -1513,7 +1513,7 @@ void GridFunc::initTrigo3d(const short bc[3], const int n[3]) double (*f2)(double); int init[3] = { 0, 0, 0 }; int end[3] = { dim_[0] + 2 * nghosts, dim_[1] + 2 * nghosts, - dim_[2] + 2 * nghosts }; + dim_[2] + 2 * nghosts }; if (bc[0] == 1) { f0 = &cos; @@ -3089,8 +3089,8 @@ void GridFunc::test_setBoundaryValues() assert(std::abs(my_dot - (grid_.gsize() - grid_.gdim(0) * grid_.gdim(1) - - grid_.gdim(0) * (grid_.gdim(2) - 1) - - (grid_.gdim(1) - 1) * (grid_.gdim(2) - 1))) + - grid_.gdim(0) * (grid_.gdim(2) - 1) + - (grid_.gdim(1) - 1) * (grid_.gdim(2) - 1))) < 1.e-8); if ((2 * (dim(0) >> 1) == dim(0)) && (2 * (dim(1) >> 1) == dim(1)) diff --git a/src/pb/GridFuncInterface.h b/src/pb/GridFuncInterface.h index b17d2307..b9e908da 100644 --- a/src/pb/GridFuncInterface.h +++ b/src/pb/GridFuncInterface.h @@ -30,7 +30,7 @@ class GridFuncInterface static Timer finishExchangeEastWest_tm_; public: - virtual ~GridFuncInterface() {} + virtual ~GridFuncInterface() { } static void printTimers(std::ostream& os) { diff --git a/src/pb/Lap.h b/src/pb/Lap.h index 460e36c1..af766ba3 100644 --- a/src/pb/Lap.h +++ b/src/pb/Lap.h @@ -26,7 +26,7 @@ class Lap : public FDoper std::string name_; public: - Lap(const Grid& mygrid) : FDoper(mygrid) {} + Lap(const Grid& mygrid) : FDoper(mygrid) { } Lap& operator=(const Lap& v) = default; @@ -50,7 +50,7 @@ class Lap : public FDoper virtual void setLowerOrderGrid(void) = 0; virtual double jacobiFactor(void) const = 0; - ~Lap() override {} + ~Lap() override { } }; } // namespace pb diff --git a/src/pb/Laph4.h b/src/pb/Laph4.h index b433a39e..c846526a 100644 --- a/src/pb/Laph4.h +++ b/src/pb/Laph4.h @@ -34,7 +34,7 @@ class Laph4 : public Lap //-1/12 16/12 -30/12 16/12 -1/12 diagEl_ = 2.5 * (Lap::inv_h2(0) + Lap::inv_h2(1) - + Lap::inv_h2(2)); // 2.5 = 30./12. + + Lap::inv_h2(2)); // 2.5 = 30./12. invDiagEl_ = 1. / diagEl_; Lap::name_ = "Classical 4th order"; diff --git a/src/pb/MGkernels.cc b/src/pb/MGkernels.cc index 253a57fb..d500ecc2 100644 --- a/src/pb/MGkernels.cc +++ b/src/pb/MGkernels.cc @@ -114,7 +114,7 @@ void MGkernelExtend3D(ScalarType* coarse_data, const Grid& coarse_grid, int izf = ify + iz; fine_data[izf] = 0.5 * (fine_data[izf + incy_fine] - + fine_data[izf - incy_fine]); + + fine_data[izf - incy_fine]); assert(izf < static_cast(nfunc * fine_grid.sizeg())); } @@ -124,9 +124,9 @@ void MGkernelExtend3D(ScalarType* coarse_data, const Grid& coarse_grid, int izf = ify + iz; fine_data[izf] = 0.25 * (fine_data[izf + 1 + incy_fine] - + fine_data[izf + 1 - incy_fine] - + fine_data[izf - 1 + incy_fine] - + fine_data[izf - 1 - incy_fine]); + + fine_data[izf + 1 - incy_fine] + + fine_data[izf - 1 + incy_fine] + + fine_data[izf - 1 - incy_fine]); assert(izf < static_cast(nfunc * fine_grid.sizeg())); } } @@ -146,9 +146,9 @@ void MGkernelExtend3D(ScalarType* coarse_data, const Grid& coarse_grid, int izf = ify + iz; fine_data[izf] = 0.25 * (fine_data[izf + incx_fine + 1] - + fine_data[izf + incx_fine - 1] - + fine_data[izf - incx_fine + 1] - + fine_data[izf - incx_fine - 1]); + + fine_data[izf + incx_fine - 1] + + fine_data[izf - incx_fine + 1] + + fine_data[izf - incx_fine - 1]); } for (int iz = nghosts_fine; iz < dimz + nghosts_fine; @@ -157,7 +157,7 @@ void MGkernelExtend3D(ScalarType* coarse_data, const Grid& coarse_grid, int izf = ify + iz; fine_data[izf] = 0.5 * (fine_data[izf + incx_fine] - + fine_data[izf - incx_fine]); + + fine_data[izf - incx_fine]); } } @@ -173,13 +173,13 @@ void MGkernelExtend3D(ScalarType* coarse_data, const Grid& coarse_grid, fine_data[izf] = 0.125 * (fine_data[izf + incx_fine + incy_fine + 1] - + fine_data[izf + incx_fine + incy_fine - 1] - + fine_data[izf + incx_fine - incy_fine + 1] - + fine_data[izf + incx_fine - incy_fine - 1] - + fine_data[izf - incx_fine + incy_fine + 1] - + fine_data[izf - incx_fine + incy_fine - 1] - + fine_data[izf - incx_fine - incy_fine + 1] - + fine_data[izf - incx_fine - incy_fine - 1]); + + fine_data[izf + incx_fine + incy_fine - 1] + + fine_data[izf + incx_fine - incy_fine + 1] + + fine_data[izf + incx_fine - incy_fine - 1] + + fine_data[izf - incx_fine + incy_fine + 1] + + fine_data[izf - incx_fine + incy_fine - 1] + + fine_data[izf - incx_fine - incy_fine + 1] + + fine_data[izf - incx_fine - incy_fine - 1]); } for (int iz = nghosts_fine; iz < dimz + nghosts_fine; @@ -189,9 +189,9 @@ void MGkernelExtend3D(ScalarType* coarse_data, const Grid& coarse_grid, fine_data[izf] = 0.25 * (fine_data[izf + incx_fine + incy_fine] - + fine_data[izf + incx_fine - incy_fine] - + fine_data[izf - incx_fine + incy_fine] - + fine_data[izf - incx_fine - incy_fine]); + + fine_data[izf + incx_fine - incy_fine] + + fine_data[izf - incx_fine + incy_fine] + + fine_data[izf - incx_fine - incy_fine]); } } } @@ -265,9 +265,10 @@ void MGkernelRestrict3D(ScalarType* fine_data, const Grid& fine_grid, + (double)umxmy[twoiz] + (double)upxmy[twoiz] + (double)umxpy[twoiz] + (double)upxpy[twoiz]; - coarse_data[iyc + iz + nghosts_coarse] = (ScalarType)( - inv64 - * (8. * u0[twoiz] + 4. * face + 2. * edge + corner)); + coarse_data[iyc + iz + nghosts_coarse] + = (ScalarType)(inv64 + * (8. * u0[twoiz] + 4. * face + 2. * edge + + corner)); } iyf += 2 * incy_fine; diff --git a/src/pb/Oper.h b/src/pb/Oper.h index 2d126a19..42dd9889 100644 --- a/src/pb/Oper.h +++ b/src/pb/Oper.h @@ -19,7 +19,7 @@ class Oper { public: - virtual ~Oper() {} + virtual ~Oper() { } }; } // namespace pb diff --git a/src/pb/PBh4MP.h b/src/pb/PBh4MP.h index 38537ccc..0f06f553 100644 --- a/src/pb/PBh4MP.h +++ b/src/pb/PBh4MP.h @@ -23,13 +23,13 @@ class PBh4MP : public PBh4M public: // constructor PBh4MP(const Grid& mygrid, DielFunc& myepsilon) - : PBh4M(mygrid, myepsilon){}; + : PBh4M(mygrid, myepsilon) {}; PBh4MP(const Grid& mygrid, const double e0, const double rho0, const double drho0) - : PBh4M(mygrid, e0, rho0, drho0){}; + : PBh4M(mygrid, e0, rho0, drho0) {}; PBh4MP(const Grid& mygrid, DielFunc& myepsilon, GridFunc& pp) - : PBh4M(mygrid, myepsilon, pp){}; + : PBh4M(mygrid, myepsilon, pp) {}; // construct a coarse grid operator PBh4MP coarseOp(const Grid& mygrid) diff --git a/src/pb/Solver.h b/src/pb/Solver.h index f82b3b6a..6754bfda 100644 --- a/src/pb/Solver.h +++ b/src/pb/Solver.h @@ -39,7 +39,7 @@ class Solver const bool gather_coarse_level = true) = 0; - virtual ~Solver() {} + virtual ~Solver() { } virtual short getNbSweeps() const = 0; virtual double getFinalResidual() const = 0; diff --git a/src/pb/SolverPB.h b/src/pb/SolverPB.h index 99eae0e3..2da00517 100644 --- a/src/pb/SolverPB.h +++ b/src/pb/SolverPB.h @@ -68,7 +68,7 @@ class SolverPB : public Solver GridFunc& gf_rhod, GridFunc& gf_vks); bool solve(GridFunc& gf_phi, const GridFunc& gf_rhs) override; - ~SolverPB() override{}; + ~SolverPB() override {}; short getNbSweeps() const override { return nb_sweeps_; } double getFinalResidual() const override { return final_residual_; } diff --git a/src/radial/RadialInter.h b/src/radial/RadialInter.h index 28e67f0c..3d7284e6 100644 --- a/src/radial/RadialInter.h +++ b/src/radial/RadialInter.h @@ -16,16 +16,16 @@ class RadialInter : public RadialMeshFunction { private: public: - RadialInter(const std::vector& x) : RadialMeshFunction(x) {} + RadialInter(const std::vector& x) : RadialMeshFunction(x) { } - RadialInter() : RadialMeshFunction() {} + RadialInter() : RadialMeshFunction() { } RadialInter(std::vector& x, std::vector>& y) : RadialMeshFunction(x, y) { } - ~RadialInter() override {} + ~RadialInter() override { } double linint(const double x, const int j = 0) const; double cubint(const double x, const int j = 0) const; diff --git a/src/radial/RadialMeshFunction.h b/src/radial/RadialMeshFunction.h index 06464e65..b8f21bfd 100644 --- a/src/radial/RadialMeshFunction.h +++ b/src/radial/RadialMeshFunction.h @@ -36,7 +36,7 @@ class RadialMeshFunction invdr_ = -1.; }; - virtual ~RadialMeshFunction(){}; + virtual ~RadialMeshFunction() {}; RadialMeshFunction(const std::vector& x) { diff --git a/src/readInput.cc b/src/readInput.cc index bcb823c4..6de74058 100644 --- a/src/readInput.cc +++ b/src/readInput.cc @@ -21,7 +21,7 @@ #include -//#define DEBUG 1 +// #define DEBUG 1 template int MGmol::readLRsFromInput(std::ifstream* tfile) diff --git a/src/read_config.cc b/src/read_config.cc index 2444fd6a..96ac1f74 100644 --- a/src/read_config.cc +++ b/src/read_config.cc @@ -7,6 +7,7 @@ // This file is part of MGmol. For details, see https://github.com/llnl/mgmol. // Please also read this link https://github.com/llnl/mgmol/LICENSE +#include "MGmol_prototypes.h" #include #include #include @@ -14,7 +15,6 @@ #include #include #include -#include "MGmol_prototypes.h" #include namespace po = boost::program_options; @@ -118,8 +118,7 @@ int read_config(int argc, char** argv, po::variables_map& vm, "Compute MLWF (apply rotation) in quench")( "Quench.num_lin_iterations", po::value()->default_value(0), "Number of iterations without potential update in quench")( - "Preconditioner.num_levels", - po::value()->default_value(2), + "Preconditioner.num_levels", po::value()->default_value(2), "Number of levels for MG preconditioner")( "Preconditioner.npresmoothing", po::value()->default_value(2), @@ -127,8 +126,7 @@ int read_config(int argc, char** argv, po::variables_map& vm, "Preconditioner.npostsmoothing", po::value()->default_value(2), "Number of postsmoothing steps i preconditioner")( - "Preconditioner.precision", - po::value()->default_value(32), + "Preconditioner.precision", po::value()->default_value(32), "Precision for MG preconditioner")("Quench.spread_penalty_damping", po::value()->default_value(0.), "Spread penalty damping factor")("Quench.spread_penalty_target", @@ -257,8 +255,9 @@ int read_config(int argc, char** argv, po::variables_map& vm, "continuum solvent: beta")("Poisson.FDtype", po::value()->default_value("Mehrstellen"), "FDtype")("Poisson.nu1", po::value()->default_value(1), - "MG pre-smoothing sweeps")("Poisson.nu2", po::value()->default_value(1), - "MG post-smoothing sweeps")("Poisson.max_steps", po::value()->default_value(20), + "MG pre-smoothing sweeps")("Poisson.nu2", + po::value()->default_value(1), "MG post-smoothing sweeps")( + "Poisson.max_steps", po::value()->default_value(20), "max. nb. steps Poisson solver")("Poisson.max_steps_initial", po::value()->default_value(20), "max. nb. steps Poisson solver in first solve")( @@ -423,32 +422,35 @@ int read_config(int argc, char** argv, po::variables_map& vm, } #ifdef MGMOL_HAS_LIBROM -void setupROMConfigOption(po::options_description &rom_cfg) +void setupROMConfigOption(po::options_description& rom_cfg) { - rom_cfg.add_options() - ("ROM.stage", po::value()->default_value("none"), - "ROM workflow stage: offline; build; online; none.") - ("ROM.offline.restart_filefmt", po::value()->default_value(""), - "File name format to read for snapshots.") - ("ROM.offline.restart_min_idx", po::value()->default_value(-1), - "Minimum index for snapshot file format.") - ("ROM.offline.restart_max_idx", po::value()->default_value(-1), - "Maximum index for snapshot file format.") - ("ROM.offline.basis_file", po::value()->default_value(""), - "File name for libROM snapshot/POD matrices.") - ("ROM.offline.save_librom_snapshot", po::value()->default_value(false), - "Save libROM snapshot file at FOM simulation.") - ("ROM.offline.librom_snapshot_freq", po::value()->default_value(-1), - "Frequency of saving libROM snapshot file at FOM simulation.") - ("ROM.offline.variable", po::value()->default_value(""), - "FOM variable to perform POD: either orbitals or potential.") - ("ROM.basis.compare_md", po::value()->default_value(false), - "Compare MD or single-step force.") - ("ROM.basis.number_of_orbital_basis", po::value()->default_value(-1), - "Number of orbital POD basis.") - ("ROM.basis.number_of_potential_basis", po::value()->default_value(-1), - "Number of potential POD basis to build Hartree potential ROM operator.") - ("ROM.potential_rom_file", po::value()->default_value(""), - "File name to save/load potential ROM operators."); + rom_cfg.add_options()("ROM.stage", + po::value()->default_value("none"), + "ROM workflow stage: offline; build; online; none.")( + "ROM.offline.restart_filefmt", + po::value()->default_value(""), + "File name format to read for snapshots.")( + "ROM.offline.restart_min_idx", po::value()->default_value(-1), + "Minimum index for snapshot file format.")( + "ROM.offline.restart_max_idx", po::value()->default_value(-1), + "Maximum index for snapshot file format.")("ROM.offline.basis_file", + po::value()->default_value(""), + "File name for libROM snapshot/POD matrices.")( + "ROM.offline.save_librom_snapshot", + po::value()->default_value(false), + "Save libROM snapshot file at FOM simulation.")( + "ROM.offline.librom_snapshot_freq", po::value()->default_value(-1), + "Frequency of saving libROM snapshot file at FOM simulation.")( + "ROM.offline.variable", po::value()->default_value(""), + "FOM variable to perform POD: either orbitals or potential.")( + "ROM.basis.compare_md", po::value()->default_value(false), + "Compare MD or single-step force.")("ROM.basis.number_of_orbital_basis", + po::value()->default_value(-1), + "Number of orbital POD basis.")("ROM.basis.number_of_potential_basis", + po::value()->default_value(-1), + "Number of potential POD basis to build Hartree potential ROM " + "operator.")("ROM.potential_rom_file", + po::value()->default_value(""), + "File name to save/load potential ROM operators."); } #endif diff --git a/src/rom.cc b/src/rom.cc index 863d76e9..00bc0ee0 100644 --- a/src/rom.cc +++ b/src/rom.cc @@ -10,23 +10,24 @@ #include "mgmol_config.h" #ifdef MGMOL_HAS_LIBROM -#include "LocGridOrbitals.h" #include "ExtendedGridOrbitals.h" +#include "LocGridOrbitals.h" #include "MGmol.h" #include "librom.h" -#include -#include #include +#include +#include #include // Save the wavefunction snapshots template -int MGmol::save_orbital_snapshot(std::string file_path, OrbitalsType& orbitals) +int MGmol::save_orbital_snapshot( + std::string file_path, OrbitalsType& orbitals) { std::string snapshot_filename = file_path; - struct stat s; + struct stat s; if (stat(file_path.c_str(), &s) == 0) { if (s.st_mode & S_IFDIR) @@ -39,16 +40,19 @@ int MGmol::save_orbital_snapshot(std::string file_path, OrbitalsTy } else { - std::cout << file_path << " exists but is not a directory or a file." << std::endl; + std::cout << file_path + << " exists but is not a directory or a file." + << std::endl; return 1; } } - const int dim = orbitals.getLocNumpt(); + const int dim = orbitals.getLocNumpt(); const int totalSamples = orbitals.chromatic_number(); CAROM::Options svd_options(dim, totalSamples, 1); - CAROM::BasisGenerator basis_generator(svd_options, false, snapshot_filename); + CAROM::BasisGenerator basis_generator( + svd_options, false, snapshot_filename); for (int i = 0; i < totalSamples; ++i) basis_generator.takeSample(orbitals.getPsi(i)); @@ -59,9 +63,10 @@ int MGmol::save_orbital_snapshot(std::string file_path, OrbitalsTy } template -void MGmol::project_orbital(std::string file_path, int rdim, OrbitalsType& orbitals) +void MGmol::project_orbital( + std::string file_path, int rdim, OrbitalsType& orbitals) { - const int dim = orbitals.getLocNumpt(); + const int dim = orbitals.getLocNumpt(); const int totalSamples = orbitals.chromatic_number(); CAROM::Options svd_options(dim, totalSamples, 1); @@ -69,17 +74,21 @@ void MGmol::project_orbital(std::string file_path, int rdim, Orbit for (int i = 0; i < totalSamples; ++i) basis_generator.takeSample(orbitals.getPsi(i)); - const CAROM::Matrix* orbital_snapshots = basis_generator.getSnapshotMatrix(); + const CAROM::Matrix* orbital_snapshots + = basis_generator.getSnapshotMatrix(); CAROM::BasisReader reader(file_path); CAROM::Matrix* orbital_basis = reader.getSpatialBasis(rdim); - CAROM::Matrix* proj_orbital_coeff = orbital_basis->transposeMult(orbital_snapshots); - CAROM::Matrix* proj_orbital_snapshots = orbital_basis->mult(proj_orbital_coeff); - + CAROM::Matrix* proj_orbital_coeff + = orbital_basis->transposeMult(orbital_snapshots); + CAROM::Matrix* proj_orbital_snapshots + = orbital_basis->mult(proj_orbital_coeff); + Control& ct = *(Control::instance()); - Mesh* mesh = Mesh::instance(); - pb::GridFunc gf_psi(mesh->grid(), ct.bcWF[0], ct.bcWF[1], ct.bcWF[2]); + Mesh* mesh = Mesh::instance(); + pb::GridFunc gf_psi( + mesh->grid(), ct.bcWF[0], ct.bcWF[1], ct.bcWF[2]); CAROM::Vector snapshot, proj_snapshot; for (int i = 0; i < totalSamples; ++i) { @@ -88,23 +97,23 @@ void MGmol::project_orbital(std::string file_path, int rdim, Orbit gf_psi.assign(proj_snapshot.getData()); orbitals.setPsi(gf_psi, i); snapshot -= proj_snapshot; - std::cout << "Error for orbital " << i << " = " << snapshot.norm() << std::endl; + std::cout << "Error for orbital " << i << " = " << snapshot.norm() + << std::endl; } } -//template -//void ExtendedGridOrbitals::set(std::string file_path, int rdim) +// template +// void ExtendedGridOrbitals::set(std::string file_path, int rdim) //{ -// const int dim = getLocNumpt(); +// const int dim = getLocNumpt(); // CAROM::BasisReader reader(file_path); // CAROM::Matrix* orbital_basis = reader.getSpatialBasis(rdim); // Control& ct = *(Control::instance()); // Mesh* mymesh = Mesh::instance(); -// pb::GridFunc gf_psi(mymesh->grid(), ct.bcWF[0], ct.bcWF[1], ct.bcWF[2]); -// CAROM::Vector psi; -// for (int i = 0; i < rdim; ++i) +// pb::GridFunc gf_psi(mymesh->grid(), ct.bcWF[0], ct.bcWF[1], +// ct.bcWF[2]); CAROM::Vector psi; for (int i = 0; i < rdim; ++i) // { // orbital_basis->getColumn(i, psi); // gf_psi.assign(psi.getData()); @@ -115,4 +124,4 @@ void MGmol::project_orbital(std::string file_path, int rdim, Orbit template class MGmol>; template class MGmol>; -#endif // MGMOL_HAS_LIBROM +#endif // MGMOL_HAS_LIBROM diff --git a/src/rom_Control.h b/src/rom_Control.h index 2d143d84..0e852522 100644 --- a/src/rom_Control.h +++ b/src/rom_Control.h @@ -20,7 +20,7 @@ enum class ROMStage { OFFLINE, ONLINE, - RESTORE, // TODO(kevin): what stage is this? + RESTORE, // TODO(kevin): what stage is this? BUILD, ONLINE_PINNED_H2O_3DOF, TEST_ORBITAL, @@ -44,20 +44,20 @@ struct ROMPrivateOptions ROMStage rom_stage = ROMStage::UNSUPPORTED; std::string restart_file_fmt = ""; - int restart_file_minidx = -1; - int restart_file_maxidx = -1; - std::string basis_file = ""; - ROMVariable variable = ROMVariable::NONE; + int restart_file_minidx = -1; + int restart_file_maxidx = -1; + std::string basis_file = ""; + ROMVariable variable = ROMVariable::NONE; /* save librom orbital snapshot matrix at FOM simulation. */ bool save_librom_snapshot = false; - int librom_snapshot_freq = -1; + int librom_snapshot_freq = -1; /* options for ROM building */ - bool compare_md = false; - int num_orbbasis = -1; - int num_potbasis = -1; + bool compare_md = false; + int num_orbbasis = -1; + int num_potbasis = -1; std::string pot_rom_file = ""; }; -#endif // ROM_CONTROL_H +#endif // ROM_CONTROL_H diff --git a/src/rom_main.cc b/src/rom_main.cc index d29419b9..f9014bee 100644 --- a/src/rom_main.cc +++ b/src/rom_main.cc @@ -88,72 +88,74 @@ int main(int argc, char** argv) /* Get ROM driver mode */ ROMPrivateOptions rom_options = ct.getROMOptions(); - ROMStage rom_stage = rom_options.rom_stage; + ROMStage rom_stage = rom_options.rom_stage; // Enter main scope { MGmolInterface* mgmol; if (ct.isLocMode()) - mgmol = new MGmol>(global_comm, *MPIdata::sout, - input_filename, lrs_filename, constraints_filename); + mgmol = new MGmol>(global_comm, + *MPIdata::sout, input_filename, lrs_filename, + constraints_filename); else - mgmol = new MGmol>(global_comm, *MPIdata::sout, - input_filename, lrs_filename, constraints_filename); + mgmol = new MGmol>(global_comm, + *MPIdata::sout, input_filename, lrs_filename, + constraints_filename); mgmol->setup(); switch (rom_stage) { - case (ROMStage::OFFLINE): - if (ct.isLocMode()) - readRestartFiles>(mgmol); - else - readRestartFiles>(mgmol); - break; - - case (ROMStage::BUILD): - if (ct.isLocMode()) - buildROMPoissonOperator>(mgmol); - else - buildROMPoissonOperator>(mgmol); - break; - - case (ROMStage::ONLINE_POISSON): - if (ct.isLocMode()) - runPoissonROM>(mgmol); - else - runPoissonROM>(mgmol); - break; - - case (ROMStage::TEST_POISSON): - if (ct.isLocMode()) - testROMPoissonOperator>(mgmol); - else - testROMPoissonOperator>(mgmol); - break; - - case (ROMStage::TEST_RHO): - if (ct.isLocMode()) - testROMRhoOperator>(mgmol); - else - testROMRhoOperator>(mgmol); - - case (ROMStage::TEST_ION): - if (ct.isLocMode()) - testROMIonDensity>(mgmol); - else - testROMIonDensity>(mgmol); - - break; - - default: - std::cerr << "rom_main error: Unknown ROM stage" << std::endl; - MPI_Abort(MPI_COMM_WORLD, 0); - break; + case (ROMStage::OFFLINE): + if (ct.isLocMode()) + readRestartFiles>(mgmol); + else + readRestartFiles>(mgmol); + break; + + case (ROMStage::BUILD): + if (ct.isLocMode()) + buildROMPoissonOperator>(mgmol); + else + buildROMPoissonOperator>( + mgmol); + break; + + case (ROMStage::ONLINE_POISSON): + if (ct.isLocMode()) + runPoissonROM>(mgmol); + else + runPoissonROM>(mgmol); + break; + + case (ROMStage::TEST_POISSON): + if (ct.isLocMode()) + testROMPoissonOperator>(mgmol); + else + testROMPoissonOperator>( + mgmol); + break; + + case (ROMStage::TEST_RHO): + if (ct.isLocMode()) + testROMRhoOperator>(mgmol); + else + testROMRhoOperator>(mgmol); + + case (ROMStage::TEST_ION): + if (ct.isLocMode()) + testROMIonDensity>(mgmol); + else + testROMIonDensity>(mgmol); + + break; + + default: + std::cerr << "rom_main error: Unknown ROM stage" << std::endl; + MPI_Abort(MPI_COMM_WORLD, 0); + break; } - - delete mgmol; } // close main scope diff --git a/src/rom_workflows.cc b/src/rom_workflows.cc index 00ffefa5..73584eec 100644 --- a/src/rom_workflows.cc +++ b/src/rom_workflows.cc @@ -10,16 +10,18 @@ #include "rom_workflows.h" #include "Electrostatic.h" #include -#include #include +#include -namespace CAROM { +namespace CAROM +{ /* This is not implemented in libROM yet. A temporary placeholder until it is merged into libROM. */ -void mult(const Matrix &A, const std::vector &sample_row_A, const Matrix &B, Matrix &AB) +void mult(const Matrix& A, const std::vector& sample_row_A, + const Matrix& B, Matrix& AB) { CAROM_VERIFY(!B.distributed()); CAROM_VERIFY(A.distributed() == AB.distributed()); @@ -32,13 +34,16 @@ void mult(const Matrix &A, const std::vector &sample_row_A, const Matrix &B // Do the multiplication. const int Acol = A.numColumns(); - for (int r = 0; r < num_rows; ++r) { - double *d_Arow = A.getData() + sample_row_A[r] * Acol; - for (int c = 0; c < num_cols; ++c) { - double *d_Aptr = d_Arow; - double *d_Bptr = B.getData() + c; + for (int r = 0; r < num_rows; ++r) + { + double* d_Arow = A.getData() + sample_row_A[r] * Acol; + for (int c = 0; c < num_cols; ++c) + { + double* d_Aptr = d_Arow; + double* d_Bptr = B.getData() + c; double result_val = 0.0; - for (int ac = 0; ac < Acol; ++ac, d_Aptr++) { + for (int ac = 0; ac < Acol; ++ac, d_Aptr++) + { result_val += (*d_Aptr) * (*d_Bptr); d_Bptr += num_cols; } @@ -47,21 +52,26 @@ void mult(const Matrix &A, const std::vector &sample_row_A, const Matrix &B } } -} +} -template -std::string string_format( const std::string& format, Args ... args ) +template +std::string string_format(const std::string& format, Args... args) { - int size_s = std::snprintf( nullptr, 0, format.c_str(), args ... ) + 1; // Extra space for '\0' - if( size_s <= 0 ){ throw std::runtime_error( "Error during formatting." ); } - auto size = static_cast( size_s ); - std::unique_ptr buf( new char[ size ] ); - std::snprintf( buf.get(), size, format.c_str(), args ... ); - return std::string( buf.get(), buf.get() + size - 1 ); // We don't want the '\0' inside + int size_s = std::snprintf(nullptr, 0, format.c_str(), args...) + + 1; // Extra space for '\0' + if (size_s <= 0) + { + throw std::runtime_error("Error during formatting."); + } + auto size = static_cast(size_s); + std::unique_ptr buf(new char[size]); + std::snprintf(buf.get(), size, format.c_str(), args...); + return std::string( + buf.get(), buf.get() + size - 1); // We don't want the '\0' inside } template -void readRestartFiles(MGmolInterface *mgmol_) +void readRestartFiles(MGmolInterface* mgmol_) { Control& ct = *(Control::instance()); Mesh* mymesh = Mesh::instance(); @@ -74,43 +84,44 @@ void readRestartFiles(MGmolInterface *mgmol_) /* number of restart files, start/end indices */ assert(rom_options.restart_file_minidx >= 0); assert(rom_options.restart_file_maxidx >= 0); - const int minidx = rom_options.restart_file_minidx; - const int maxidx = rom_options.restart_file_maxidx; + const int minidx = rom_options.restart_file_minidx; + const int maxidx = rom_options.restart_file_maxidx; const int num_restart = maxidx - minidx + 1; - MGmol *mgmol = static_cast *>(mgmol_); - OrbitalsType *orbitals = mgmol->getOrbitals(); - Potentials& pot = mgmol->getHamiltonian()->potential(); + MGmol* mgmol = static_cast*>(mgmol_); + OrbitalsType* orbitals = mgmol->getOrbitals(); + Potentials& pot = mgmol->getHamiltonian()->potential(); std::string filename; /* Determine basis prefix, dimension, and sample size */ std::string basis_prefix = rom_options.basis_file; int dim; - int totalSamples = num_restart; + int totalSamples = num_restart; const int chrom_num = orbitals->chromatic_number(); switch (rom_var) { - case ROMVariable::ORBITALS: - basis_prefix += "_orbitals"; - dim = orbitals->getLocNumpt(); - /* if orbitals, each sample have chromatic number of wave functions */ - totalSamples *= orbitals->chromatic_number(); - break; - - case ROMVariable::POTENTIAL: - basis_prefix += "_potential"; - dim = pot.size(); - break; - - default: - ct.global_exit(); - break; + case ROMVariable::ORBITALS: + basis_prefix += "_orbitals"; + dim = orbitals->getLocNumpt(); + /* if orbitals, each sample have chromatic number of wave functions + */ + totalSamples *= orbitals->chromatic_number(); + break; + + case ROMVariable::POTENTIAL: + basis_prefix += "_potential"; + dim = pot.size(); + break; + + default: + ct.global_exit(); + break; } /* Initialize libROM classes */ CAROM::Options svd_options(dim, totalSamples, 1); CAROM::BasisGenerator basis_generator(svd_options, false, basis_prefix); - + /* Collect the restart files */ for (int k = minidx; k <= maxidx; k++) { @@ -121,17 +132,18 @@ void readRestartFiles(MGmolInterface *mgmol_) switch (rom_var) { - case ROMVariable::ORBITALS: - for (int i = 0; i < chrom_num; ++i) - basis_generator.takeSample(orbitals->getPsi(i)); - break; - - case ROMVariable::POTENTIAL: - basis_prefix += "_potential"; - /* we save hartree potential */ - // TODO: consider create overloaded takeSample with const - basis_generator.takeSample(const_cast(pot.vh_rho().data())); - break; + case ROMVariable::ORBITALS: + for (int i = 0; i < chrom_num; ++i) + basis_generator.takeSample(orbitals->getPsi(i)); + break; + + case ROMVariable::POTENTIAL: + basis_prefix += "_potential"; + /* we save hartree potential */ + // TODO: consider create overloaded takeSample with const + basis_generator.takeSample( + const_cast(pot.vh_rho().data())); + break; } } basis_generator.writeSnapshot(); @@ -139,7 +151,7 @@ void readRestartFiles(MGmolInterface *mgmol_) } template -void buildROMPoissonOperator(MGmolInterface *mgmol_) +void buildROMPoissonOperator(MGmolInterface* mgmol_) { Control& ct = *(Control::instance()); Mesh* mymesh = Mesh::instance(); @@ -153,22 +165,24 @@ void buildROMPoissonOperator(MGmolInterface *mgmol_) ROMVariable rom_var = rom_options.variable; if (rom_var != ROMVariable::POTENTIAL) { - std::cerr << "buildROMPoissonOperator error: ROM variable must be POTENTIAL to run this stage!\n" << std::endl; + std::cerr << "buildROMPoissonOperator error: ROM variable must be " + "POTENTIAL to run this stage!\n" + << std::endl; MPI_Abort(MPI_COMM_WORLD, 0); } /* Load Hartree potential basis matrix */ - std::string basis_file = rom_options.basis_file; + std::string basis_file = rom_options.basis_file; const int num_pot_basis = rom_options.num_potbasis; CAROM::BasisReader basis_reader(basis_file); - CAROM::Matrix *pot_basis = basis_reader.getSpatialBasis(num_pot_basis); + CAROM::Matrix* pot_basis = basis_reader.getSpatialBasis(num_pot_basis); /* Load PoissonSolver pointer */ - MGmol *mgmol = static_cast *>(mgmol_); - Poisson *poisson = mgmol->electrostat_->getPoissonSolver(); + MGmol* mgmol = static_cast*>(mgmol_); + Poisson* poisson = mgmol->electrostat_->getPoissonSolver(); /* GridFunc initialization inputs */ - const pb::Grid &grid(poisson->vh().grid()); + const pb::Grid& grid(poisson->vh().grid()); short bc[3]; for (int d = 0; d < 3; d++) bc[d] = poisson->vh().bc(d); @@ -183,7 +197,7 @@ void buildROMPoissonOperator(MGmolInterface *mgmol_) for (int c = 0; c < num_pot_basis; c++) { /* copy c-th column librom vector to GridFunc gf_col */ - CAROM::Vector *col = pot_basis->getColumn(c); + CAROM::Vector* col = pot_basis->getColumn(c); col_gf.assign(col->getData(), 'd'); /* apply Laplace operator */ @@ -201,20 +215,21 @@ void buildROMPoissonOperator(MGmolInterface *mgmol_) pot_rom(r, c) = rom_col(r); delete col; - } // for (int c = 0; c < num_pot_basis; c++) + } // for (int c = 0; c < num_pot_basis; c++) /* DEIM hyperreduction */ CAROM::Matrix pot_rhs_rom(num_pot_basis, num_pot_basis, false); - std::vector global_sampled_row(num_pot_basis), sampled_rows_per_proc(nprocs); + std::vector global_sampled_row(num_pot_basis), + sampled_rows_per_proc(nprocs); DEIM(pot_basis, num_pot_basis, global_sampled_row, sampled_rows_per_proc, - pot_rhs_rom, rank, nprocs); + pot_rhs_rom, rank, nprocs); /* ROM rescaleTotalCharge operator */ CAROM::Vector fom_ones(pot_basis->numRows(), true); CAROM::Vector rom_ones(num_pot_basis, false); fom_ones = mymesh->grid().vel(); // volume element pot_basis->transposeMult(fom_ones, rom_ones); - + /* Save ROM operator */ // write the file from PE0 only if (MPIdata::onpe0) @@ -224,27 +239,27 @@ void buildROMPoissonOperator(MGmolInterface *mgmol_) h5_helper.create(rom_oper); h5_helper.putInteger("number_of_potential_basis", num_pot_basis); h5_helper.putDoubleArray("potential_rom_operator", pot_rom.getData(), - num_pot_basis * num_pot_basis, false); + num_pot_basis * num_pot_basis, false); /* save the inverse as well */ pot_rom.inverse(); h5_helper.putDoubleArray("potential_rom_inverse", pot_rom.getData(), - num_pot_basis * num_pot_basis, false); + num_pot_basis * num_pot_basis, false); /* save right-hand side hyper-reduction operator */ - h5_helper.putDoubleArray("potential_rhs_rom_inverse", pot_rhs_rom.getData(), - num_pot_basis * num_pot_basis, false); + h5_helper.putDoubleArray("potential_rhs_rom_inverse", + pot_rhs_rom.getData(), num_pot_basis * num_pot_basis, false); /* save right-hand side rescaling operator */ - h5_helper.putDoubleArray("potential_rhs_rescaler", rom_ones.getData(), - num_pot_basis, false); + h5_helper.putDoubleArray( + "potential_rhs_rescaler", rom_ones.getData(), num_pot_basis, false); h5_helper.close(); } } template -void runPoissonROM(MGmolInterface *mgmol_) +void runPoissonROM(MGmolInterface* mgmol_) { Control& ct = *(Control::instance()); Mesh* mymesh = Mesh::instance(); @@ -258,22 +273,24 @@ void runPoissonROM(MGmolInterface *mgmol_) ROMVariable rom_var = rom_options.variable; if (rom_var != ROMVariable::POTENTIAL) { - std::cerr << "runPoissonROM error: ROM variable must be POTENTIAL to run this stage!\n" << std::endl; + std::cerr << "runPoissonROM error: ROM variable must be POTENTIAL to " + "run this stage!\n" + << std::endl; MPI_Abort(MPI_COMM_WORLD, 0); } /* Load Hartree potential basis matrix */ - std::string basis_file = rom_options.basis_file; + std::string basis_file = rom_options.basis_file; const int num_pot_basis = rom_options.num_potbasis; CAROM::BasisReader basis_reader(basis_file); - CAROM::Matrix *pot_basis = basis_reader.getSpatialBasis(num_pot_basis); + CAROM::Matrix* pot_basis = basis_reader.getSpatialBasis(num_pot_basis); /* initialize rom operator variables */ CAROM::Matrix pot_rom(num_pot_basis, num_pot_basis, false); CAROM::Matrix pot_rom_inv(num_pot_basis, num_pot_basis, false); CAROM::Matrix pot_rhs_rom(num_pot_basis, num_pot_basis, false); CAROM::Vector pot_rhs_rescaler(num_pot_basis, false); - + /* Load ROM operator */ // read the file from PE0 only if (MPIdata::onpe0) @@ -286,19 +303,19 @@ void runPoissonROM(MGmolInterface *mgmol_) CAROM_VERIFY(num_pot_basis_file == num_pot_basis); h5_helper.getDoubleArray("potential_rom_operator", pot_rom.getData(), - num_pot_basis * num_pot_basis, false); + num_pot_basis * num_pot_basis, false); /* load the inverse as well */ h5_helper.getDoubleArray("potential_rom_inverse", pot_rom_inv.getData(), - num_pot_basis * num_pot_basis, false); + num_pot_basis * num_pot_basis, false); /* load right-hand side hyper-reduction operator */ - h5_helper.getDoubleArray("potential_rhs_rom_inverse", pot_rhs_rom.getData(), - num_pot_basis * num_pot_basis, false); + h5_helper.getDoubleArray("potential_rhs_rom_inverse", + pot_rhs_rom.getData(), num_pot_basis * num_pot_basis, false); /* load right-hand side rescaling operator */ - h5_helper.getDoubleArray("potential_rhs_rescaler", pot_rhs_rescaler.getData(), - num_pot_basis, false); + h5_helper.getDoubleArray("potential_rhs_rescaler", + pot_rhs_rescaler.getData(), num_pot_basis, false); h5_helper.close(); } @@ -307,7 +324,7 @@ void runPoissonROM(MGmolInterface *mgmol_) /* test routines */ template -void testROMPoissonOperator(MGmolInterface *mgmol_) +void testROMPoissonOperator(MGmolInterface* mgmol_) { Control& ct = *(Control::instance()); Mesh* mymesh = Mesh::instance(); @@ -316,14 +333,14 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) ROMPrivateOptions rom_options = ct.getROMOptions(); /* Load MGmol pointer and Potentials */ - MGmol *mgmol = static_cast *>(mgmol_); - Poisson *poisson = mgmol->electrostat_->getPoissonSolver(); - Potentials& pot = mgmol->getHamiltonian()->potential(); - const int dim = pot.size(); + MGmol* mgmol = static_cast*>(mgmol_); + Poisson* poisson = mgmol->electrostat_->getPoissonSolver(); + Potentials& pot = mgmol->getHamiltonian()->potential(); + const int dim = pot.size(); printf("pot size: %d\n", dim); /* GridFunc initialization inputs */ - const pb::Grid &grid(poisson->vh().grid()); + const pb::Grid& grid(poisson->vh().grid()); short bc[3]; for (int d = 0; d < 3; d++) bc[d] = poisson->vh().bc(d); @@ -342,7 +359,7 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) rhs[s].resize(dim); for (int d = 0; d < dim; d++) rhs[s][d] = ran0(); - + /* average out for periodic bc */ pb::GridFunc rhs_gf(grid, bc[0], bc[1], bc[2]); rhs_gf.assign(rhs[s].data(), 'd'); @@ -364,7 +381,7 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) /* apply Laplace operator */ poisson->applyOperator(sol_gf, res); /* FD operator scales rhs by 4pi */ - res.axpy(- 4. * M_PI, rhs_gf); + res.axpy(-4. * M_PI, rhs_gf); printf("FOM res norm: %.3e\n", res.norm2()); } @@ -379,20 +396,21 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) basis_generator.endSamples(); /* Load POD basis. We use the maximum number of basis vectors. */ - const CAROM::Matrix *pot_basis = basis_generator.getSpatialBasis(); + const CAROM::Matrix* pot_basis = basis_generator.getSpatialBasis(); /* Check if full projection preserves FOM solution */ for (int c = 0; c < nsnapshot; c++) { - CAROM::Vector *fom_sol_vec = nullptr; + CAROM::Vector* fom_sol_vec = nullptr; /* get librom view-vector of fom_sol */ - fom_sol_vec = new CAROM::Vector(fom_sol[c].data(), pot_basis->numRows(), true, false); + fom_sol_vec = new CAROM::Vector( + fom_sol[c].data(), pot_basis->numRows(), true, false); - CAROM::Vector *rom_proj = pot_basis->transposeMult(*fom_sol_vec); - CAROM::Vector *reconstruct = pot_basis->mult(*rom_proj); + CAROM::Vector* rom_proj = pot_basis->transposeMult(*fom_sol_vec); + CAROM::Vector* reconstruct = pot_basis->mult(*rom_proj); /* error on libROM side */ - CAROM::Vector *librom_error = reconstruct->minus(fom_sol_vec); + CAROM::Vector* librom_error = reconstruct->minus(fom_sol_vec); printf("librom reconstruction error: %.3e\n", librom_error->norm()); /* error on mgmol side */ @@ -440,12 +458,12 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) printf("ROM res norm: %.3e\n", rom_res.norm()); /* FOM residual: FD operator scales rhs by 4pi */ - res.axpy(- 4. * M_PI, mgmol_rhs); + res.axpy(-4. * M_PI, mgmol_rhs); printf("FOM res norm: %.3e\n", res.norm2()); /* projection of the residual */ res.init_vect(fom_res.getData(), 'd'); - CAROM::Vector *res_proj = pot_basis->transposeMult(fom_res); + CAROM::Vector* res_proj = pot_basis->transposeMult(fom_res); printf("FOM res projection norm: %.3e\n", res_proj->norm()); delete res_proj; @@ -462,7 +480,7 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) CAROM::Vector op_col(pot_basis->numRows(), true); /* copy c-th column librom vector to GridFunc gf_col */ - CAROM::Vector *col = pot_basis->getColumn(c); + CAROM::Vector* col = pot_basis->getColumn(c); col_gf.assign(col->getData(), 'd'); /* apply Laplace operator */ @@ -473,7 +491,7 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) /* Compute basis projection of the column */ /* Resulting vector is undistributed */ - CAROM::Vector *rom_col = pot_basis->transposeMult(op_col); + CAROM::Vector* rom_col = pot_basis->transposeMult(op_col); /* libROM matrix is row-major, so data copy is necessary */ for (int r = 0; r < nsnapshot; r++) @@ -481,14 +499,14 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) delete col; delete rom_col; - } // for (int c = 0; c < num_pot_basis; c++) + } // for (int c = 0; c < num_pot_basis; c++) /* Inverse of the projection ROM matrix */ CAROM::Matrix pot_rom_inv(pot_rom); pot_rom_inv.inverse(); /* Check the inverse */ - CAROM::Matrix *identity = pot_rom_inv.mult(pot_rom); + CAROM::Matrix* identity = pot_rom_inv.mult(pot_rom); printf("pot_rom_inv * pot_rom = identity\n"); for (int i = 0; i < nsnapshot; i++) { @@ -498,8 +516,9 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) } delete identity; - /* Test with sample RHS. ROM must be able to 100% reproduce the FOM solution. */ - std::vector rom_sol(0), rom_rhs(0); + /* Test with sample RHS. ROM must be able to 100% reproduce the FOM + * solution. */ + std::vector rom_sol(0), rom_rhs(0); std::vector> test_sol(nsnapshot); for (int s = 0; s < nsnapshot; s++) { @@ -516,7 +535,7 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) rom_sol.push_back(pot_rom_inv.mult(*rom_rhs.back())); /* check ROM solution */ - CAROM::Vector &res(*pot_rom.mult(*rom_sol.back())); + CAROM::Vector& res(*pot_rom.mult(*rom_sol.back())); res -= *rom_rhs.back(); printf("rom res norm: %.3e\n", res.norm()); @@ -540,8 +559,7 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) double rel_error = testsol_gf.norm2() / fomsol_gf.norm2(); printf("%d-th sample relative error: %.3e\n", s, rel_error); - if (rel_error > 1.0e-9) - abort(); + if (rel_error > 1.0e-9) abort(); } /* clean up pointers */ @@ -553,15 +571,15 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) } template -void testROMRhoOperator(MGmolInterface *mgmol_) +void testROMRhoOperator(MGmolInterface* mgmol_) { Control& ct = *(Control::instance()); Mesh* mymesh = Mesh::instance(); - const int subdivx = mymesh->subdivx(); + const int subdivx = mymesh->subdivx(); const pb::PEenv& myPEenv = mymesh->peenv(); - MGmol_MPI& mmpi = *(MGmol_MPI::instance()); - const int rank = mmpi.mypeGlobal(); - const int nprocs = mmpi.size(); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + const int rank = mmpi.mypeGlobal(); + const int nprocs = mmpi.size(); // if (ct.isLocMode()) // printf("LocMode is On!\n"); @@ -571,11 +589,11 @@ void testROMRhoOperator(MGmolInterface *mgmol_) ROMPrivateOptions rom_options = ct.getROMOptions(); /* Load MGmol pointer and Potentials */ - MGmol *mgmol = static_cast *>(mgmol_); - Poisson *poisson = mgmol->electrostat_->getPoissonSolver(); - Potentials& pot = mgmol->getHamiltonian()->potential(); + MGmol* mgmol = static_cast*>(mgmol_); + Poisson* poisson = mgmol->electrostat_->getPoissonSolver(); + Potentials& pot = mgmol->getHamiltonian()->potential(); std::shared_ptr> rho = NULL; // mgmol->getRho(); - const OrthoType ortho_type = rho->getOrthoType(); + const OrthoType ortho_type = rho->getOrthoType(); assert(ortho_type == OrthoType::Nonorthogonal); /* potential should have the same size as rho */ @@ -584,13 +602,13 @@ void testROMRhoOperator(MGmolInterface *mgmol_) /* number of restart files, start/end indices */ assert(rom_options.restart_file_minidx >= 0); assert(rom_options.restart_file_maxidx >= 0); - const int minidx = rom_options.restart_file_minidx; - const int maxidx = rom_options.restart_file_maxidx; + const int minidx = rom_options.restart_file_minidx; + const int maxidx = rom_options.restart_file_maxidx; const int num_restart = maxidx - minidx + 1; - + /* Initialize libROM classes */ /* - In practice, we do not produce rho POD basis. + In practice, we do not produce rho POD basis. This rho POD basis is for the sake of verification. */ CAROM::Options svd_options(dim, num_restart, 1); @@ -604,19 +622,20 @@ void testROMRhoOperator(MGmolInterface *mgmol_) filename = string_format(rom_options.restart_file_fmt, k); mgmol->loadRestartFile(filename); basis_generator.takeSample(&rho->rho_[0][0]); - } + } // basis_generator.writeSnapshot(); const CAROM::Matrix rho_snapshots(*basis_generator.getSnapshotMatrix()); basis_generator.endSamples(); - const CAROM::Matrix *rho_basis = basis_generator.getSpatialBasis(); - CAROM::Matrix *proj_rho = rho_basis->transposeMult(rho_snapshots); + const CAROM::Matrix* rho_basis = basis_generator.getSpatialBasis(); + CAROM::Matrix* proj_rho = rho_basis->transposeMult(rho_snapshots); /* DEIM hyperreduction */ CAROM::Matrix rho_basis_inv(num_restart, num_restart, false); - std::vector global_sampled_row(num_restart), sampled_rows_per_proc(nprocs); + std::vector global_sampled_row(num_restart), + sampled_rows_per_proc(nprocs); DEIM(rho_basis, num_restart, global_sampled_row, sampled_rows_per_proc, - rho_basis_inv, rank, nprocs); + rho_basis_inv, rank, nprocs); if (rank == 0) { int num_sample_rows = 0; @@ -624,11 +643,12 @@ void testROMRhoOperator(MGmolInterface *mgmol_) num_sample_rows += sampled_rows_per_proc[k]; printf("number of sampled row: %d\n", num_sample_rows); } - + /* get local sampled row */ std::vector offsets, sampled_row(sampled_rows_per_proc[rank]); - int num_global_sample = CAROM::get_global_offsets(sampled_rows_per_proc[rank], offsets); - for (int s = 0, gs = offsets[rank]; gs < offsets[rank+1]; gs++, s++) + int num_global_sample + = CAROM::get_global_offsets(sampled_rows_per_proc[rank], offsets); + for (int s = 0, gs = offsets[rank]; gs < offsets[rank + 1]; gs++, s++) sampled_row[s] = global_sampled_row[gs]; /* load only the first restart file for now */ @@ -639,22 +659,25 @@ void testROMRhoOperator(MGmolInterface *mgmol_) currently, this does not update rho. computeRhoOnSamplePts computes with the new density matrix, while mgmol rho remains the same as the initial condition. - Commenting line 516 gives a consistent result, both for the initial condition. + Commenting line 516 gives a consistent result, both for the initial + condition. */ mgmol->loadRestartFile(filename); const int nrows = mymesh->locNumpt(); // printf("mesh::locNumpt: %d\n", nrows); - - OrbitalsType *orbitals = mgmol->getOrbitals(); + + OrbitalsType* orbitals = mgmol->getOrbitals(); // printf("orbitals::locNumpt: %d\n", orbitals->getLocNumpt()); /* NOTE(kevin): we assume we only use ProjectedMatrices class */ // ProjectedMatrices> *proj_matrices = - // static_cast> *>(orbitals->getProjMatrices()); - ProjectedMatricesInterface *proj_matrices = orbitals->getProjMatrices(); + // static_cast> + // *>(orbitals->getProjMatrices()); + ProjectedMatricesInterface* proj_matrices = orbitals->getProjMatrices(); proj_matrices->updateSubMatX(); - SquareLocalMatrices& localX(proj_matrices->getLocalX()); + SquareLocalMatrices& localX( + proj_matrices->getLocalX()); // printf("localX nmat: %d\n", localX.nmat()); // printf("localX n: %d\n", localX.n()); @@ -664,14 +687,14 @@ void testROMRhoOperator(MGmolInterface *mgmol_) assert(!dm_distributed); /* copy density matrix */ - CAROM::Matrix dm(localX.getRawPtr(), localX.m(), localX.n(), dm_distributed, true); + CAROM::Matrix dm( + localX.getRawPtr(), localX.m(), localX.n(), dm_distributed, true); // /* random density matrix */ - // /* NOTE(kevin): Due to rescaleTotalCharge, the result is slightly inconsistent. */ - // CAROM::Matrix dm(localX.m(), localX.n(), dm_distributed, true); - // dm = 0.0; - // dist_matrix::DistMatrix dm_mgmol(proj_matrices->dm()); - // for (int i = 0; i < dm_mgmol.m(); i++) + // /* NOTE(kevin): Due to rescaleTotalCharge, the result is slightly + // inconsistent. */ CAROM::Matrix dm(localX.m(), localX.n(), dm_distributed, + // true); dm = 0.0; dist_matrix::DistMatrix + // dm_mgmol(proj_matrices->dm()); for (int i = 0; i < dm_mgmol.m(); i++) // { // for (int j = 0; j < dm_mgmol.m(); j++) // { @@ -691,8 +714,8 @@ void testROMRhoOperator(MGmolInterface *mgmol_) CAROM::Matrix psi(dim, chrom_num, true); for (int c = 0; c < chrom_num; c++) { - ORBDTYPE *d_psi = orbitals->getPsi(c); - for (int d = 0; d < dim ; d++) + ORBDTYPE* d_psi = orbitals->getPsi(c); + for (int d = 0; d < dim; d++) psi.item(d, c) = *(d_psi + d); } @@ -710,22 +733,27 @@ void testROMRhoOperator(MGmolInterface *mgmol_) { const double error = abs(rho->rho_[0][sampled_row[s]] - sample_rho(s)); if (error > 1.0e-4) - printf("rank %d, rho[%d]: %.5e, sample_rho: %.5e, librom_snapshot: %.5e\n", - rank, sampled_row[s], rho->rho_[0][sampled_row[s]], sample_rho(s), rho_snapshots(sampled_row[s], test_idx)); + printf("rank %d, rho[%d]: %.5e, sample_rho: %.5e, librom_snapshot: " + "%.5e\n", + rank, sampled_row[s], rho->rho_[0][sampled_row[s]], + sample_rho(s), rho_snapshots(sampled_row[s], test_idx)); CAROM_VERIFY(error < 1.0e-4); } sample_rho.gather(); - CAROM::Vector *rom_rho = rho_basis_inv.mult(sample_rho); + CAROM::Vector* rom_rho = rho_basis_inv.mult(sample_rho); for (int d = 0; d < rom_rho->dim(); d++) { - if ((rank == 0) && (abs(proj_rho->item(d, test_idx) - rom_rho->item(d)) > 1.0e-3)) - printf("rom_rho error: %.3e\n", abs(proj_rho->item(d, test_idx) - rom_rho->item(d))); - CAROM_VERIFY(abs(proj_rho->item(d, test_idx) - rom_rho->item(d)) < 1.0e-3); + if ((rank == 0) + && (abs(proj_rho->item(d, test_idx) - rom_rho->item(d)) > 1.0e-3)) + printf("rom_rho error: %.3e\n", + abs(proj_rho->item(d, test_idx) - rom_rho->item(d))); + CAROM_VERIFY( + abs(proj_rho->item(d, test_idx) - rom_rho->item(d)) < 1.0e-3); } - CAROM::Vector *fom_rho = rho_basis->mult(*rom_rho); + CAROM::Vector* fom_rho = rho_basis->mult(*rom_rho); CAROM_VERIFY(fom_rho->dim() == rho->rho_[0].size()); for (int d = 0; d < fom_rho->dim(); d++) @@ -738,13 +766,13 @@ void testROMRhoOperator(MGmolInterface *mgmol_) /* dm: density matrix converted to CAROM::Matrix phi_basis: POD basis matrix of orbitals, or orbitals themselves - rom_phi: ROM coefficients of POD Basis. If phi_basis is orbitals themselves, then simply an identity. - local_idx: sampled local grid indices on this processor. - sampled_rho: FOM density on sampled grid points. + rom_phi: ROM coefficients of POD Basis. If phi_basis is orbitals themselves, + then simply an identity. local_idx: sampled local grid indices on this + processor. sampled_rho: FOM density on sampled grid points. */ -void computeRhoOnSamplePts(const CAROM::Matrix &dm, - const CAROM::Matrix &phi_basis, const CAROM::Matrix &rom_phi, - const std::vector &local_idx, CAROM::Vector &sampled_rho) +void computeRhoOnSamplePts(const CAROM::Matrix& dm, + const CAROM::Matrix& phi_basis, const CAROM::Matrix& rom_phi, + const std::vector& local_idx, CAROM::Vector& sampled_rho) { assert(!dm.distributed()); assert(sampled_rho.distributed() == phi_basis.distributed()); @@ -759,8 +787,8 @@ void computeRhoOnSamplePts(const CAROM::Matrix &dm, sampled_phi.mult(dm, product); sampled_rho.setSize(sampled_phi.numRows()); - double *d_product = product.getData(); - double *d_phi = sampled_phi.getData(); + double* d_product = product.getData(); + double* d_phi = sampled_phi.getData(); for (int d = 0; d < sampled_rho.dim(); d++) { double val = 0.0; @@ -773,7 +801,8 @@ void computeRhoOnSamplePts(const CAROM::Matrix &dm, /* TODO(kevin): need to figure out what these functions do. - and probably need to make ROM-equivalent functions with another hyper-reduction? + and probably need to make ROM-equivalent functions with another + hyper-reduction? */ // gatherSpin(); @@ -784,38 +813,42 @@ void computeRhoOnSamplePts(const CAROM::Matrix &dm, } template -void testROMIonDensity(MGmolInterface *mgmol_) +void testROMIonDensity(MGmolInterface* mgmol_) { /* random number generator */ - static std::random_device rd; // Will be used to obtain a seed for the random number engine - static std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd(){} + static std::random_device + rd; // Will be used to obtain a seed for the random number engine + static std::mt19937 gen( + rd()); // Standard mersenne_twister_engine seeded with rd(){} static std::uniform_real_distribution<> dis(0.0, 1.0); Control& ct = *(Control::instance()); Mesh* mymesh = Mesh::instance(); - const pb::Grid& mygrid = mymesh->grid(); - const int subdivx = mymesh->subdivx(); + const pb::Grid& mygrid = mymesh->grid(); + const int subdivx = mymesh->subdivx(); const pb::PEenv& myPEenv = mymesh->peenv(); - MGmol_MPI& mmpi = *(MGmol_MPI::instance()); - const int rank = mmpi.mypeGlobal(); - const int nprocs = mmpi.size(); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + const int rank = mmpi.mypeGlobal(); + const int nprocs = mmpi.size(); ROMPrivateOptions rom_options = ct.getROMOptions(); /* Load MGmol pointer and Potentials */ - MGmol *mgmol = static_cast *>(mgmol_); - Poisson *poisson = mgmol->electrostat_->getPoissonSolver(); - Potentials& pot = mgmol->getHamiltonian()->potential(); - const int dim = pot.size(); - std::shared_ptr ions = mgmol->getIons(); + MGmol* mgmol = static_cast*>(mgmol_); + Poisson* poisson = mgmol->electrostat_->getPoissonSolver(); + Potentials& pot = mgmol->getHamiltonian()->potential(); + const int dim = pot.size(); + std::shared_ptr ions = mgmol->getIons(); /* get the extent of global domain */ - const double origin[3] = { mygrid.origin(0), mygrid.origin(1), mygrid.origin(2) }; + const double origin[3] + = { mygrid.origin(0), mygrid.origin(1), mygrid.origin(2) }; const double lattice[3] = { mygrid.ll(0), mygrid.ll(1), mygrid.ll(2) }; if (rank == 0) { printf("origin: (%.3e, %.3e, %.3e)\n", origin[0], origin[1], origin[2]); - printf("lattice: (%.3e, %.3e, %.3e)\n", lattice[0], lattice[1], lattice[2]); + printf("lattice: (%.3e, %.3e, %.3e)\n", lattice[0], lattice[1], + lattice[2]); } /* get global atomic numbers */ @@ -837,7 +870,7 @@ void testROMIonDensity(MGmolInterface *mgmol_) for (int k = 0; k < num_ions; k++) for (int d = 0; d < 3; d++) cfgs[idx][3 * k + d] = origin[d] + lattice[d] * dis(gen); - + mmpi.bcastGlobal(cfgs[idx].data(), 3 * num_ions, 0); } @@ -858,7 +891,8 @@ void testROMIonDensity(MGmolInterface *mgmol_) { fom_overlap_ions[idx][k].resize(3); for (int d = 0; d < 3; d++) - fom_overlap_ions[idx][k][d] = new_ions->overlappingVL_ions()[k]->position(d); + fom_overlap_ions[idx][k][d] + = new_ions->overlappingVL_ions()[k]->position(d); } /* compute resulting ion density */ @@ -873,7 +907,7 @@ void testROMIonDensity(MGmolInterface *mgmol_) /* Initialize libROM classes */ /* - In practice, we do not produce rhoc POD basis. + In practice, we do not produce rhoc POD basis. This rhoc POD basis is for the sake of verification. */ CAROM::Options svd_options(dim, num_snap, 1); @@ -885,14 +919,15 @@ void testROMIonDensity(MGmolInterface *mgmol_) const CAROM::Matrix rhoc_snapshots(*basis_generator.getSnapshotMatrix()); basis_generator.endSamples(); - const CAROM::Matrix *rhoc_basis = basis_generator.getSpatialBasis(); - CAROM::Matrix *proj_rhoc = rhoc_basis->transposeMult(rhoc_snapshots); + const CAROM::Matrix* rhoc_basis = basis_generator.getSpatialBasis(); + CAROM::Matrix* proj_rhoc = rhoc_basis->transposeMult(rhoc_snapshots); /* DEIM hyperreduction */ CAROM::Matrix rhoc_basis_inv(num_snap, num_snap, false); - std::vector global_sampled_row(num_snap), sampled_rows_per_proc(nprocs); + std::vector global_sampled_row(num_snap), + sampled_rows_per_proc(nprocs); DEIM(rhoc_basis, num_snap, global_sampled_row, sampled_rows_per_proc, - rhoc_basis_inv, rank, nprocs); + rhoc_basis_inv, rank, nprocs); if (rank == 0) { int num_sample_rows = 0; @@ -900,15 +935,16 @@ void testROMIonDensity(MGmolInterface *mgmol_) num_sample_rows += sampled_rows_per_proc[k]; printf("number of sampled row: %d\n", num_sample_rows); } - + /* get local sampled row */ std::vector offsets, sampled_row(sampled_rows_per_proc[rank]); - int num_global_sample = CAROM::get_global_offsets(sampled_rows_per_proc[rank], offsets); - for (int s = 0, gs = offsets[rank]; gs < offsets[rank+1]; gs++, s++) + int num_global_sample + = CAROM::get_global_offsets(sampled_rows_per_proc[rank], offsets); + for (int s = 0, gs = offsets[rank]; gs < offsets[rank + 1]; gs++, s++) sampled_row[s] = global_sampled_row[gs]; /* test one solution */ - std::uniform_int_distribution<> distrib(0, num_snap-1); + std::uniform_int_distribution<> distrib(0, num_snap - 1); int test_idx = distrib(gen); mmpi.bcastGlobal(&test_idx); if (rank == 0) printf("test index: %d\n", test_idx); @@ -917,10 +953,13 @@ void testROMIonDensity(MGmolInterface *mgmol_) new_ions = new Ions(cfgs[test_idx], atnumbers, lattice, ions->getSpecies()); /* Sanity check for overlapping ions */ - CAROM_VERIFY(fom_overlap_ions[test_idx].size() == new_ions->overlappingVL_ions().size()); + CAROM_VERIFY(fom_overlap_ions[test_idx].size() + == new_ions->overlappingVL_ions().size()); for (int k = 0; k < new_ions->overlappingVL_ions().size(); k++) for (int d = 0; d < 3; d++) - CAROM_VERIFY(abs(fom_overlap_ions[test_idx][k][d] - new_ions->overlappingVL_ions()[k]->position(d)) < 1.0e-12); + CAROM_VERIFY(abs(fom_overlap_ions[test_idx][k][d] + - new_ions->overlappingVL_ions()[k]->position(d)) + < 1.0e-12); /* set up potentials */ mgmol->setupPotentials(*new_ions); @@ -929,33 +968,48 @@ void testROMIonDensity(MGmolInterface *mgmol_) std::vector sampled_rhoc(sampled_row.size()); pot.evalIonDensityOnSamplePts(*new_ions, sampled_row, sampled_rhoc); - // For now, we relax the threshold to allow the slight difference of the values at the sampled indices. - // This is because the rescaling procedure after getting the radial data on mesh requires global information, - // thus cannot be done in the ROM level only with sampled values. - // After we have the DEIM reconstruction, we will do rescaling and revisit this. + // For now, we relax the threshold to allow the slight difference of the + // values at the sampled indices. This is because the rescaling procedure + // after getting the radial data on mesh requires global information, thus + // cannot be done in the ROM level only with sampled values. After we have + // the DEIM reconstruction, we will do rescaling and revisit this. for (int d = 0; d < sampled_row.size(); d++) { - printf("rank %d, fom rhoc[%d]: %.3e, rom rhoc: %.3e\n", rank, sampled_row[d], fom_rhoc[test_idx][sampled_row[d]], sampled_rhoc[d]); - CAROM_VERIFY(abs(fom_rhoc[test_idx][sampled_row[d]] - sampled_rhoc[d]) < 1.0e-2); + printf("rank %d, fom rhoc[%d]: %.3e, rom rhoc: %.3e\n", rank, + sampled_row[d], fom_rhoc[test_idx][sampled_row[d]], + sampled_rhoc[d]); + CAROM_VERIFY( + abs(fom_rhoc[test_idx][sampled_row[d]] - sampled_rhoc[d]) < 1.0e-2); } delete new_ions; } -template void readRestartFiles>(MGmolInterface *mgmol_); -template void readRestartFiles>(MGmolInterface *mgmol_); - -template void buildROMPoissonOperator>(MGmolInterface *mgmol_); -template void buildROMPoissonOperator>(MGmolInterface *mgmol_); - -template void runPoissonROM>(MGmolInterface *mgmol_); -template void runPoissonROM>(MGmolInterface *mgmol_); - -template void testROMPoissonOperator>(MGmolInterface *mgmol_); -template void testROMPoissonOperator>(MGmolInterface *mgmol_); - -template void testROMRhoOperator>(MGmolInterface *mgmol_); -template void testROMRhoOperator>(MGmolInterface *mgmol_); - -template void testROMIonDensity>(MGmolInterface *mgmol_); -template void testROMIonDensity>(MGmolInterface *mgmol_); +template void readRestartFiles>( + MGmolInterface* mgmol_); +template void readRestartFiles>( + MGmolInterface* mgmol_); + +template void buildROMPoissonOperator>( + MGmolInterface* mgmol_); +template void buildROMPoissonOperator>( + MGmolInterface* mgmol_); + +template void runPoissonROM>(MGmolInterface* mgmol_); +template void runPoissonROM>( + MGmolInterface* mgmol_); + +template void testROMPoissonOperator>( + MGmolInterface* mgmol_); +template void testROMPoissonOperator>( + MGmolInterface* mgmol_); + +template void testROMRhoOperator>( + MGmolInterface* mgmol_); +template void testROMRhoOperator>( + MGmolInterface* mgmol_); + +template void testROMIonDensity>( + MGmolInterface* mgmol_); +template void testROMIonDensity>( + MGmolInterface* mgmol_); diff --git a/src/rom_workflows.h b/src/rom_workflows.h index 9222249f..3dcaf9ba 100644 --- a/src/rom_workflows.h +++ b/src/rom_workflows.h @@ -12,13 +12,13 @@ #include "Control.h" #include "ExtendedGridOrbitals.h" -#include "ProjectedMatrices.h" #include "LocGridOrbitals.h" -#include "Potentials.h" #include "MGmol.h" #include "MGmol_MPI.h" #include "MPIdata.h" #include "Mesh.h" +#include "Potentials.h" +#include "ProjectedMatrices.h" #include "mgmol_run.h" #include "tools.h" @@ -40,25 +40,25 @@ namespace po = boost::program_options; #include "utils/mpi_utils.h" template -void readRestartFiles(MGmolInterface *mgmol_); +void readRestartFiles(MGmolInterface* mgmol_); template -void buildROMPoissonOperator(MGmolInterface *mgmol_); +void buildROMPoissonOperator(MGmolInterface* mgmol_); template -void runPoissonROM(MGmolInterface *mgmol_); +void runPoissonROM(MGmolInterface* mgmol_); template -void testROMPoissonOperator(MGmolInterface *mgmol_); +void testROMPoissonOperator(MGmolInterface* mgmol_); template -void testROMRhoOperator(MGmolInterface *mgmol_); +void testROMRhoOperator(MGmolInterface* mgmol_); template -void testROMIonDensity(MGmolInterface *mgmol_); +void testROMIonDensity(MGmolInterface* mgmol_); -void computeRhoOnSamplePts(const CAROM::Matrix &dm, - const CAROM::Matrix &phi_basis, const CAROM::Matrix &rom_phi, - const std::vector &local_idx, CAROM::Vector &sampled_rho); +void computeRhoOnSamplePts(const CAROM::Matrix& dm, + const CAROM::Matrix& phi_basis, const CAROM::Matrix& rom_phi, + const std::vector& local_idx, CAROM::Vector& sampled_rho); -#endif // ROM_WORKFLOWS_H +#endif // ROM_WORKFLOWS_H diff --git a/src/sparse_linear_algebra/DataDistribution.cc b/src/sparse_linear_algebra/DataDistribution.cc index 67bb1319..ae613b71 100644 --- a/src/sparse_linear_algebra/DataDistribution.cc +++ b/src/sparse_linear_algebra/DataDistribution.cc @@ -277,7 +277,7 @@ void DataDistribution::distributeLocalDataWithCommOvlp(const int nsteps, if (nsteps > 0) { #ifndef NDEBUG // check if receive buffers are large enough - //#if 1 + // #if 1 int remote_size; MPI_Irecv(&remote_size, 1, MPI_INT, source, 0, cart_comm_, &request[0]); MPI_Isend(&siz, 1, MPI_INT, dest, 0, cart_comm_, &request[1]); @@ -401,7 +401,7 @@ void DataDistribution::distributeLocalData(const int nsteps, const int dir, if ((step == nsteps - 1) && (bcflag - == true)) // this is the last step - apply boundary condition + == true)) // this is the last step - apply boundary condition { packed_buffer.updateMatrixEntriesWithRecvBuf( amat, packed_buffer.recvBuffer()); diff --git a/src/sparse_linear_algebra/PreconILU.cc b/src/sparse_linear_algebra/PreconILU.cc index bf46aec5..a5158c02 100644 --- a/src/sparse_linear_algebra/PreconILU.cc +++ b/src/sparse_linear_algebra/PreconILU.cc @@ -353,7 +353,7 @@ int PreconILU::lofC(LinearSolverMatrix& csmat_) { int pos = row - (*L_).getColumnIterator( - 0); /* get the position of row */ + 0); /* get the position of row */ int it = llvl[pos] + *uptr + 1; if (it > lof_) continue; int ip = iw[*row]; @@ -887,7 +887,7 @@ int PreconILU::milut(LinearSolverMatrix& csmat_) { alpha = sr * (std::sqrt(tval) - + droptol_ * (gamma - std::sqrt(tval))); + + droptol_ * (gamma - std::sqrt(tval))); } } else diff --git a/src/sparse_linear_algebra/SparseRow.h b/src/sparse_linear_algebra/SparseRow.h index 325841d0..af36bd35 100644 --- a/src/sparse_linear_algebra/SparseRow.h +++ b/src/sparse_linear_algebra/SparseRow.h @@ -29,7 +29,7 @@ typedef enum INSERTMODE } INSERTMODE; /* define maximum local matrix size */ -//#define MAX_MAT_SIZE 15000 +// #define MAX_MAT_SIZE 15000 /* define default tolerance for pruning matrix entries */ #define MAT_TOL 1.0e-14 /* define maximum number of print rows */ diff --git a/src/sparse_linear_algebra/Table.h b/src/sparse_linear_algebra/Table.h index 3ed9e39d..14a668c5 100644 --- a/src/sparse_linear_algebra/Table.h +++ b/src/sparse_linear_algebra/Table.h @@ -21,7 +21,7 @@ Adapted from pARMS code. (original version by Z. Li) #include #define USE_POWERS2 1 -//#define COUNT_LINKS 1 +// #define COUNT_LINKS 1 typedef struct Slot { diff --git a/src/sparse_linear_algebra/VariableSizeMatrixInterface.h b/src/sparse_linear_algebra/VariableSizeMatrixInterface.h index 98b9f8ab..275c27df 100644 --- a/src/sparse_linear_algebra/VariableSizeMatrixInterface.h +++ b/src/sparse_linear_algebra/VariableSizeMatrixInterface.h @@ -38,7 +38,7 @@ class VariableSizeMatrixInterface static Timer AmultSymB_tm_; public: - virtual ~VariableSizeMatrixInterface() {} + virtual ~VariableSizeMatrixInterface() { } static void printTimers(std::ostream& os) { diff --git a/src/tools/PinnedH2O.cc b/src/tools/PinnedH2O.cc index 22f4677a..540babb8 100644 --- a/src/tools/PinnedH2O.cc +++ b/src/tools/PinnedH2O.cc @@ -17,32 +17,42 @@ PinnedH2O::PinnedH2O() H2_idx(-1), planar_rotation_angle(0.0) { - for (int i = 0; i < 3; ++i) { - for (int j = 0; j < 3; ++j) { + for (int i = 0; i < 3; ++i) + { + for (int j = 0; j < 3; ++j) + { out_of_plane_rotation_matrix[i][j] = 0.0; } } } -double PinnedH2O::calculate_bondlength(const double atom1[3], const double atom2[3]) const +double PinnedH2O::calculate_bondlength( + const double atom1[3], const double atom2[3]) const { - return sqrt(pow(atom1[0] - atom2[0], 2) + pow(atom1[1] - atom2[1], 2) + pow(atom1[2] - atom2[2], 2)); + return sqrt(pow(atom1[0] - atom2[0], 2) + pow(atom1[1] - atom2[1], 2) + + pow(atom1[2] - atom2[2], 2)); } -double PinnedH2O::calculate_bondangle(const double atom1[3], const double atom2[3], const double atom3[3]) const +double PinnedH2O::calculate_bondangle( + const double atom1[3], const double atom2[3], const double atom3[3]) const { - double vector1[3] = {atom1[0] - atom2[0], atom1[1] - atom2[1], atom1[2] - atom2[2]}; - double vector2[3] = {atom3[0] - atom2[0], atom3[1] - atom2[1], atom3[2] - atom2[2]}; - - double dot_product = vector1[0] * vector2[0] + vector1[1] * vector2[1] + vector1[2] * vector2[2]; - double magnitude_product = sqrt(pow(vector1[0], 2) + pow(vector1[1], 2) + pow(vector1[2], 2)) * - sqrt(pow(vector2[0], 2) + pow(vector2[1], 2) + pow(vector2[2], 2)); + double vector1[3] + = { atom1[0] - atom2[0], atom1[1] - atom2[1], atom1[2] - atom2[2] }; + double vector2[3] + = { atom3[0] - atom2[0], atom3[1] - atom2[1], atom3[2] - atom2[2] }; + + double dot_product = vector1[0] * vector2[0] + vector1[1] * vector2[1] + + vector1[2] * vector2[2]; + double magnitude_product + = sqrt(pow(vector1[0], 2) + pow(vector1[1], 2) + pow(vector1[2], 2)) + * sqrt(pow(vector2[0], 2) + pow(vector2[1], 2) + pow(vector2[2], 2)); double angle = acos(dot_product / magnitude_product); return angle; } -void PinnedH2O::rotation_matrix(const double axis[3], double angle, double matrix[3][3]) const +void PinnedH2O::rotation_matrix( + const double axis[3], double angle, double matrix[3][3]) const { double cos_theta = cos(angle); double sin_theta = sin(angle); @@ -69,31 +79,44 @@ void PinnedH2O::normalize(double vec[3]) const vec[2] /= norm; } -void PinnedH2O::cross(const double a[3], const double b[3], double result[3]) const +void PinnedH2O::cross( + const double a[3], const double b[3], double result[3]) const { result[0] = a[1] * b[2] - a[2] * b[1]; result[1] = a[2] * b[0] - a[0] * b[2]; result[2] = a[0] * b[1] - a[1] * b[0]; } -void PinnedH2O::apply_rotation(const double matrix[3][3], const double vec[3], double result[3]) const +void PinnedH2O::apply_rotation( + const double matrix[3][3], const double vec[3], double result[3]) const { - result[0] = matrix[0][0] * vec[0] + matrix[0][1] * vec[1] + matrix[0][2] * vec[2]; - result[1] = matrix[1][0] * vec[0] + matrix[1][1] * vec[1] + matrix[1][2] * vec[2]; - result[2] = matrix[2][0] * vec[0] + matrix[2][1] * vec[1] + matrix[2][2] * vec[2]; + result[0] + = matrix[0][0] * vec[0] + matrix[0][1] * vec[1] + matrix[0][2] * vec[2]; + result[1] + = matrix[1][0] * vec[0] + matrix[1][1] * vec[1] + matrix[1][2] * vec[2]; + result[2] + = matrix[2][0] * vec[0] + matrix[2][1] * vec[1] + matrix[2][2] * vec[2]; } -void PinnedH2O::apply_transpose_rotation(const double matrix[3][3], const double vec[3], double result[3]) const +void PinnedH2O::apply_transpose_rotation( + const double matrix[3][3], const double vec[3], double result[3]) const { - result[0] = matrix[0][0] * vec[0] + matrix[1][0] * vec[1] + matrix[2][0] * vec[2]; - result[1] = matrix[0][1] * vec[0] + matrix[1][1] * vec[1] + matrix[2][1] * vec[2]; - result[2] = matrix[0][2] * vec[0] + matrix[1][2] * vec[1] + matrix[2][2] * vec[2]; + result[0] + = matrix[0][0] * vec[0] + matrix[1][0] * vec[1] + matrix[2][0] * vec[2]; + result[1] + = matrix[0][1] * vec[0] + matrix[1][1] * vec[1] + matrix[2][1] * vec[2]; + result[2] + = matrix[0][2] * vec[0] + matrix[1][2] * vec[1] + matrix[2][2] * vec[2]; } -void PinnedH2O::rotate(std::vector& positions, std::vector& anumbers) +void PinnedH2O::rotate( + std::vector& positions, std::vector& anumbers) { - for (int i = 0; i < 3; i++) { - if (positions[3 * i] == 0.0 && positions[3 * i + 1] == 0.0 && positions[3 * i + 2] == 0.0) { + for (int i = 0; i < 3; i++) + { + if (positions[3 * i] == 0.0 && positions[3 * i + 1] == 0.0 + && positions[3 * i + 2] == 0.0) + { O1_idx = i; break; } @@ -103,13 +126,16 @@ void PinnedH2O::rotate(std::vector& positions, std::vector& anumb H1_idx = (O1_idx + 1) % 3; H2_idx = (O1_idx + 2) % 3; - double O1[3] = {positions[3 * O1_idx], positions[3 * O1_idx + 1], positions[3 * O1_idx + 2]}; - double H1[3] = {positions[3 * H1_idx], positions[3 * H1_idx + 1], positions[3 * H1_idx + 2]}; - double H2[3] = {positions[3 * H2_idx], positions[3 * H2_idx + 1], positions[3 * H2_idx + 2]}; + double O1[3] = { positions[3 * O1_idx], positions[3 * O1_idx + 1], + positions[3 * O1_idx + 2] }; + double H1[3] = { positions[3 * H1_idx], positions[3 * H1_idx + 1], + positions[3 * H1_idx + 2] }; + double H2[3] = { positions[3 * H2_idx], positions[3 * H2_idx + 1], + positions[3 * H2_idx + 2] }; bondlength1 = calculate_bondlength(H1, O1); bondlength2 = calculate_bondlength(H2, O1); - bondangle = calculate_bondangle(H1, O1, H2); + bondangle = calculate_bondangle(H1, O1, H2); double H1_temp[3], H2_temp[3]; double H1_rotated[3], H2_rotated[3]; @@ -117,11 +143,12 @@ void PinnedH2O::rotate(std::vector& positions, std::vector& anumb double plane_normal[3]; cross(H2, H1, plane_normal); normalize(plane_normal); - double target_plane_normal[3] = {0, 0, 1}; - double dot_product = plane_normal[0] * target_plane_normal[0] + - plane_normal[1] * target_plane_normal[1] + - plane_normal[2] * target_plane_normal[2]; - double angle_to_align = std::acos(std::min(std::max(dot_product, -1.0), 1.0)); + double target_plane_normal[3] = { 0, 0, 1 }; + double dot_product = plane_normal[0] * target_plane_normal[0] + + plane_normal[1] * target_plane_normal[1] + + plane_normal[2] * target_plane_normal[2]; + double angle_to_align + = std::acos(std::min(std::max(dot_product, -1.0), 1.0)); double axis_to_align[3]; if (abs(dot_product) > 1.0 - 1e-8) { @@ -135,14 +162,16 @@ void PinnedH2O::rotate(std::vector& positions, std::vector& anumb normalize(axis_to_align); } - rotation_matrix(axis_to_align, angle_to_align, out_of_plane_rotation_matrix); + rotation_matrix( + axis_to_align, angle_to_align, out_of_plane_rotation_matrix); apply_rotation(out_of_plane_rotation_matrix, H1, H1_temp); apply_rotation(out_of_plane_rotation_matrix, H2, H2_temp); - double theta1 = std::atan2(H1_temp[1], H1_temp[0]); + double theta1 = std::atan2(H1_temp[1], H1_temp[0]); planar_rotation_angle = -theta1 + bondangle / 2.0; double planar_rotation_matrix[3][3]; - rotation_matrix(target_plane_normal, planar_rotation_angle, planar_rotation_matrix); + rotation_matrix( + target_plane_normal, planar_rotation_angle, planar_rotation_matrix); apply_rotation(planar_rotation_matrix, H1_temp, H1_rotated); apply_rotation(planar_rotation_matrix, H2_temp, H2_rotated); @@ -170,15 +199,16 @@ void PinnedH2O::rotate(std::vector& positions, std::vector& anumb anumbers[2] = 1; } -void PinnedH2O::transpose_rotate(std::vector& positions, std::vector& anumbers, std::vector& forces) +void PinnedH2O::transpose_rotate(std::vector& positions, + std::vector& anumbers, std::vector& forces) { - double H2_rotated[3] = {positions[0], positions[1], positions[2]}; - double O1_rotated[3] = {positions[3], positions[4], positions[5]}; - double H1_rotated[3] = {positions[6], positions[7], positions[8]}; + double H2_rotated[3] = { positions[0], positions[1], positions[2] }; + double O1_rotated[3] = { positions[3], positions[4], positions[5] }; + double H1_rotated[3] = { positions[6], positions[7], positions[8] }; - double f_H2_rotated[3] = {forces[0], forces[1], forces[2]}; - double f_O1_rotated[3] = {forces[3], forces[4], forces[5]}; - double f_H1_rotated[3] = {forces[6], forces[7], forces[8]}; + double f_H2_rotated[3] = { forces[0], forces[1], forces[2] }; + double f_O1_rotated[3] = { forces[3], forces[4], forces[5] }; + double f_H1_rotated[3] = { forces[6], forces[7], forces[8] }; if (flipped_bond) { @@ -190,8 +220,9 @@ void PinnedH2O::transpose_rotate(std::vector& positions, std::vector& positions, std::vector -#include #include +#include #include +#include class PinnedH2O { @@ -23,21 +23,29 @@ class PinnedH2O ~PinnedH2O() = default; void rotate(std::vector& positions, std::vector& anumbers); - void transpose_rotate(std::vector& positions, std::vector& anumbers, std::vector& forces); + void transpose_rotate(std::vector& positions, + std::vector& anumbers, std::vector& forces); void print(std::ostream& os) { - os << "Bondlengths = " << bondlength1 << ", " << bondlength2 << " Bohrs; " - << "Bondangle = " << bondangle * 180.0 / M_PI << " degrees." << std::endl; + os << "Bondlengths = " << bondlength1 << ", " << bondlength2 + << " Bohrs; " + << "Bondangle = " << bondangle * 180.0 / M_PI << " degrees." + << std::endl; } private: - double calculate_bondlength(const double atom1[3], const double atom2[3]) const; - double calculate_bondangle(const double atom1[3], const double atom2[3], const double atom3[3]) const; - void rotation_matrix(const double axis[3], double angle, double matrix[3][3]) const; + double calculate_bondlength( + const double atom1[3], const double atom2[3]) const; + double calculate_bondangle(const double atom1[3], const double atom2[3], + const double atom3[3]) const; + void rotation_matrix( + const double axis[3], double angle, double matrix[3][3]) const; void normalize(double vec[3]) const; void cross(const double a[3], const double b[3], double result[3]) const; - void apply_rotation(const double matrix[3][3], const double vec[3], double result[3]) const; - void apply_transpose_rotation(const double matrix[3][3], const double vec[3], double result[3]) const; + void apply_rotation( + const double matrix[3][3], const double vec[3], double result[3]) const; + void apply_transpose_rotation( + const double matrix[3][3], const double vec[3], double result[3]) const; double bondlength1; double bondlength2; diff --git a/src/tools/Timer.h b/src/tools/Timer.h index d9cf86a8..02cdd35a 100644 --- a/src/tools/Timer.h +++ b/src/tools/Timer.h @@ -67,7 +67,7 @@ class Timer total_real_(0.0), running_(false), ncalls_(0), - comm_(comm){}; + comm_(comm) {}; void reset() { diff --git a/src/tools/memory_space.h b/src/tools/memory_space.h index 767ecffb..958f05ba 100644 --- a/src/tools/memory_space.h +++ b/src/tools/memory_space.h @@ -227,7 +227,7 @@ struct Memory ptr = nullptr; } - static void free_host_view(T* /*ptr*/) {} + static void free_host_view(T* /*ptr*/) { } static void copy(T const* in, unsigned int size, T* out) { diff --git a/src/tools/mgmol_mpi_tools.cc b/src/tools/mgmol_mpi_tools.cc index 282fd8cb..b926611e 100644 --- a/src/tools/mgmol_mpi_tools.cc +++ b/src/tools/mgmol_mpi_tools.cc @@ -147,7 +147,7 @@ int gatherV(std::vector& sendbuf, } char* recvdata = new char[totchars]; mpi_err = MPI_Gatherv(&charStr[0], sendcount, MPI_CHAR, &recvdata[0], - recvcounts, displs, MPI_CHAR, root, comm); + recvcounts, displs, MPI_CHAR, root, comm); if (mpi_err != MPI_SUCCESS) { std::cerr << "ERROR in MPI_Gatherv in MGmol_MPI::GatherV() !!!" diff --git a/src/unused/KBPsiMatrix.cc b/src/unused/KBPsiMatrix.cc index 6b96a5cb..84e211dd 100644 --- a/src/unused/KBPsiMatrix.cc +++ b/src/unused/KBPsiMatrix.cc @@ -15,7 +15,7 @@ #include -//#define USE_OLD_ALGO 1 +// #define USE_OLD_ALGO 1 Timer KBPsiMatrix::allreduce_tm_("KBPsiMatrix::allreduce"); Timer KBPsiMatrix::global_sum_tm_("KBPsiMatrix::global_sum"); diff --git a/src/unused/radial_functions.cc b/src/unused/radial_functions.cc index 2ff38769..b474f3bd 100644 --- a/src/unused/radial_functions.cc +++ b/src/unused/radial_functions.cc @@ -167,20 +167,20 @@ double rho2(const double r, const double a) { const double ar2 = (a - r) * (a - r); - alpha = (4 * (-8 + exp(3 * r2)) - * (1 - + (-1 + e0) - * (1 - - (1 + pow(a, 2) - 2 * a * r + r2) - / exp(ar2)))) - / exp(4 * r2) - - ((-1 + e0) - * ((2 * (a - r)) / exp(ar2) - + (2 * (1 + ar2 + 0. * pow(a - r, 4)) * (-a + r)) - / exp(ar2)) - * ((2 * (-2 + exp(3 * r2)) * r) / exp(4 * r2) - + sqrt(M_PI) * (erf(2. * r) - erf(r)))) - / r2; + alpha + = (4 * (-8 + exp(3 * r2)) + * (1 + + (-1 + e0) + * (1 + - (1 + pow(a, 2) - 2 * a * r + r2) / exp(ar2)))) + / exp(4 * r2) + - ((-1 + e0) + * ((2 * (a - r)) / exp(ar2) + + (2 * (1 + ar2 + 0. * pow(a - r, 4)) * (-a + r)) + / exp(ar2)) + * ((2 * (-2 + exp(3 * r2)) * r) / exp(4 * r2) + + sqrt(M_PI) * (erf(2. * r) - erf(r)))) + / r2; } else { @@ -211,8 +211,8 @@ double rho4(const double r, const double a) / exp(4 * r2) + exp(-a * a + 2 * a * r - 5 * r2) * (-1 + e0) * ar5 * (2 * (-2 + exp(3 * r2)) * r - + exp(4 * r2) * sqrt(M_PI) - * (erf(2. * r) - erf(r)))) + + exp(4 * r2) * sqrt(M_PI) + * (erf(2. * r) - erf(r)))) / (4. * pow(M_PI, 1.5) * r2); } else diff --git a/tests/PinnedH2O_3DOF/aux/rotation_test.cc b/tests/PinnedH2O_3DOF/aux/rotation_test.cc index 21cfd93a..3f2ae1c6 100644 --- a/tests/PinnedH2O_3DOF/aux/rotation_test.cc +++ b/tests/PinnedH2O_3DOF/aux/rotation_test.cc @@ -1,26 +1,33 @@ -#include -#include #include +#include #include +#include using namespace std; double calculate_bondlength(const double atom1[3], const double atom2[3]) { - return sqrt(pow(atom1[0] - atom2[0], 2) + pow(atom1[1] - atom2[1], 2) + pow(atom1[2] - atom2[2], 2)); + return sqrt(pow(atom1[0] - atom2[0], 2) + pow(atom1[1] - atom2[1], 2) + + pow(atom1[2] - atom2[2], 2)); } -double calculate_bondangle(const double atom1[3], const double atom2[3], const double atom3[3], bool radian) +double calculate_bondangle(const double atom1[3], const double atom2[3], + const double atom3[3], bool radian) { - double vector1[3] = {atom1[0] - atom2[0], atom1[1] - atom2[1], atom1[2] - atom2[2]}; - double vector2[3] = {atom3[0] - atom2[0], atom3[1] - atom2[1], atom3[2] - atom2[2]}; - - double dot_product = vector1[0] * vector2[0] + vector1[1] * vector2[1] + vector1[2] * vector2[2]; - double magnitude_product = sqrt(pow(vector1[0], 2) + pow(vector1[1], 2) + pow(vector1[2], 2)) * - sqrt(pow(vector2[0], 2) + pow(vector2[1], 2) + pow(vector2[2], 2)); + double vector1[3] + = { atom1[0] - atom2[0], atom1[1] - atom2[1], atom1[2] - atom2[2] }; + double vector2[3] + = { atom3[0] - atom2[0], atom3[1] - atom2[1], atom3[2] - atom2[2] }; + + double dot_product = vector1[0] * vector2[0] + vector1[1] * vector2[1] + + vector1[2] * vector2[2]; + double magnitude_product + = sqrt(pow(vector1[0], 2) + pow(vector1[1], 2) + pow(vector1[2], 2)) + * sqrt(pow(vector2[0], 2) + pow(vector2[1], 2) + pow(vector2[2], 2)); double angle = acos(dot_product / magnitude_product); - if (!radian) { + if (!radian) + { angle = angle * 180.0 / M_PI; } return angle; @@ -60,29 +67,37 @@ void cross(const double a[3], const double b[3], double result[3]) result[2] = a[0] * b[1] - a[1] * b[0]; } -void apply_rotation(const double matrix[3][3], const double vec[3], double result[3]) +void apply_rotation( + const double matrix[3][3], const double vec[3], double result[3]) { - result[0] = matrix[0][0] * vec[0] + matrix[0][1] * vec[1] + matrix[0][2] * vec[2]; - result[1] = matrix[1][0] * vec[0] + matrix[1][1] * vec[1] + matrix[1][2] * vec[2]; - result[2] = matrix[2][0] * vec[0] + matrix[2][1] * vec[1] + matrix[2][2] * vec[2]; + result[0] + = matrix[0][0] * vec[0] + matrix[0][1] * vec[1] + matrix[0][2] * vec[2]; + result[1] + = matrix[1][0] * vec[0] + matrix[1][1] * vec[1] + matrix[1][2] * vec[2]; + result[2] + = matrix[2][0] * vec[0] + matrix[2][1] * vec[1] + matrix[2][2] * vec[2]; } -void apply_transpose_rotation(const double matrix[3][3], const double vec[3], double result[3]) +void apply_transpose_rotation( + const double matrix[3][3], const double vec[3], double result[3]) { - result[0] = matrix[0][0] * vec[0] + matrix[1][0] * vec[1] + matrix[2][0] * vec[2]; - result[1] = matrix[0][1] * vec[0] + matrix[1][1] * vec[1] + matrix[2][1] * vec[2]; - result[2] = matrix[0][2] * vec[0] + matrix[1][2] * vec[1] + matrix[2][2] * vec[2]; + result[0] + = matrix[0][0] * vec[0] + matrix[1][0] * vec[1] + matrix[2][0] * vec[2]; + result[1] + = matrix[0][1] * vec[0] + matrix[1][1] * vec[1] + matrix[2][1] * vec[2]; + result[2] + = matrix[0][2] * vec[0] + matrix[1][2] * vec[1] + matrix[2][2] * vec[2]; } int main() { - double O1[3] = {0.00, 0.00, 0.00}; - double H1[3] = {-0.45, -1.48, -0.97}; - double H2[3] = {-0.45, 1.42, -1.07}; + double O1[3] = { 0.00, 0.00, 0.00 }; + double H1[3] = { -0.45, -1.48, -0.97 }; + double H2[3] = { -0.45, 1.42, -1.07 }; double bondlength1 = calculate_bondlength(H1, O1); double bondlength2 = calculate_bondlength(H2, O1); - double bondangle = calculate_bondangle(H1, O1, H2, false); + double bondangle = calculate_bondangle(H1, O1, H2, false); cout << "Original system" << endl; cout << "H1 = (" << H1[0] << ", " << H1[1] << ", " << H1[2] << ")" << endl; @@ -97,23 +112,25 @@ int main() double plane_normal[3]; cross(H2, H1, plane_normal); normalize(plane_normal); - double target_plane_normal[3] = {0, 0, 1}; + double target_plane_normal[3] = { 0, 0, 1 }; double axis_to_align[3]; cross(plane_normal, target_plane_normal, axis_to_align); normalize(axis_to_align); - double dot_product = plane_normal[0] * target_plane_normal[0] + - plane_normal[1] * target_plane_normal[1] + - plane_normal[2] * target_plane_normal[2]; + double dot_product = plane_normal[0] * target_plane_normal[0] + + plane_normal[1] * target_plane_normal[1] + + plane_normal[2] * target_plane_normal[2]; double angle_to_align = acos(min(max(dot_product, -1.0), 1.0)); double out_of_plane_rotation_matrix[3][3]; - rotation_matrix(axis_to_align, angle_to_align, out_of_plane_rotation_matrix); + rotation_matrix( + axis_to_align, angle_to_align, out_of_plane_rotation_matrix); apply_rotation(out_of_plane_rotation_matrix, H1, H1_temp); apply_rotation(out_of_plane_rotation_matrix, H2, H2_temp); - double theta1 = atan2(H1_temp[1], H1_temp[0]); + double theta1 = atan2(H1_temp[1], H1_temp[0]); double planar_rotation_angle = -theta1 + bondangle / 360.0 * M_PI; double planar_rotation_matrix[3][3]; - rotation_matrix(target_plane_normal, planar_rotation_angle, planar_rotation_matrix); + rotation_matrix( + target_plane_normal, planar_rotation_angle, planar_rotation_matrix); apply_rotation(planar_rotation_matrix, H1_temp, H1_rotated); apply_rotation(planar_rotation_matrix, H2_temp, H2_rotated); @@ -126,11 +143,16 @@ int main() double bondlength1_rotated = calculate_bondlength(H1_rotated, O1); double bondlength2_rotated = calculate_bondlength(H2_rotated, O1); - double bondangle_rotated = calculate_bondangle(H1_rotated, O1, H2_rotated, false); - - cout << "Reference system (in z=0 plane, symmetric about x=0 axis, with longer bondlength in Q1)" << endl; - cout << "H1 = (" << H1_rotated[0] << ", " << H1_rotated[1] << ", " << H1_rotated[2] << ")" << endl; - cout << "H2 = (" << H2_rotated[0] << ", " << H2_rotated[1] << ", " << H2_rotated[2] << ")" << endl; + double bondangle_rotated + = calculate_bondangle(H1_rotated, O1, H2_rotated, false); + + cout << "Reference system (in z=0 plane, symmetric about x=0 axis, with " + "longer bondlength in Q1)" + << endl; + cout << "H1 = (" << H1_rotated[0] << ", " << H1_rotated[1] << ", " + << H1_rotated[2] << ")" << endl; + cout << "H2 = (" << H2_rotated[0] << ", " << H2_rotated[1] << ", " + << H2_rotated[2] << ")" << endl; cout << "Bondlength of O1-H1 = " << bondlength1_rotated << endl; cout << "Bondlength of O1-H2 = " << bondlength2_rotated << endl; cout << "Angle between O1-H1 and O1-H2 = " << bondangle_rotated << endl; @@ -147,16 +169,21 @@ int main() apply_transpose_rotation(planar_rotation_matrix, H1_rotated, H1_temp); apply_transpose_rotation(planar_rotation_matrix, H2_rotated, H2_temp); - apply_transpose_rotation(out_of_plane_rotation_matrix, H1_temp, H1_restored); - apply_transpose_rotation(out_of_plane_rotation_matrix, H2_temp, H2_restored); + apply_transpose_rotation( + out_of_plane_rotation_matrix, H1_temp, H1_restored); + apply_transpose_rotation( + out_of_plane_rotation_matrix, H2_temp, H2_restored); double bondlength1_restored = calculate_bondlength(H1_restored, O1); double bondlength2_restored = calculate_bondlength(H2_restored, O1); - double bondangle_restored = calculate_bondangle(H1_restored, O1, H2_restored, false); + double bondangle_restored + = calculate_bondangle(H1_restored, O1, H2_restored, false); cout << "Restored system" << endl; - cout << "H1 = (" << H1_restored[0] << ", " << H1_restored[1] << ", " << H1_restored[2] << ")" << endl; - cout << "H2 = (" << H2_restored[0] << ", " << H2_restored[1] << ", " << H2_restored[2] << ")" << endl; + cout << "H1 = (" << H1_restored[0] << ", " << H1_restored[1] << ", " + << H1_restored[2] << ")" << endl; + cout << "H2 = (" << H2_restored[0] << ", " << H2_restored[1] << ", " + << H2_restored[2] << ")" << endl; cout << "Bondlength of O1-H1 = " << bondlength1_restored << endl; cout << "Bondlength of O1-H2 = " << bondlength2_restored << endl; cout << "Angle between O1-H1 and O1-H2 = " << bondangle_restored << endl; diff --git a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc index ad9ef743..8d84309d 100644 --- a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc +++ b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc @@ -13,19 +13,19 @@ #include "MGmol.h" #include "MGmol_MPI.h" #include "MPIdata.h" -#include "mgmol_run.h" #include "PinnedH2O.h" +#include "mgmol_run.h" #ifdef MGMOL_HAS_LIBROM #include "librom.h" +#include #include +#include +#include #include #include #include -#include -#include -#include #include namespace po = boost::program_options; @@ -92,8 +92,9 @@ int main(int argc, char** argv) std::cout << "-------------------------" << std::endl; } - MGmolInterface* mgmol = new MGmol>(global_comm, - *MPIdata::sout, input_filename, lrs_filename, constraints_filename); + MGmolInterface* mgmol = new MGmol>( + global_comm, *MPIdata::sout, input_filename, lrs_filename, + constraints_filename); if (MPIdata::onpe0) { @@ -139,7 +140,8 @@ int main(int argc, char** argv) const pb::Grid& mygrid = mymesh->grid(); const pb::PEenv& myPEenv = mymesh->peenv(); - // compute energy and forces again with projected problem onto ROM subspace + // compute energy and forces again with projected problem onto ROM + // subspace const int rdim = ct.getROMOptions().num_orbbasis; if (rdim != ct.numst) { @@ -152,11 +154,11 @@ int main(int argc, char** argv) std::shared_ptr projmatrices = mgmol->getProjectedMatrices(); - ExtendedGridOrbitals orbitals("new_orbitals", mygrid, mymesh->subdivx(), - ct.numst, ct.bcWF, projmatrices.get(), nullptr, nullptr, nullptr, - nullptr); + ExtendedGridOrbitals orbitals("new_orbitals", mygrid, + mymesh->subdivx(), ct.numst, ct.bcWF, projmatrices.get(), nullptr, + nullptr, nullptr, nullptr); - orbitals.set(ct.getROMOptions().basis_file, ct.numst); + orbitals.set(ct.getROMOptions().basis_file, ct.numst); orbitals.orthonormalizeLoewdin(); orbitals.setDataWithGhosts(true); @@ -253,4 +255,4 @@ int main(int argc, char** argv) return 0; } -#endif // MGMOL_HAS_LIBROM +#endif // MGMOL_HAS_LIBROM diff --git a/tests/testSetGhostValues.cc b/tests/testSetGhostValues.cc index 90aafde0..086bf4c6 100644 --- a/tests/testSetGhostValues.cc +++ b/tests/testSetGhostValues.cc @@ -27,13 +27,13 @@ TEST_CASE("Set ghost values", "[set ghosts") // prepare 3 mesh sizes to test with std::vector> meshes; { - std::array ngpts1{ { 16, 24, 20 } }; + std::array ngpts1 { { 16, 24, 20 } }; meshes.push_back(ngpts1); - std::array ngpts2{ { 20, 32, 24 } }; + std::array ngpts2 { { 20, 32, 24 } }; meshes.push_back(ngpts2); - std::array ngpts3{ { 24, 20, 32 } }; + std::array ngpts3 { { 24, 20, 32 } }; meshes.push_back(ngpts3); } diff --git a/tests/testTradeGhostValues.cc b/tests/testTradeGhostValues.cc index 8294aaec..00586d12 100644 --- a/tests/testTradeGhostValues.cc +++ b/tests/testTradeGhostValues.cc @@ -41,13 +41,13 @@ TEST_CASE("Trade ghost values", "[trade]") // prepare 3 mesh sizes to test with std::vector> meshes; { - std::array ngpts1{ { 32, 24, 20 } }; + std::array ngpts1 { { 32, 24, 20 } }; meshes.push_back(ngpts1); - std::array ngpts2{ { 20, 32, 24 } }; + std::array ngpts2 { { 20, 32, 24 } }; meshes.push_back(ngpts2); - std::array ngpts3{ { 24, 20, 32 } }; + std::array ngpts3 { { 24, 20, 32 } }; meshes.push_back(ngpts3); } @@ -187,8 +187,8 @@ TEST_CASE("Trade ghost values", "[trade]") double ref_val = (i + 1) * (uvalue - + cos3(ix - nghosts, iy - nghosts, - iz - nghosts, nx, ny, nz)); + + cos3(ix - nghosts, iy - nghosts, + iz - nghosts, nx, ny, nz)); CHECK(uu[iiy + iz] == Approx(ref_val).epsilon(1.e-8));