From 9080ac27c4e54a6cab4d64522a6adc840a01abb0 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Mon, 3 Jun 2024 22:33:59 -0400 Subject: [PATCH 01/17] Reenable testShortSighted test (#248) --- tests/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4c68adb7..fe7b571b 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -482,7 +482,7 @@ add_test(NAME ChebyshevMVP ${CMAKE_CURRENT_SOURCE_DIR}/Chebyshev/coords.in ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) -if(NOT ${MAGMA_FOUND}) +if(NOT ${MGMOL_WITH_MAGMA}) add_test(NAME testShortSighted COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/ShortSighted/test.py ${MPIEXEC} --oversubscribe ${MPIEXEC_NUMPROC_FLAG} 5 ${MPIEXEC_PREFLAGS} From f13236f18decd2e7eb62d2a20f93bca86e0e0de3 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Mon, 3 Jun 2024 22:34:35 -0400 Subject: [PATCH 02/17] Make new driver to check input (#247) * clean up/reorganize main.cc * use shared_ptr in class MGmol --- CMakeLists.txt | 4 +- drivers/CMakeLists.txt | 9 ++ drivers/check_input.cc | 93 ++++++++++++++++++ src/CMakeLists.txt | 17 ++-- src/Electrostatic.cc | 4 +- src/MGmol.cc | 208 +++++++++++++++++------------------------ src/MGmol.h | 48 +++++----- src/MGmol_prototypes.h | 2 +- src/computeHij.cc | 11 ++- src/lbfgsrlx.cc | 7 +- src/main.cc | 62 +++--------- src/md.cc | 19 ++-- src/mlwf.cc | 6 +- src/quench.cc | 47 +++++----- src/readInput.cc | 4 +- src/read_config.cc | 7 +- src/runfire.cc | 11 +-- src/setup.cc | 14 ++- 18 files changed, 310 insertions(+), 263 deletions(-) create mode 100644 drivers/CMakeLists.txt create mode 100644 drivers/check_input.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index aee0ddb5..753a455d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -162,7 +162,7 @@ if(${MGMOL_WITH_CLANG_FORMAT}) find_package(CLANG_FORMAT) if(${CLANG_FORMAT_FOUND}) message(STATUS "Indent with clang-format") - file(GLOB_RECURSE FORMAT_SOURCES src/*.cc src/*.h tests/*.cc tests/*.h) + file(GLOB_RECURSE FORMAT_SOURCES src/*.cc src/*.h tests/*.cc tests/*.h drivers/*.cc) add_custom_target(format COMMAND ${CLANG_FORMAT_EXECUTABLE} -i -style=file ${FORMAT_SOURCES} DEPENDS ${FORMAT_SOURCES}) @@ -252,5 +252,7 @@ include_directories("${PROJECT_SOURCE_DIR}/src") # add subdirectories for source files, tests add_subdirectory(src) +add_subdirectory(drivers) + add_subdirectory(tests) diff --git a/drivers/CMakeLists.txt b/drivers/CMakeLists.txt new file mode 100644 index 00000000..d0797ba0 --- /dev/null +++ b/drivers/CMakeLists.txt @@ -0,0 +1,9 @@ +include_directories( ${CMAKE_SOURCE_DIR}/src ) + +add_executable(check_input + check_input.cc + ) + +target_include_directories(check_input PRIVATE ${Boost_INCLUDE_DIRS}) + +target_link_libraries(check_input mgmol_src) diff --git a/drivers/check_input.cc b/drivers/check_input.cc new file mode 100644 index 00000000..14dbaa09 --- /dev/null +++ b/drivers/check_input.cc @@ -0,0 +1,93 @@ +// Copyright (c) 2017, Lawrence Livermore National Security, LLC and +// UT-Battelle, LLC. +// Produced at the Lawrence Livermore National Laboratory and the Oak Ridge +// National Laboratory. +// LLNL-CODE-743438 +// All rights reserved. +// 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 "Control.h" +#include "ExtendedGridOrbitals.h" +#include "LocGridOrbitals.h" +#include "MGmol.h" +#include "MGmol_MPI.h" +#include "MPIdata.h" +#include "mgmol_run.h" + +#include +#include + +#include +namespace po = boost::program_options; + +int main(int argc, char** argv) +{ + int mpirc = MPI_Init(&argc, &argv); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Initialization failed!!!" << std::endl; + MPI_Abort(MPI_COMM_WORLD, 0); + } + + MPI_Comm comm = MPI_COMM_WORLD; + + mgmol_init(comm); + + // read runtime parameters + std::string input_filename(""); + std::string lrs_filename; + std::string constraints_filename(""); + + float total_spin = 0.; + bool with_spin = false; + + po::variables_map vm; + + // read options from PE0 only + if (MPIdata::onpe0) + { + read_config(argc, argv, vm, input_filename, lrs_filename, + constraints_filename, total_spin, with_spin); + } + + MGmol_MPI::setup(comm, std::cout, with_spin); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + MPI_Comm global_comm = mmpi.commGlobal(); + + Control::setup(global_comm, with_spin, total_spin); + Control& ct = *(Control::instance()); + + ct.setOptions(vm); + + int ret = ct.checkOptions(); + if (ret < 0) return ret; + + mmpi.bcastGlobal(input_filename); + mmpi.bcastGlobal(lrs_filename); + + // Enter main scope + { + MGmolInterface* mgmol; + if (ct.isLocMode()) + 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); + + *MPIdata::sout << " Input parameters OK\n"; + + delete mgmol; + } // close main scope + + mgmol_finalize(); + + mpirc = MPI_Finalize(); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Finalize failed!!!" << std::endl; + } + + return 0; +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b7aec664..3e1188fa 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -177,13 +177,15 @@ install(TARGETS mgmol_src DESTINATION lib) add_executable(mgmol-opt main.cc) target_include_directories (mgmol-opt PRIVATE ${Boost_INCLUDE_DIRS}) -target_link_libraries(mgmol-opt mgmol_src ${link_libs}) -target_link_libraries(mgmol-opt ${SCALAPACK_LIBRARIES}) -target_link_libraries(mgmol-opt ${HDF5_LIBRARIES}) -target_link_libraries(mgmol-opt ${HDF5_HL_LIBRARIES}) -target_link_libraries(mgmol-opt ${BLAS_LIBRARIES}) -target_link_libraries(mgmol-opt ${LAPACK_LIBRARIES}) -target_link_libraries(mgmol-opt ${Boost_LIBRARIES}) +target_link_libraries(mgmol_src ${link_libs}) +target_link_libraries(mgmol_src ${SCALAPACK_LIBRARIES}) +target_link_libraries(mgmol_src ${LAPACK_LIBRARIES}) +target_link_libraries(mgmol_src ${BLAS_LIBRARIES}) +target_link_libraries(mgmol_src ${HDF5_LIBRARIES}) +target_link_libraries(mgmol_src ${HDF5_HL_LIBRARIES}) +target_link_libraries(mgmol_src ${Boost_LIBRARIES}) + +target_link_libraries(mgmol-opt mgmol_src) if (${OPENMP_CXX_FOUND}) target_link_libraries(mgmol-opt OpenMP::OpenMP_CXX) endif() @@ -192,4 +194,3 @@ if(${MGMOL_WITH_LIBXC}) endif (${MGMOL_WITH_LIBXC}) install(TARGETS mgmol-opt DESTINATION bin) - diff --git a/src/Electrostatic.cc b/src/Electrostatic.cc index b0afeb76..80cb754b 100644 --- a/src/Electrostatic.cc +++ b/src/Electrostatic.cc @@ -36,7 +36,7 @@ Timer Electrostatic::solve_tm_("Electrostatic::solve"); Electrostatic::Electrostatic(PoissonFDtype lap_type, const short bcPoisson[3], const double screening_const) - : laptype_(lap_type) + : laptype_(lap_type), poisson_solver_(nullptr) { assert(bcPoisson[0] >= 0); assert(bcPoisson[1] >= 0); @@ -166,6 +166,8 @@ Electrostatic::Electrostatic(PoissonFDtype lap_type, const short bcPoisson[3], Electrostatic::~Electrostatic() { + assert(poisson_solver_ != nullptr); + delete poisson_solver_; if (grhod_ != nullptr) delete grhod_; if (grhoc_ != nullptr) delete grhoc_; diff --git a/src/MGmol.cc b/src/MGmol.cc index b6a95777..85acc78d 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -113,73 +113,29 @@ template MGmol::MGmol(MPI_Comm comm, std::ostream& os, std::string input_filename, std::string lrs_filename, std::string constraints_filename) - : os_(os) + : os_(os), comm_(comm) { - comm_ = comm; - - constraints_ = new ConstraintSet(); + constraints_.reset(new ConstraintSet()); mgmol::out = &os; - geom_optimizer_ = nullptr; - local_cluster_ = nullptr; - proj_matrices_ = nullptr; - dm_strategy_ = nullptr; - - h5f_file_ = nullptr; - - aomm_ = nullptr; - - spreadf_ = nullptr; - - spread_penalty_ = nullptr; - - orbitals_precond_ = nullptr; - - forces_ = nullptr; - - energy_ = nullptr; + current_orbitals_ = nullptr; setupFromInput(input_filename); + /* + * Extra setup if using localization regions + */ setupLRs(lrs_filename); - setupConstraintsFromInput(constraints_filename); - - setup(); + if (!constraints_filename.empty()) + setupConstraintsFromInput(constraints_filename); } template MGmol::~MGmol() { - delete electrostat_; - delete rho_; - delete constraints_; - delete xcongrid_; - delete energy_; - if (hamiltonian_ != nullptr) delete hamiltonian_; - if (geom_optimizer_ != nullptr) delete geom_optimizer_; - - if (currentMasks_ != nullptr) delete currentMasks_; - if (corrMasks_ != nullptr) delete corrMasks_; - - if (aomm_ != nullptr) delete aomm_; - - delete current_orbitals_; - delete ions_; - delete g_kbpsi_; - - delete proj_matrices_; - if (local_cluster_ != nullptr) delete local_cluster_; - - if (h5f_file_ != nullptr) delete h5f_file_; - - if (spreadf_ != nullptr) delete spreadf_; - - if (spread_penalty_ != nullptr) delete spread_penalty_; - - delete forces_; - if (dm_strategy_ != nullptr) delete dm_strategy_; + if (current_orbitals_) delete current_orbitals_; } template <> @@ -192,18 +148,17 @@ void MGmol::initialMasks() if (ct.verbose > 0) printWithTimeStamp("MGmol::initialMasks()...", os_); - currentMasks_ = new MasksSet(false, ct.getMGlevels()); + currentMasks_ + = std::shared_ptr(new MasksSet(false, ct.getMGlevels())); currentMasks_->setup(lrs_); - corrMasks_ = new MasksSet(true, 0); + corrMasks_ = std::shared_ptr(new MasksSet(true, 0)); corrMasks_->setup(lrs_); } template <> void MGmol::initialMasks() { - currentMasks_ = nullptr; - corrMasks_ = nullptr; } template @@ -227,7 +182,7 @@ int MGmol::initial() pb::Lap* lapop = ct.Mehrstellen() ? hamiltonian_->lapOper() : nullptr; - g_kbpsi_ = new KBPsiMatrixSparse(lapop); + g_kbpsi_.reset(new KBPsiMatrixSparse(lapop)); check_anisotropy(); @@ -269,7 +224,7 @@ int MGmol::initial() // set number of iterations to 10. if (ct.load_balancing_alpha > 0.0) { - local_cluster_ = new ClusterOrbitals(lrs_); + local_cluster_.reset(new ClusterOrbitals(lrs_)); local_cluster_->setup(); local_cluster_->computeClusters(ct.load_balancing_max_iterations); } @@ -289,29 +244,31 @@ int MGmol::initial() { #ifdef HAVE_MAGMA if (use_replicated_matrix) - proj_matrices_ = new ProjectedMatricesMehrstellen( - ct.numst, with_spin, ct.occ_width); + proj_matrices_.reset( + new ProjectedMatricesMehrstellen( + ct.numst, with_spin, ct.occ_width)); else #endif - proj_matrices_ = new ProjectedMatricesMehrstellen< + proj_matrices_.reset(new ProjectedMatricesMehrstellen< dist_matrix::DistMatrix>( - ct.numst, with_spin, ct.occ_width); + ct.numst, with_spin, ct.occ_width)); } else if (ct.short_sighted) - proj_matrices_ = new ProjectedMatricesSparse( - ct.numst, ct.occ_width, lrs_, local_cluster_); + proj_matrices_.reset(new ProjectedMatricesSparse( + ct.numst, ct.occ_width, lrs_, local_cluster_.get())); else #ifdef HAVE_MAGMA if (use_replicated_matrix) - proj_matrices_ = new ProjectedMatrices( - ct.numst, with_spin, ct.occ_width); + proj_matrices_.reset(new ProjectedMatrices( + ct.numst, with_spin, ct.occ_width)); else #endif - proj_matrices_ - = new ProjectedMatrices>( - ct.numst, with_spin, ct.occ_width); + proj_matrices_.reset( + new ProjectedMatrices>( + ct.numst, with_spin, ct.occ_width)); - forces_ = new Forces(hamiltonian_, rho_, proj_matrices_); + forces_.reset(new Forces( + hamiltonian_.get(), rho_.get(), proj_matrices_.get())); if (ct.verbose > 0) printWithTimeStamp( @@ -326,8 +283,8 @@ int MGmol::initial() printWithTimeStamp("MGmol::initial(), create T...", os_); current_orbitals_ = new OrbitalsType("Primary", mygrid, mymesh->subdivx(), - ct.numst, ct.bcWF, proj_matrices_, lrs_, currentMasks_, corrMasks_, - local_cluster_, true); + ct.numst, ct.bcWF, proj_matrices_.get(), lrs_, currentMasks_.get(), + corrMasks_.get(), local_cluster_.get(), true); increaseMemorySlotsForOrbitals(); @@ -372,7 +329,7 @@ int MGmol::initial() "MGmol::initial(), init wf and masks...", os_); // Make temp mask for initial random wave functions - if (ct.init_loc == 1 && currentMasks_ != nullptr) + if (ct.init_loc == 1 && currentMasks_) { float cut_init = ct.initRadius(); assert(cut_init > 0.); @@ -383,7 +340,7 @@ int MGmol::initial() current_orbitals_->initWF(lrs_); // initialize masks again - if (ct.init_loc == 1 && currentMasks_ != nullptr) + if (ct.init_loc == 1 && currentMasks_) { currentMasks_->initialize(lrs_, 0); } @@ -427,9 +384,10 @@ int MGmol::initial() } if (ct.verbose > 0) printWithTimeStamp("Initialize XC functional...", os_); - xcongrid_ = XCfunctionalFactory::create( - ct.xctype, mmpi.nspin(), *rho_, pot); - assert(xcongrid_ != nullptr); + xcongrid_ + = std::shared_ptr(XCfunctionalFactory::create( + ct.xctype, mmpi.nspin(), *rho_, pot)); + assert(xcongrid_); // initialize nl potentials with restart values if possible // if( ct.restart_info>1 ) @@ -473,7 +431,7 @@ int MGmol::initial() { Vector3D origin(mygrid.origin(0), mygrid.origin(1), mygrid.origin(2)); Vector3D ll(mygrid.ll(0), mygrid.ll(1), mygrid.ll(2)); - spreadf_ = new SpreadsAndCenters(origin, ll); + spreadf_.reset(new SpreadsAndCenters(origin, ll)); } bool energy_with_spread_penalty = false; @@ -481,26 +439,30 @@ int MGmol::initial() { if (ct.isSpreadFunctionalVolume()) { - spread_penalty_ = new SpreadPenaltyVolume(spreadf_, - ct.spreadPenaltyTarget(), ct.spreadPenaltyAlphaFactor(), - ct.spreadPenaltyDampingFactor()); + spread_penalty_.reset( + new SpreadPenaltyVolume(spreadf_.get(), + ct.spreadPenaltyTarget(), ct.spreadPenaltyAlphaFactor(), + ct.spreadPenaltyDampingFactor())); } else if (ct.isSpreadFunctionalEnergy()) { energy_with_spread_penalty = true; - spread_penalty_ = new EnergySpreadPenalty(spreadf_, - ct.spreadPenaltyTarget(), ct.spreadPenaltyAlphaFactor()); + spread_penalty_.reset( + new EnergySpreadPenalty(spreadf_.get(), + ct.spreadPenaltyTarget(), ct.spreadPenaltyAlphaFactor())); } else - spread_penalty_ = new SpreadPenalty(spreadf_, - ct.spreadPenaltyTarget(), ct.spreadPenaltyAlphaFactor(), - ct.spreadPenaltyDampingFactor()); + spread_penalty_.reset( + new SpreadPenalty(spreadf_.get(), + ct.spreadPenaltyTarget(), ct.spreadPenaltyAlphaFactor(), + ct.spreadPenaltyDampingFactor())); } - SpreadPenaltyInterface* spread_penalty + std::shared_ptr> spread_penalty = energy_with_spread_penalty ? spread_penalty_ : nullptr; - energy_ = new Energy( - mygrid, *ions_, pot, *electrostat_, *rho_, *xcongrid_, spread_penalty); + energy_ = std::shared_ptr>( + new Energy(mygrid, *ions_, pot, *electrostat_, *rho_, + *xcongrid_, spread_penalty.get())); if (ct.verbose > 0) printWithTimeStamp("Setup matrices...", os_); @@ -509,15 +471,16 @@ int MGmol::initial() // HMVP algorithm requires that H is initialized #ifdef HAVE_MAGMA if (use_replicated_matrix) - dm_strategy_ - = DMStrategyFactory::create(comm_, - os_, *ions_, rho_, energy_, electrostat_, this, proj_matrices_, - current_orbitals_); + dm_strategy_.reset( + DMStrategyFactory::create(comm_, + os_, *ions_, rho_.get(), energy_.get(), electrostat_.get(), + this, proj_matrices_.get(), current_orbitals_)); else #endif - dm_strategy_ = DMStrategyFactory>::create(comm_, os_, *ions_, rho_, - energy_, electrostat_, this, proj_matrices_, current_orbitals_); + dm_strategy_.reset(DMStrategyFactory>::create(comm_, os_, *ions_, + rho_.get(), energy_.get(), electrostat_.get(), this, + proj_matrices_.get(), current_orbitals_)); // theta = invB * Hij proj_matrices_->updateThetaAndHB(); @@ -609,7 +572,7 @@ void MGmol::finalEnergy() // Get the total energy const double ts = 0.5 * proj_matrices_->computeEntropy(); // in [Ha] total_energy_ = energy_->evaluateTotal( - ts, proj_matrices_, *current_orbitals_, 2, os_); + ts, proj_matrices_.get(), *current_orbitals_, 2, os_); } template @@ -624,7 +587,7 @@ void MGmol::printMM() ProjectedMatrices>* projmatrices = dynamic_cast< ProjectedMatrices>*>( - proj_matrices_); + proj_matrices_.get()); assert(projmatrices != nullptr); projmatrices->printHamiltonianMM(tfileh); } @@ -697,7 +660,7 @@ void MGmol::write_header() ct.printPoissonOptions(os_); } // onpe0 - if (current_orbitals_ != nullptr && ct.verbose > 0) + if (current_orbitals_ && ct.verbose > 0) { current_orbitals_->printNumst(os_); current_orbitals_->printChromaticNumber(os_); @@ -802,9 +765,9 @@ void MGmol::printEigAndOcc() #ifdef HAVE_MAGMA // try with ReplicatedMatrix first { - ProjectedMatrices* projmatrices - = dynamic_cast*>( - proj_matrices_); + std::shared_ptr> projmatrices + = std::dynamic_pointer_cast< + ProjectedMatrices>(proj_matrices_); if (projmatrices) { projmatrices->printEigenvalues(os_); @@ -815,10 +778,10 @@ void MGmol::printEigAndOcc() #endif if (!printflag) { - ProjectedMatrices>* - projmatrices - = dynamic_cast< - ProjectedMatrices>*>( + std::shared_ptr< + ProjectedMatrices>> + projmatrices = std::dynamic_pointer_cast< + ProjectedMatrices>>( proj_matrices_); assert(projmatrices); @@ -1062,21 +1025,23 @@ double MGmol::get_evnl(const Ions& ions) double val; if (ct.short_sighted) { - ProjectedMatricesSparse* projmatrices - = dynamic_cast(proj_matrices_); + std::shared_ptr projmatrices + = std::dynamic_pointer_cast( + proj_matrices_); assert(projmatrices); - val = g_kbpsi_->getEvnl(ions, projmatrices); + val = g_kbpsi_->getEvnl(ions, projmatrices.get()); } else { - ProjectedMatrices>* projmatrices - = dynamic_cast< - ProjectedMatrices>*>( + std::shared_ptr< + ProjectedMatrices>> + projmatrices = std::dynamic_pointer_cast< + ProjectedMatrices>>( proj_matrices_); assert(projmatrices); - val = g_kbpsi_->getEvnl(ions, projmatrices); + val = g_kbpsi_->getEvnl(ions, projmatrices.get()); } evnl_tm_.stop(); @@ -1102,11 +1067,11 @@ void MGmol::setup() printWithTimeStamp("MGmol::setup()...", os_); if (ct.verbose > 0) printWithTimeStamp("Setup VH...", os_); - electrostat_ = new Electrostatic( - ct.getPoissonFDtype(), ct.bcPoisson, ct.screening_const); + electrostat_ = std::shared_ptr(new Electrostatic( + ct.getPoissonFDtype(), ct.bcPoisson, ct.screening_const)); electrostat_->setup(ct.vh_init); - rho_ = new Rho(); + rho_ = std::shared_ptr>(new Rho()); rho_->setVerbosityLevel(ct.verbose); #ifdef HAVE_MAGMA @@ -1193,17 +1158,18 @@ template void MGmol::setGamma( const pb::Lap& lapOper, const Potentials& pot) { - assert(orbitals_precond_ != nullptr); + assert(orbitals_precond_); Control& ct = *(Control::instance()); - orbitals_precond_->setGamma(lapOper, pot, ct.getMGlevels(), proj_matrices_); + orbitals_precond_->setGamma( + lapOper, pot, ct.getMGlevels(), proj_matrices_.get()); } template void MGmol::precond_mg(OrbitalsType& phi) { - assert(orbitals_precond_ != nullptr); + assert(orbitals_precond_); orbitals_precond_->precond_mg(phi); } @@ -1434,7 +1400,7 @@ template void MGmol::addResidualSpreadPenalty( OrbitalsType& phi, OrbitalsType& res) { - assert(spread_penalty_ != nullptr); + assert(spread_penalty_); spread_penalty_->addResidual(phi, res); } diff --git a/src/MGmol.h b/src/MGmol.h index a44a37f8..ceec7390 100644 --- a/src/MGmol.h +++ b/src/MGmol.h @@ -42,6 +42,7 @@ class IonicAlgorithm; #include "AOMMprojector.h" #include "ClusterOrbitals.h" #include "DMStrategy.h" +#include "Energy.h" #include "ExtendedGridOrbitals.h" #include "Forces.h" #include "Ions.h" @@ -53,6 +54,8 @@ class IonicAlgorithm; #include "SpreadPenaltyInterface.h" #include "SpreadsAndCenters.h" +#include + template class MGmol : public MGmolInterface { @@ -61,48 +64,49 @@ class MGmol : public MGmolInterface MPI_Comm comm_; - XConGrid* xcongrid_; + std::shared_ptr xcongrid_; OrbitalsType* current_orbitals_; - AOMMprojector* aomm_; + std::shared_ptr aomm_; - Ions* ions_; + std::shared_ptr ions_; - Rho* rho_; + std::shared_ptr> rho_; - Energy* energy_; + std::shared_ptr> energy_; - Hamiltonian* hamiltonian_; + std::shared_ptr> hamiltonian_; - Forces* forces_; + std::shared_ptr> forces_; - MasksSet* currentMasks_; - MasksSet* corrMasks_; + std::shared_ptr currentMasks_; + std::shared_ptr corrMasks_; - // ProjectedMatrices* proj_matrices_; - ProjectedMatricesInterface* proj_matrices_; + std::shared_ptr proj_matrices_; - IonicAlgorithm* geom_optimizer_; + std::shared_ptr> geom_optimizer_; std::shared_ptr lrs_; - ClusterOrbitals* local_cluster_; + std::shared_ptr local_cluster_; - KBPsiMatrixSparse* g_kbpsi_; + std::shared_ptr g_kbpsi_; - SpreadsAndCenters* spreadf_; + std::shared_ptr> spreadf_; - SpreadPenaltyInterface* spread_penalty_; + std::shared_ptr> spread_penalty_; - DMStrategy* dm_strategy_; + std::shared_ptr> dm_strategy_; - HDFrestart* h5f_file_; + std::shared_ptr h5f_file_; + + std::shared_ptr> orbitals_precond_; double total_energy_; - ConstraintSet* constraints_; + std::shared_ptr constraints_; - OrbitalsExtrapolation* orbitals_extrapol_ = nullptr; + std::shared_ptr> orbitals_extrapol_; float md_time_; int md_iteration_; @@ -166,10 +170,8 @@ class MGmol : public MGmolInterface static Timer comp_res_tm_; static Timer init_nuc_tm_; - OrbitalsPreconditioning* orbitals_precond_; - public: - Electrostatic* electrostat_; + std::shared_ptr electrostat_; MGmol(MPI_Comm comm, std::ostream& os, std::string input_filename, std::string lrs_filename, std::string constraints_filename); diff --git a/src/MGmol_prototypes.h b/src/MGmol_prototypes.h index 8c47c388..cedd514a 100644 --- a/src/MGmol_prototypes.h +++ b/src/MGmol_prototypes.h @@ -23,5 +23,5 @@ double getLAeigen(const double tol, const int maxit, Ions& ions); int read_config(int argc, char** argv, boost::program_options::variables_map& vm, std::string& input_file, std::string& lrs_filename, std::string& constraints_filename, - float& total_spin, bool& with_spin, bool& tcheck); + float& total_spin, bool& with_spin); #endif diff --git a/src/computeHij.cc b/src/computeHij.cc index 4942959a..75b8ef86 100644 --- a/src/computeHij.cc +++ b/src/computeHij.cc @@ -280,13 +280,14 @@ template void MGmol::getKBPsiAndHij(OrbitalsType& orbitals, Ions& ions, KBPsiMatrixSparse* kbpsi, dist_matrix::DistMatrix& hij) { - getKBPsiAndHij(orbitals, orbitals, ions, kbpsi, proj_matrices_, hij); + getKBPsiAndHij(orbitals, orbitals, ions, kbpsi, proj_matrices_.get(), hij); } template void MGmol::getKBPsiAndHij(OrbitalsType& orbitals, Ions& ions) { - getKBPsiAndHij(orbitals, orbitals, ions, g_kbpsi_, proj_matrices_); + getKBPsiAndHij( + orbitals, orbitals, ions, g_kbpsi_.get(), proj_matrices_.get()); } template @@ -396,7 +397,7 @@ template void MGmol::getHpsiAndTheta( Ions& ions, OrbitalsType& phi, OrbitalsType& hphi) { - getHpsiAndTheta(ions, phi, hphi, g_kbpsi_); + getHpsiAndTheta(ions, phi, hphi, g_kbpsi_.get()); } template @@ -434,10 +435,10 @@ void MGmol::getHpsiAndTheta(Ions& ions, OrbitalsType& phi, proj_matrices_->clearSparseH(); // Compute the contribution of the non-local potential into sh - kbpsi->computeHvnlMatrix(ions, proj_matrices_); + kbpsi->computeHvnlMatrix(ions, proj_matrices_.get()); // add local part of H to sh - hamiltonian_->addHlocalij(phi, proj_matrices_); + hamiltonian_->addHlocalij(phi, proj_matrices_.get()); energy_->saveVofRho(); diff --git a/src/lbfgsrlx.cc b/src/lbfgsrlx.cc index 3881abab..b2aff795 100644 --- a/src/lbfgsrlx.cc +++ b/src/lbfgsrlx.cc @@ -36,15 +36,14 @@ void MGmol::lbfgsrlx(OrbitalsType** orbitals, Ions& ions) Control& ct = *(Control::instance()); LBFGS lbfgs(orbitals, ions, *rho_, *constraints_, lrs_, - local_cluster_, *currentMasks_, *corrMasks_, *electrostat_, ct.dt, + local_cluster_.get(), *currentMasks_, *corrMasks_, *electrostat_, ct.dt, *this); DFTsolver::resetItCount(); - lbfgs.init(h5f_file_); + lbfgs.init(h5f_file_.get()); - delete h5f_file_; - h5f_file_ = nullptr; + h5f_file_.reset(); // additional quench to compensate random start if (ct.restart_info < 3) diff --git a/src/main.cc b/src/main.cc index f6f396d2..3fcbf791 100644 --- a/src/main.cc +++ b/src/main.cc @@ -29,39 +29,16 @@ #include "MGmol.h" #include "MGmol_MPI.h" #include "MPIdata.h" -#include "Mesh.h" #include "mgmol_run.h" -#include "tools.h" #include -#include #include -#include -#include #include #include -#include - #include namespace po = boost::program_options; -//#include "MemTrack.h" - -/* -void trapfpe () { - feenableexcept(FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW); -} -*/ - -// A helper function -template -std::ostream& operator<<(std::ostream& os, const std::vector& v) -{ - copy(v.begin(), v.end(), std::ostream_iterator(std::cout, " ")); - return os; -} - int main(int argc, char** argv) { int mpirc = MPI_Init(&argc, &argv); @@ -73,33 +50,37 @@ int main(int argc, char** argv) MPI_Comm comm = MPI_COMM_WORLD; + /* + * Initialize general things, like magma, openmp, IO, ... + */ mgmol_init(comm); - // read runtime parameters + /* + * read runtime parameters + */ std::string input_filename(""); std::string lrs_filename; std::string constraints_filename(""); - bool tcheck = false; float total_spin = 0.; bool with_spin = false; po::variables_map vm; - // use configure file if it can be found - // std::string config_filename("mgmol.cfg"); - - // read options from PE0 only + // read from PE0 only if (MPIdata::onpe0) { read_config(argc, argv, vm, input_filename, lrs_filename, - constraints_filename, total_spin, with_spin, tcheck); + constraints_filename, total_spin, with_spin); } MGmol_MPI::setup(comm, std::cout, with_spin); MGmol_MPI& mmpi = *(MGmol_MPI::instance()); MPI_Comm global_comm = mmpi.commGlobal(); + /* + * Setup control struct with run time parameters + */ Control::setup(global_comm, with_spin, total_spin); Control& ct = *(Control::instance()); @@ -108,11 +89,6 @@ int main(int argc, char** argv) int ret = ct.checkOptions(); if (ret < 0) return ret; - unsigned ngpts[3] = { ct.ngpts_[0], ct.ngpts_[1], ct.ngpts_[2] }; - double origin[3] = { ct.ox_, ct.oy_, ct.oz_ }; - const double cell[3] = { ct.lx_, ct.ly_, ct.lz_ }; - Mesh::setup(mmpi.commSpin(), ngpts, origin, cell, ct.lap_type); - mmpi.bcastGlobal(input_filename); mmpi.bcastGlobal(lrs_filename); @@ -126,14 +102,10 @@ int main(int argc, char** argv) mgmol = new MGmol(global_comm, *MPIdata::sout, input_filename, lrs_filename, constraints_filename); - if (!tcheck) - { - mgmol->run(); - } - else - { - *MPIdata::sout << " Input parameters OK\n"; - } + mgmol->setup(); + + mgmol->run(); + delete mgmol; } // close main scope @@ -150,9 +122,5 @@ int main(int argc, char** argv) time(&tt); if (onpe0) std::cout << " Run ended at " << ctime(&tt) << std::endl; - // MemTrack::TrackDumpBlocks(); - - // MemTrack::TrackListMemoryUsage(); - return 0; } diff --git a/src/md.cc b/src/md.cc index c01589e5..f1011565 100644 --- a/src/md.cc +++ b/src/md.cc @@ -137,8 +137,8 @@ OrbitalsType* MGmol::new_orbitals_with_current_LRs(bool setup) // need to build new orbitals as masks have changed OrbitalsType* new_orbitals = new OrbitalsType("NewMasks", mygrid, - mymesh->subdivx(), ct.numst, ct.bcWF, proj_matrices_, lrs_, - currentMasks_, corrMasks_, local_cluster_, setup); + mymesh->subdivx(), ct.numst, ct.bcWF, proj_matrices_.get(), lrs_, + currentMasks_.get(), corrMasks_.get(), local_cluster_.get(), setup); return new_orbitals; } @@ -317,8 +317,8 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) int size_tau = (int)tau0.size(); DFTsolver::resetItCount(); - orbitals_extrapol_ = OrbitalsExtrapolationFactory::create( - ct.WFExtrapolation()); + orbitals_extrapol_.reset(OrbitalsExtrapolationFactory::create( + ct.WFExtrapolation())); MD_IonicStepper* stepper = new MD_IonicStepper( ct.dt, atmove, tau0, taup, taum, fion, pmass, rand_states); @@ -376,8 +376,8 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) if (onpe0) os_ << "Create new orbitals_minus1..." << std::endl; orbitals_extrapol_->setupPreviousOrbitals(¤t_orbitals_, - proj_matrices_, lrs_, local_cluster_, currentMasks_, - corrMasks_, *h5f_file_); + proj_matrices_.get(), lrs_, local_cluster_.get(), + currentMasks_.get(), corrMasks_.get(), *h5f_file_); // need to reset a few things as we just read new orbitals (*orbitals)->computeGramAndInvS(); @@ -395,8 +395,7 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) moveVnuc(ions); } - delete h5f_file_; - h5f_file_ = nullptr; + h5f_file_.reset(); } // additional SC steps to compensate random start @@ -437,7 +436,7 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) if (ct.adaptiveLRs()) { assert(lrs_); - adaptLR(spreadf_, nullptr); + adaptLR(spreadf_.get(), nullptr); last_move_is_small = lrs_->moveIsSmall(); @@ -670,7 +669,7 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) } delete stepper; - delete orbitals_extrapol_; + orbitals_extrapol_.reset(); } template diff --git a/src/mlwf.cc b/src/mlwf.cc index cc4fdeb8..3efdfeb2 100644 --- a/src/mlwf.cc +++ b/src/mlwf.cc @@ -279,9 +279,9 @@ int MGmol::get_NOLMO(NOLMOTransform& noot, OrbitalsType& orbitals, } sincos.clear(); - ProjectedMatrices>* projmatrices - = dynamic_cast< - ProjectedMatrices>*>( + std::shared_ptr>> + projmatrices = std::dynamic_pointer_cast< + ProjectedMatrices>>( proj_matrices_); assert(projmatrices); diff --git a/src/quench.cc b/src/quench.cc index de01dad3..614aa636 100644 --- a/src/quench.cc +++ b/src/quench.cc @@ -326,8 +326,9 @@ bool MGmol::rotateStatesPairsOverlap( spreadf2st.print(os_); orbitals.computeGramAndInvS(); - ProjectedMatricesSparse* projmatrices - = dynamic_cast(proj_matrices_); + std::shared_ptr projmatrices + = std::dynamic_pointer_cast( + proj_matrices_); if (onpe0) std::cout << "Gram Matrix for pair of states after transformation:" << std::endl; @@ -398,7 +399,7 @@ void MGmol::disentangleOrbitals(OrbitalsType& orbitals, template <> void MGmol::applyAOMMprojection(LocGridOrbitals& orbitals) { - aomm_ = new AOMMprojector(orbitals, lrs_); + aomm_.reset(new AOMMprojector(orbitals, lrs_)); aomm_->projectOut(orbitals); } @@ -422,8 +423,9 @@ int MGmol::outerSolve(LocGridOrbitals& orbitals, case OuterSolverType::ABPG: case OuterSolverType::NLCG: { - DFTsolver solver(hamiltonian_, proj_matrices_, - energy_, electrostat_, this, ions, rho_, dm_strategy_, os_); + DFTsolver solver(hamiltonian_.get(), + proj_matrices_.get(), energy_.get(), electrostat_.get(), this, + ions, rho_.get(), dm_strategy_.get(), os_); retval = solver.solve( orbitals, work_orbitals, ions, max_steps, iprint, last_eks); @@ -433,9 +435,9 @@ int MGmol::outerSolve(LocGridOrbitals& orbitals, case OuterSolverType::PolakRibiere: { - PolakRibiereSolver solver(hamiltonian_, - proj_matrices_, energy_, electrostat_, this, ions, rho_, - dm_strategy_, os_); + PolakRibiereSolver solver(hamiltonian_.get(), + proj_matrices_.get(), energy_.get(), electrostat_.get(), this, + ions, rho_.get(), dm_strategy_.get(), os_); retval = solver.solve( orbitals, work_orbitals, ions, max_steps, iprint, last_eks); @@ -466,8 +468,9 @@ int MGmol::outerSolve(OrbitalsType& orbitals, case OuterSolverType::ABPG: case OuterSolverType::NLCG: { - DFTsolver solver(hamiltonian_, proj_matrices_, - energy_, electrostat_, this, ions, rho_, dm_strategy_, os_); + DFTsolver solver(hamiltonian_.get(), + proj_matrices_.get(), energy_.get(), electrostat_.get(), this, + ions, rho_.get(), dm_strategy_.get(), os_); retval = solver.solve( orbitals, work_orbitals, ions, max_steps, iprint, last_eks); @@ -477,9 +480,9 @@ int MGmol::outerSolve(OrbitalsType& orbitals, case OuterSolverType::PolakRibiere: { - PolakRibiereSolver solver(hamiltonian_, - proj_matrices_, energy_, electrostat_, this, ions, rho_, - dm_strategy_, os_); + PolakRibiereSolver solver(hamiltonian_.get(), + proj_matrices_.get(), energy_.get(), electrostat_.get(), this, + ions, rho_.get(), dm_strategy_.get(), os_); retval = solver.solve( orbitals, work_orbitals, ions, max_steps, iprint, last_eks); @@ -499,8 +502,9 @@ int MGmol::outerSolve(OrbitalsType& orbitals, #else DavidsonSolver> #endif - solver(os_, *ions_, hamiltonian_, rho_, energy_, electrostat_, - this, gids, ct.dm_mix, with_spin); + solver(os_, *ions_, hamiltonian_.get(), rho_.get(), + energy_.get(), electrostat_.get(), this, gids, ct.dm_mix, + with_spin); retval = solver.solve(orbitals, work_orbitals); break; @@ -558,9 +562,9 @@ int MGmol::quench(OrbitalsType* orbitals, Ions& ions, applyAOMMprojection(*orbitals); } - orbitals_precond_ = new OrbitalsPreconditioning(); + orbitals_precond_.reset(new OrbitalsPreconditioning()); orbitals_precond_->setup( - *orbitals, ct.getMGlevels(), ct.lap_type, currentMasks_, lrs_); + *orbitals, ct.getMGlevels(), ct.lap_type, currentMasks_.get(), lrs_); // solve electronic structure problem // (inner iterations) @@ -570,11 +574,9 @@ int MGmol::quench(OrbitalsType* orbitals, Ions& ions, if (ct.use_kernel_functions) { - delete aomm_; - aomm_ = nullptr; + aomm_.reset(); } - delete orbitals_precond_; - orbitals_precond_ = nullptr; + orbitals_precond_.reset(); // Get the n.l. energy // TODO: Fix bug where energy vs. time output is incorrect if get_evnl is @@ -599,7 +601,8 @@ int MGmol::quench(OrbitalsType* orbitals, Ions& ions, << " TS [Ha] = " << ts << std::endl; } } - last_eks = energy_->evaluateTotal(ts, proj_matrices_, *orbitals, 2, os_); + last_eks + = energy_->evaluateTotal(ts, proj_matrices_.get(), *orbitals, 2, os_); if (ct.computeCondGramMD()) { diff --git a/src/readInput.cc b/src/readInput.cc index b0a47411..79fbecf1 100644 --- a/src/readInput.cc +++ b/src/readInput.cc @@ -175,7 +175,7 @@ int MGmol::readCoordinates( // setup ions const std::vector& sp(ct.getSpecies()); - ions_ = new Ions(lattice, sp); + ions_.reset(new Ions(lattice, sp)); if (ct.restart_info > 0 && ct.override_restart == 0) // read restart ionic positions @@ -217,7 +217,7 @@ int MGmol::readCoordinates( // setup ions const std::vector& sp(ct.getSpecies()); - ions_ = new Ions(lattice, sp); + ions_.reset(new Ions(lattice, sp)); if (ct.restart_info > 0 && ct.override_restart == 0) // read restart ionic positions diff --git a/src/read_config.cc b/src/read_config.cc index 534981d8..642b1b99 100644 --- a/src/read_config.cc +++ b/src/read_config.cc @@ -20,8 +20,7 @@ namespace po = boost::program_options; int read_config(int argc, char** argv, po::variables_map& vm, std::string& input_file, std::string& lrs_filename, - std::string& constraints_filename, float& total_spin, bool& with_spin, - bool& tcheck) + std::string& constraints_filename, float& total_spin, bool& with_spin) { // use configure file if it can be found // std::string config_filename("mgmol.cfg"); @@ -373,10 +372,6 @@ int read_config(int argc, char** argv, po::variables_map& vm, return 0; } - if (vm.count("check")) - { - tcheck = true; - } if (vm.count("spin")) { total_spin = vm["spin"].as(); diff --git a/src/runfire.cc b/src/runfire.cc index d71c52a8..3239e1b1 100644 --- a/src/runfire.cc +++ b/src/runfire.cc @@ -33,13 +33,12 @@ void MGmol::runfire(OrbitalsType** orbitals, Ions& ions) DFTsolver::resetItCount(); - orbitals_extrapol_ = OrbitalsExtrapolationFactory::create( - ct.WFExtrapolation()); + orbitals_extrapol_.reset(OrbitalsExtrapolationFactory::create( + ct.WFExtrapolation())); - fire.init(h5f_file_); + fire.init(h5f_file_.get()); - delete h5f_file_; - h5f_file_ = nullptr; + h5f_file_.reset(); // additional quench to compensate random start if (ct.restart_info < 3) @@ -129,7 +128,7 @@ void MGmol::runfire(OrbitalsType** orbitals, Ions& ions) } // end for steps - delete orbitals_extrapol_; + orbitals_extrapol_.reset(); // final dump if (ct.out_restart_info > 0) diff --git a/src/setup.cc b/src/setup.cc index 62d9f835..f022678c 100644 --- a/src/setup.cc +++ b/src/setup.cc @@ -28,7 +28,15 @@ int MGmol::setupFromInput(const std::string filename) MGmol_MPI& mmpi = *(MGmol_MPI::instance()); - hamiltonian_ = new Hamiltonian(); + /* + * Setup global mesh for calculations + */ + unsigned ngpts[3] = { ct.ngpts_[0], ct.ngpts_[1], ct.ngpts_[2] }; + double origin[3] = { ct.ox_, ct.oy_, ct.oz_ }; + const double cell[3] = { ct.lx_, ct.ly_, ct.lz_ }; + Mesh::setup(mmpi.commSpin(), ngpts, origin, cell, ct.lap_type); + + hamiltonian_.reset(new Hamiltonian()); Potentials& pot = hamiltonian_->potential(); ct.registerPotentials(pot); @@ -43,8 +51,8 @@ int MGmol::setupFromInput(const std::string filename) const pb::PEenv& myPEenv = mymesh->peenv(); if (ct.restart_info > 0) - h5f_file_ - = new HDFrestart(ct.restart_file, myPEenv, ct.restart_file_type); + h5f_file_.reset( + new HDFrestart(ct.restart_file, myPEenv, ct.restart_file_type)); int status = readCoordinates(filename, false); if (status == -1) return -1; From bad80d13fcb1a525aa65351787ec86f75d2daedc Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Tue, 4 Jun 2024 09:44:27 -0400 Subject: [PATCH 03/17] Add possible periodic dimensions to xyz2in.py (#249) * Add possible periodic dimensions to xyz2in.py --- util/xyz2in.py | 46 ++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 38 insertions(+), 8 deletions(-) diff --git a/util/xyz2in.py b/util/xyz2in.py index 3555e138..bb2e022d 100644 --- a/util/xyz2in.py +++ b/util/xyz2in.py @@ -5,9 +5,11 @@ # 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 # -# Python program to convert .xyz file into mgmol input +# Python program to convert .xyz file into mgmol input coordinates file +# Optional arguments: [lx ly lz] to define box size and map all atoms into +# box (0,0,0)-(lx,ly,lz) using periodic boundary conditions # -# use: python coords.xyz > coords.in +# use: python coords.xyz [lx ly lz] > coords.in #------------------------------------------------------------------------------- import sys, string @@ -15,6 +17,16 @@ #read file ifile=open(sys.argv[1],'r') + +lx = 0. +ly = 0. +lz = 0. +if( len(sys.argv) > 2 ): + lx = eval(sys.argv[2]) + ly = eval(sys.argv[3]) + lz = eval(sys.argv[4]) + + lines=ifile.readlines() count=0 @@ -22,13 +34,31 @@ dummy=0 #unused flag set to 0 for line in lines: ## loop over lines of file - words=line.split() - if len(words)>1: - if words[0][0:1]!='#': - name=words[0]+str(count) + if count>1: + words=line.split() + if len(words)>1: + name=words[0]+str(count-2) x=eval(words[1])*ang2bohr y=eval(words[2])*ang2bohr z=eval(words[3])*ang2bohr - + + if lx > 0.: + if x<0: + x = x +lx + if x>lx: + x = x -lx + + if ly > 0.: + if y<0: + y = y +ly + if y>ly: + y = y -ly + + if lz > 0.: + if z<0: + z = z +lz + if z>lz: + z = z -lz + print(name,'\t',dummy,'\t',x,'\t',y,'\t',z,'\t',movable) - count=count+1 + count=count+1 From 7aac379a2bec8f532772a57ba3a0aebe5b8d9cff Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Tue, 4 Jun 2024 10:18:29 -0400 Subject: [PATCH 04/17] Remove unused/untested option extrapolateH (#250) --- src/OrbitalsExtrapolation.h | 2 -- src/md.cc | 5 ----- 2 files changed, 7 deletions(-) diff --git a/src/OrbitalsExtrapolation.h b/src/OrbitalsExtrapolation.h index 47c1a21b..e2d0ff0c 100644 --- a/src/OrbitalsExtrapolation.h +++ b/src/OrbitalsExtrapolation.h @@ -26,8 +26,6 @@ class OrbitalsExtrapolation virtual ~OrbitalsExtrapolation(); - virtual bool extrapolatedH() const { return false; } - virtual void clearOldOrbitals(); bool getRestartData(OrbitalsType& orbitals); virtual void setupPreviousOrbitals(OrbitalsType** orbitals, diff --git a/src/md.cc b/src/md.cc index f1011565..41dcd39b 100644 --- a/src/md.cc +++ b/src/md.cc @@ -91,11 +91,6 @@ void MGmol::postWFextrapolation(OrbitalsType* orbitals) else proj_matrices_->setDMto2InvS(); } - else - { - if (orbitals_extrapol_->extrapolatedH()) - dm_strategy_->update(*orbitals); - } } template From bffc514d1941346cca10ee5a9dc5b5d46cfe97c0 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Thu, 6 Jun 2024 10:34:05 -0400 Subject: [PATCH 05/17] Exit with failure if density off by more than 2% (#251) * Exit with failure if density off by more than 2% * adapt SiH4 test to catch that * fix bug in DFTsolver that was leading to wrong density --- src/DFTsolver.cc | 4 +++- src/Rho.cc | 6 ++++++ tests/SiH4/testSiH4.py | 13 +++++++++---- 3 files changed, 18 insertions(+), 5 deletions(-) diff --git a/src/DFTsolver.cc b/src/DFTsolver.cc index c1484160..91c4a0f6 100644 --- a/src/DFTsolver.cc +++ b/src/DFTsolver.cc @@ -339,13 +339,15 @@ int DFTsolver::solve(OrbitalsType& orbitals, else { bool updateDM = false; + bool updatedS = false; if (!ct.fullyOccupied()) { orbitals.computeGramAndInvS(); dm_strategy_->dressDM(); updateDM = true; + updatedS = true; } - orbitals.orthonormalizeLoewdin(true, nullptr, updateDM); + orbitals.orthonormalizeLoewdin(updatedS, nullptr, updateDM); orbitals_stepper_->restartMixing(); } diff --git a/src/Rho.cc b/src/Rho.cc index b456f6c1..ca0f20ea 100644 --- a/src/Rho.cc +++ b/src/Rho.cc @@ -175,6 +175,12 @@ void Rho::rescaleTotalCharge() (*MPIdata::sout) << " Rescaling factor: " << t1 << std::endl; (*MPIdata::sout) << " Num. electrons: " << nel << std::endl; } + if (fabs(t1 - 1.) > 0.02) + { + std::cerr << "Error on density too large to continue. Abort." + << std::endl; + mmpi.abort(); + } for (int ispin = 0; ispin < nspin; ispin++) LinearAlgebraUtils::MPscal( diff --git a/tests/SiH4/testSiH4.py b/tests/SiH4/testSiH4.py index 1a0965ac..c7f81bdc 100755 --- a/tests/SiH4/testSiH4.py +++ b/tests/SiH4/testSiH4.py @@ -44,16 +44,21 @@ lines=output.split(b'\n') tol = 5.e-3 +found_forces = False for line in lines: - num_matches = line.count(b'%%') - if num_matches: + if line.count(b'%%'): print(line) - num_matches = line.count(b'##') - if num_matches: + if line.count(b'##'): words=line.split() if len(words)==8: print(line) + found_forces = True for i in range(5,8): if abs(eval(words[i]))>tol: sys.exit(1) +if (not found_forces): + print("no forces found") + sys.exit(1) + +sys.exit(0) From de1be3b23a3d5ac0f8e19bb64bc71f8b5c48ef3e Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Thu, 6 Jun 2024 11:08:38 -0400 Subject: [PATCH 06/17] Example driver (#252) * add example driver, showing use of MGmol as a force/energy computational engine * clean up related functions in class Ions --- drivers/CMakeLists.txt | 8 +- drivers/example1.cc | 167 +++++++++++++++++++++++++++++++++++ src/Ions.cc | 191 ++++++++++++++++++++--------------------- src/Ions.h | 28 +++--- src/MGmol.cc | 33 +++++++ src/MGmol.h | 15 +++- src/MGmolInterface.h | 6 ++ src/md.cc | 12 ++- 8 files changed, 340 insertions(+), 120 deletions(-) create mode 100644 drivers/example1.cc diff --git a/drivers/CMakeLists.txt b/drivers/CMakeLists.txt index d0797ba0..4eb1c741 100644 --- a/drivers/CMakeLists.txt +++ b/drivers/CMakeLists.txt @@ -1,9 +1,11 @@ include_directories( ${CMAKE_SOURCE_DIR}/src ) -add_executable(check_input - check_input.cc - ) +add_executable(check_input check_input.cc) + +add_executable(example1 example1.cc) target_include_directories(check_input PRIVATE ${Boost_INCLUDE_DIRS}) +target_include_directories(example1 PRIVATE ${Boost_INCLUDE_DIRS}) target_link_libraries(check_input mgmol_src) +target_link_libraries(example1 mgmol_src) diff --git a/drivers/example1.cc b/drivers/example1.cc new file mode 100644 index 00000000..561f5e08 --- /dev/null +++ b/drivers/example1.cc @@ -0,0 +1,167 @@ +// Copyright (c) 2017, Lawrence Livermore National Security, LLC and +// UT-Battelle, LLC. +// Produced at the Lawrence Livermore National Laboratory and the Oak Ridge +// National Laboratory. +// LLNL-CODE-743438 +// All rights reserved. +// 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 "Control.h" +#include "ExtendedGridOrbitals.h" +#include "LocGridOrbitals.h" +#include "MGmol.h" +#include "MGmol_MPI.h" +#include "MPIdata.h" +#include "mgmol_run.h" + +#include +#include +#include +#include + +#include +namespace po = boost::program_options; + +int main(int argc, char** argv) +{ + int mpirc = MPI_Init(&argc, &argv); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Initialization failed!!!" << std::endl; + MPI_Abort(MPI_COMM_WORLD, 0); + } + + MPI_Comm comm = MPI_COMM_WORLD; + + /* + * Initialize general things, like magma, openmp, IO, ... + */ + mgmol_init(comm); + + /* + * read runtime parameters + */ + std::string input_filename(""); + std::string lrs_filename; + std::string constraints_filename(""); + + float total_spin = 0.; + bool with_spin = false; + + po::variables_map vm; + + // read from PE0 only + if (MPIdata::onpe0) + { + read_config(argc, argv, vm, input_filename, lrs_filename, + constraints_filename, total_spin, with_spin); + } + + MGmol_MPI::setup(comm, std::cout, with_spin); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + MPI_Comm global_comm = mmpi.commGlobal(); + + /* + * Setup control struct with run time parameters + */ + Control::setup(global_comm, with_spin, total_spin); + Control& ct = *(Control::instance()); + + ct.setOptions(vm); + + int ret = ct.checkOptions(); + if (ret < 0) return ret; + + mmpi.bcastGlobal(input_filename); + mmpi.bcastGlobal(lrs_filename); + + // Enter main scope + { + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "Construct MGmol object..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + + MGmolInterface* mgmol; + if (ct.isLocMode()) + 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); + + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "MGmol setup..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + mgmol->setup(); + + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "Setup done..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + + // here we just use the atomic positions read in and used + // to initialize MGmol + std::vector positions; + mgmol->getAtomicPositions(positions); + std::vector anumbers; + mgmol->getAtomicNumbers(anumbers); + if (MPIdata::onpe0) + { + std::cout << "Positions:" << std::endl; + std::vector::iterator ita = anumbers.begin(); + for (std::vector::iterator it = positions.begin(); + it != positions.end(); it += 3) + { + std::cout << *ita; + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + ita++; + } + } + + // compute energy and forces using all MPI tasks + // expect positions to be replicated on all MPI tasks + std::vector forces; + mgmol->evaluateEnergyAndForces(positions, anumbers, forces); + + // print out results + if (MPIdata::onpe0) + { + std::cout << "Forces:" << std::endl; + for (std::vector::iterator it = forces.begin(); + it != forces.end(); it += 3) + { + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + } + } + + delete mgmol; + + } // close main scope + + mgmol_finalize(); + + mpirc = MPI_Finalize(); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Finalize failed!!!" << std::endl; + } + + time_t tt; + time(&tt); + if (onpe0) std::cout << " Run ended at " << ctime(&tt) << std::endl; + + return 0; +} diff --git a/src/Ions.cc b/src/Ions.cc index 53ee0cf4..940e12e8 100644 --- a/src/Ions.cc +++ b/src/Ions.cc @@ -10,6 +10,7 @@ #include "Ions.h" #include "Control.h" #include "HDFrestart.h" +#include "MGmol_MPI.h" #include "MGmol_blas1.h" #include "MPIdata.h" #include "Mesh.h" @@ -169,12 +170,9 @@ void Ions::computeMaxNumProjs() << " initialized" << std::endl; #endif - std::vector::iterator ion = list_ions_.begin(); - while (ion != list_ions_.end()) + for (auto& ion : list_ions_) { - max_num_proj_ = std::max(max_num_proj_, (*ion)->nProjectors()); - - ion++; + max_num_proj_ = std::max(max_num_proj_, ion->nProjectors()); } MGmol_MPI& mmpi(*(MGmol_MPI::instance())); @@ -210,11 +208,9 @@ void Ions::setup() if (ct.verbose > 0) printWithTimeStamp("Ions::setup()... individual ions...", std::cout); - std::vector::iterator ion = list_ions_.begin(); - while (ion != list_ions_.end()) + for (auto& ion : list_ions_) { - (*ion)->setup(); - ion++; + ion->setup(); } setMapVL(); @@ -243,7 +239,7 @@ Ions::~Ions() void Ions::setupListOverlappingIons() { Control& ct = *(Control::instance()); - if (ct.verbose > 0) + if (ct.verbose > 1) printWithTimeStamp("Ions::setupListOverlappingIons()...", std::cout); overlappingNL_ions_.clear(); @@ -269,7 +265,7 @@ void Ions::setupListOverlappingIons() void Ions::setupInteractingIons() { Control& ct = *(Control::instance()); - if (ct.verbose > 0) + if (ct.verbose > 1) printWithTimeStamp("Ions::setupInteractingIons()...", std::cout); ions_setupInteractingIons_tm.start(); @@ -1569,7 +1565,7 @@ Ion* Ions::findLocalIon(const int index) const return nullptr; } -void Ions::setPositions(const std::vector& tau) +void Ions::setLocalPositions(const std::vector& tau) { assert(tau.size() == 3 * local_ions_.size()); @@ -1586,69 +1582,6 @@ void Ions::setPositions(const std::vector& tau) setup_ = false; } -void Ions::get_positions(std::vector>& rr) const -{ - assert(rr.size() == local_ions_.size()); - if (local_ions_.empty()) return; - std::vector tau(3); - int i = 0; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) - { - tau[0] = (*ion)->position(0); - tau[1] = (*ion)->position(1); - tau[2] = (*ion)->position(2); - rr[i] = tau; - ion++; - i++; - } -} -void Ions::set_positions(const std::vector>& rr) -{ - assert(rr.size() == local_ions_.size()); - - if (local_ions_.empty()) return; - std::vector::iterator ion = local_ions_.begin(); - int i = 0; - while (ion != local_ions_.end()) - { - assert(rr[i].size() == 3); - (*ion)->setPosition(rr[i][0], rr[i][1], rr[i][2]); - ion++; - i++; - } -} -void Ions::get_forces(std::vector>& ff) const -{ - assert(ff.size() == local_ions_.size()); - - if (local_ions_.empty()) return; - int i = 0; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) - { - ff[i][0] = (*ion)->force(0); - ff[i][1] = (*ion)->force(1); - ff[i][2] = (*ion)->force(2); - ion++; - i++; - } -} -void Ions::set_forces(const std::vector>& ff) -{ - assert(ff.size() == local_ions_.size()); - - if (local_ions_.empty()) return; - std::vector::iterator ion = local_ions_.begin(); - int i = 0; - while (ion != local_ions_.end()) - { - (*ion)->setForce(ff[i][0], ff[i][1], ff[i][2]); - ion++; - i++; - } -} - int Ions::readAtoms(const std::string& filename, const bool cell_relative) { @@ -1785,6 +1718,7 @@ int Ions::readAtomsFromXYZ( } ++count; } + delete tfile; } mmpi.bcastGlobal(&count, 1); @@ -1793,6 +1727,17 @@ int Ions::readAtomsFromXYZ( mmpi.bcastGlobal(&crds[0], 3 * natoms); mmpi.bcastGlobal(&spec[0], natoms); + return setAtoms(crds, spec); +} + +int Ions::setAtoms( + const std::vector& crds, const std::vector& spec) +{ + MGmol_MPI& mmpi(*(MGmol_MPI::instance())); + Control& ct(*(Control::instance())); + + const int natoms = crds.size() / 3; + double velocity[3] = { 0., 0., 0. }; bool locked = false; for (int ia = 0; ia < natoms; ++ia) @@ -1820,7 +1765,7 @@ int Ions::readAtomsFromXYZ( } if (spname.compare("") == 0) { - (*MPIdata::serr) << "Ions::readAtomsFromXYZ() --- ERROR: unknown " + (*MPIdata::serr) << "Ions::setAtoms() --- ERROR: unknown " "species for atomic number " << spec[ia] << std::endl; return -1; @@ -1841,7 +1786,7 @@ int Ions::readAtomsFromXYZ( // Populate list_ions_ list // std::cout<<"crds: "<set_here(true); + new_ion->set_here(true); local_ions_.push_back(new_ion); } else - (new_ion)->set_here(false); + new_ion->set_here(false); } else { @@ -2193,7 +2134,7 @@ void Ions::setVelocities(const std::vector& tau0, } } -void Ions::getPositions(std::vector& tau) const +void Ions::getLocalPositions(std::vector& tau) const { assert(tau.size() == 3 * local_ions_.size()); @@ -2207,6 +2148,42 @@ void Ions::getPositions(std::vector& tau) const } } +void Ions::getPositions(std::vector& tau) +{ + std::vector tau_local(3 * local_ions_.size()); + + getLocalPositions(tau_local); + + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + mmpi.allGatherV(tau_local, tau); +} + +void Ions::getAtomicNumbers(std::vector& atnumbers) +{ + std::vector local_atnumbers; + + for (auto& ion : local_ions_) + { + local_atnumbers.push_back(ion->atomic_number()); + } + + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + mmpi.allGatherV(local_atnumbers, atnumbers); +} + +void Ions::getForces(std::vector& forces) +{ + std::vector forces_local(3 * local_ions_.size()); + + getLocalForces(forces_local); + + int n = getNumIons(); + forces.resize(3 * n); + + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + mmpi.allGatherV(forces_local, forces); +} + void Ions::setTau0() { assert(tau0_.size() == 3 * local_ions_.size()); @@ -2236,6 +2213,20 @@ void Ions::setPositionsToTau0() } } +void Ions::setPositions( + const std::vector& tau, const std::vector& anumbers) +{ + assert(tau.size() == anumbers.size() * 3); + + // clear previous data + clearLists(); + + num_ions_ = setAtoms(tau, anumbers); + + // setup required after updating local ions positions + setup(); +} + void Ions::setVelocitiesToVel() { assert(velocity_.size() == 3 * local_ions_.size()); @@ -2249,7 +2240,7 @@ void Ions::setVelocitiesToVel() } } -void Ions::getForces(std::vector& tau) const +void Ions::getLocalForces(std::vector& tau) const { assert(tau.size() == 3 * local_ions_.size()); @@ -2380,12 +2371,10 @@ double Ions::computeMaxNLprojRadius() const { double radius = 0.; - std::vector::const_iterator iion = local_ions_.begin(); - while (iion != local_ions_.end()) + for (auto& iion : local_ions_) { - double r = (*iion)->radiusNLproj(); + double r = iion->radiusNLproj(); radius = r > radius ? r : radius; - iion++; } MGmol_MPI& mmpi(*(MGmol_MPI::instance())); @@ -3061,7 +3050,7 @@ void Ions::updateForcesInteractingIons() MGmol_MPI& mmpi(*(MGmol_MPI::instance())); // get computed forces into fion_ - getForces(fion_); + getLocalForces(fion_); // initialize with local names and forces DistributedIonicData forces_data(local_names_, fion_); @@ -3136,6 +3125,18 @@ void Ions::updateTaupInteractingIons() } } +void Ions::clearLists() +{ + local_ions_.clear(); + std::vector::iterator ion = list_ions_.begin(); + while (ion != list_ions_.end()) + { + delete *ion; + ion++; + } + list_ions_.clear(); +} + // update list of local ions void Ions::updateListIons() { @@ -3160,15 +3161,7 @@ void Ions::updateListIons() // Note: this is based on data from MD std::vectors // First cleanup list_ions_ - local_ions_.clear(); - // delete current ions from list - std::vector::iterator ion = list_ions_.begin(); - while (ion != list_ions_.end()) - { - delete *ion; - ion++; - } - list_ions_.clear(); + clearLists(); // Update list starting with local ion data. // This enables overlapping data accumulation with communication. diff --git a/src/Ions.h b/src/Ions.h index 2351d678..772c51d1 100644 --- a/src/Ions.h +++ b/src/Ions.h @@ -16,7 +16,6 @@ #include #include -#include "DataDistribution.h" #include "DistributedIonicData.h" #include "Ion.h" #include "hdf5.h" @@ -37,7 +36,12 @@ class Ions const std::vector& species_; std::vector list_ions_; - std::vector local_ions_; // centered in local sub-domain + + /* + * ions located in local sub-domain + */ + std::vector local_ions_; + std::vector interacting_ions_; // for ion-ion interactions std::vector overlappingNL_ions_; // with projectors overlapping local sub-domain @@ -62,7 +66,6 @@ class Ions void readRestartPositions(HDFrestart& h5_file); int read1atom(std::ifstream* tfile, const bool cell_relative); - // void associate2PE(); void setupInteractingIons(); void setupListOverlappingIons(); void setMapVL(); @@ -158,6 +161,7 @@ class Ions void gatherForces( std::vector& forces, const int root, const MPI_Comm comm) const; bool hasLockedAtoms() const; + void clearLists(); public: Ions(const double lat[3], const std::vector& sp); @@ -247,11 +251,7 @@ class Ions // check if ion is in list of local ions bool isLocal(const std::string& ion_name) const; - void setPositions(const std::vector& tau); - void get_positions(std::vector>& r) const; - void set_positions(const std::vector>& r); - void get_forces(std::vector>& f) const; - void set_forces(const std::vector>& f); + void setLocalPositions(const std::vector& tau); void lockAtom(const std::string& name); @@ -260,6 +260,8 @@ class Ions int readAtoms(const std::string& filename, const bool cell_relative); int readAtoms(std::ifstream* tfile, const bool cell_relative); void initFromRestartFile(HDFrestart& h5_file); + int setAtoms( + const std::vector& crds, const std::vector& spec); int getNValenceElectrons() const; void syncForces(); @@ -268,9 +270,15 @@ class Ions void setTau0(); void setPositionsToTau0(); void setVelocitiesToVel(); + void setPositions( + const std::vector& tau, const std::vector& anumbers); + + void getLocalPositions(std::vector& tau) const; + void getPositions(std::vector& tau); + void getAtomicNumbers(std::vector& atnumbers); - void getPositions(std::vector& tau) const; - void getForces(std::vector& tau) const; + void getForces(std::vector& forces); + void getLocalForces(std::vector& tau) const; void syncData(const std::vector& sp); // void syncNames(const int nions, std::vector& local_names, // std::vector& names); diff --git a/src/MGmol.cc b/src/MGmol.cc index 85acc78d..e1f65d21 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -1405,6 +1405,39 @@ void MGmol::addResidualSpreadPenalty( spread_penalty_->addResidual(phi, res); } +template +void MGmol::getAtomicPositions(std::vector& tau) +{ + ions_->getPositions(tau); +} + +template +void MGmol::getAtomicNumbers(std::vector& an) +{ + ions_->getAtomicNumbers(an); +} + +template +double MGmol::evaluateEnergyAndForces( + const std::vector& tau, std::vector& atnumbers, + std::vector& forces) +{ + assert(tau.size() == 3 * atnumbers.size()); + + Control& ct = *(Control::instance()); + + ions_->setPositions(tau, atnumbers); + + double eks = 0.; + quench(current_orbitals_, *ions_, ct.max_electronic_steps, 20, eks); + + force(*current_orbitals_, *ions_); + + ions_->getForces(forces); + + return eks; +} + template class MGmol; template class MGmol; template int MGmol::initial(); diff --git a/src/MGmol.h b/src/MGmol.h index ceec7390..cb5c54f3 100644 --- a/src/MGmol.h +++ b/src/MGmol.h @@ -179,6 +179,17 @@ class MGmol : public MGmolInterface ~MGmol() override; void run() override; + + double evaluateEnergyAndForces(const std::vector& tau, + std::vector& atnumbers, std::vector& forces); + + /* + * get internal atomic positions + */ + void getAtomicPositions(std::vector& tau); + + void getAtomicNumbers(std::vector& an); + void initNuc(Ions& ions); void initKBR(); @@ -276,10 +287,6 @@ class MGmol : public MGmolInterface double get_evnl(const Ions& ions); void sebprintPositions(); void sebprintForces(); - void get_positions(std::vector>& r); - void set_positions(std::vector>& r); - void get_forces(std::vector>& f); - void set_forces(std::vector>& f); int nions() { return ions_->getNumIons(); } double getTotalEnergy(); void cleanup(); diff --git a/src/MGmolInterface.h b/src/MGmolInterface.h index 9932dd9f..046459d6 100644 --- a/src/MGmolInterface.h +++ b/src/MGmolInterface.h @@ -11,6 +11,7 @@ #define MGMOLINTERFACE_H #include +#include class MGmolInterface { @@ -24,6 +25,11 @@ class MGmolInterface virtual int setupConstraintsFromInput(const std::string input_file) = 0; virtual void setup() = 0; virtual void run() = 0; + virtual double evaluateEnergyAndForces(const std::vector& tau, + std::vector& atnumbers, std::vector& forces) + = 0; + virtual void getAtomicPositions(std::vector& tau) = 0; + virtual void getAtomicNumbers(std::vector& an) = 0; }; #endif diff --git a/src/md.cc b/src/md.cc index 41dcd39b..5f690e82 100644 --- a/src/md.cc +++ b/src/md.cc @@ -344,8 +344,9 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) constraints_->enforceConstraints(20); stepper->updateTau(); - ions.setPositions(tau0); + ions.setLocalPositions(tau0); + // setup required after updating local ions positions ions.setup(); } @@ -489,7 +490,7 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) force(**orbitals, ions); // set fion - ions.getForces(fion); + ions.getLocalForces(fion); // constraints need to be added again and setup as atoms // may have moved and local atoms are not the same anymore @@ -717,12 +718,15 @@ OrbitalsType* MGmol::loadOrbitalFromRestartFile( while the extrapolated function is set as the current orbitals. This is how the restart file is saved via dumprestartFile. - For now, we just enforce to not use the restart files with extrapolation. + For now, we just enforce to not use the restart files with + extrapolation. */ if (flag_extrapolated_data) { if (onpe0) - (*MPIdata::serr) << "loadRestartFile: does not support restart files with extrapolation." << std::endl; + (*MPIdata::serr) << "loadRestartFile: does not support restart " + "files with extrapolation." + << std::endl; global_exit(0); } From 13b86427defbd1412c219cc4092f5eda571865c7 Mon Sep 17 00:00:00 2001 From: "Kevin\" Seung Whan Chung" Date: Thu, 13 Jun 2024 18:34:29 -0700 Subject: [PATCH 07/17] loadOrbitalsFromRestartFile -> loadRestartFile (#253) --- src/MGmol.h | 6 ++++- src/md.cc | 67 +++++++---------------------------------------------- 2 files changed, 13 insertions(+), 60 deletions(-) diff --git a/src/MGmol.h b/src/MGmol.h index cb5c54f3..12096a5b 100644 --- a/src/MGmol.h +++ b/src/MGmol.h @@ -178,6 +178,10 @@ class MGmol : public MGmolInterface ~MGmol() override; + /* access functions */ + OrbitalsType* getOrbitals() { return current_orbitals_; } + std::shared_ptr> getHamiltonian() { return hamiltonian_; } + void run() override; double evaluateEnergyAndForces(const std::vector& tau, @@ -313,7 +317,7 @@ class MGmol : public MGmolInterface forces_->force(orbitals, ions); } - OrbitalsType* loadOrbitalFromRestartFile(const std::string filename); + void loadRestartFile(const std::string filename); }; // Instantiate static variables here to avoid clang warnings template diff --git a/src/md.cc b/src/md.cc index 5f690e82..5e016941 100644 --- a/src/md.cc +++ b/src/md.cc @@ -669,7 +669,7 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) } template -OrbitalsType* MGmol::loadOrbitalFromRestartFile( +void MGmol::loadRestartFile( const std::string filename) { MGmol_MPI& mmpi(*(MGmol_MPI::instance())); @@ -683,56 +683,15 @@ OrbitalsType* MGmol::loadOrbitalFromRestartFile( || (ct.AtomsDynamic() == AtomsDynamicType::Quench)); HDFrestart h5file(filename, myPEenv, ct.out_restart_file_type); - int ierr; - OrbitalsType* restart_orbitals = nullptr; - - /* This corresponds to MGmol::initial */ - { - // Copy from current orbital, instead of constructing brand-new one - restart_orbitals - = new OrbitalsType("ForLoading", *current_orbitals_, false); - - /* This corresponds to MGmol::read_restart_data */ - { - ierr = restart_orbitals->read_func_hdf5(h5file); - } // read_restart_data - } // initial() - - /* This corresponds to MGmol::md */ + int ierr = read_restart_data(h5file, *rho_, *current_orbitals_); + if (ierr < 0) { - int flag_extrapolated_data = 0; if (onpe0) - { - flag_extrapolated_data - = h5file.dset_exists("ExtrapolatedFunction0000"); - if (flag_extrapolated_data == 0) - flag_extrapolated_data - = h5file.dset_exists("ExtrapolatedFunction0"); - } - MPI_Bcast(&flag_extrapolated_data, 1, MPI_INT, 0, comm_); - - /* - If extrapolated function exists, - then function is set as previous orbitals, - while the extrapolated function is set as the current orbitals. - This is how the restart file is saved via dumprestartFile. - - For now, we just enforce to not use the restart files with - extrapolation. - */ - if (flag_extrapolated_data) - { - if (onpe0) - (*MPIdata::serr) << "loadRestartFile: does not support restart " - "files with extrapolation." - << std::endl; + (*MPIdata::serr) << "loadRestartFile: failed to read the restart file." << std::endl; - global_exit(0); - } - - /* main workflow delete h5f_file_ here, meaning the loading is over. */ - } // md() + global_exit(0); + } ierr = h5file.close(); mmpi.allreduce(&ierr, 1, MPI_MIN); @@ -741,20 +700,10 @@ OrbitalsType* MGmol::loadOrbitalFromRestartFile( if (onpe0) (*MPIdata::serr) << "loadRestartFile: cannot close file..." << std::endl; - return nullptr; + return; } - /* - In returning restart_orbitals, - we hope that the wavefunctions in restart_orbitals are all set. - At least the following functions should return proper data loaded from - the file: - - restart_orbitals->getLocNumpt() - restart_orbitals->chromatic_number() - restart_orbitals->getPsi(idx) (for int idx) - */ - return restart_orbitals; + return; } template class MGmol; From 3eaf5a6a021374474569950c0c2eedcefa796b38 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Sat, 20 Jul 2024 10:41:40 -0400 Subject: [PATCH 08/17] Add SG15 PBE potential for N (#258) --- potentials/pseudo.N_ONCV_PBE_SG15 | 1863 +++++++++++++++++++++++++++++ 1 file changed, 1863 insertions(+) create mode 100644 potentials/pseudo.N_ONCV_PBE_SG15 diff --git a/potentials/pseudo.N_ONCV_PBE_SG15 b/potentials/pseudo.N_ONCV_PBE_SG15 new file mode 100644 index 00000000..958c0c8b --- /dev/null +++ b/potentials/pseudo.N_ONCV_PBE_SG15 @@ -0,0 +1,1863 @@ +# This pseudopotential file has been produced using the code +# ONCVPSP (Optimized Norm-Conservinng Vanderbilt PSeudopotential) +# scalar-relativistic version 2.1.1, 03/26/2014 by D. R. Hamann +# as distributed by Mat-Sim Research at www.mat-simresearch.com. +# Reference: +# D. R. Hamann, Phys. Rev. B 88, 085117 (2013) +# DOI:https://doi.org/10.1103/PhysRevB.88.085117 +# +# Pseudopotential generation was done by M. Schlipf and F. Gygi using +# the input parameters obtained from their optimization algorithm. +# Reference: +# M. Schlipf and F. Gygi, Computer Physics Communications (2015) +# DOI: 10.1016/j.cpc.2015.05.011 +# http://dx.doi.org/10.1016/j.cpc.2015.05.011 +# http://www.quantum-simulation.org +# License: +# This work is licensed under the Creative Commons Attribution-ShareAlike +# 4.0 International License. To view a copy of this license, visit +# http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to +# Creative Commons, PO Box 1866, Mountain View, CA 94042, USA. +# +# Format: +# The format of this file has been adapted to match MGmol pseudopotential +# parser using the script convertONCVPSP_SG15.py +# +N_ONCV_PBE-1 +# +color +#radii of balls and covalent bonds +-1. -1. +# Nlcc flag +0 +# Atomic number +7 +# Atomic mass +14.00699997 +# Number of valence electrons +5 +#Gaussian core charge parameter rc +1.0 +# Number of potentials +3 +# l-value for state which is local, then type of potential format +2 3 +# Local potential radius +3.2 +# Non-local potential radius +3.2 +# number of points in radial grid +602 +# VANDERBILT-KLEINMAN-BYLANDER PROJECTORs +# l, nproj +0 2 0.7070940539E+01 0.9781080649E+00 +1 2 -0.4781783452E+01 -0.1328201236E+01 +# l= 0 +0.0 -8.3833937254 3.4232195615 +0.01 -8.3799978758 3.4202071705 +0.02 -8.3697860653 3.4111530231 +0.03 -8.3527825341 3.3960738696 +0.04 -8.3290115228 3.3749864603 +0.05 -8.2985069887 3.3479144681 +0.06 -8.261312626 3.3148887067 +0.07 -8.2174818893 3.2759474053 +0.08 -8.1670780196 3.231136536 +0.09 -8.1101740718 3.1805101875 +0.1 -8.0468529415 3.1241309807 +0.11 -7.9772073885 3.0620705184 +0.12 -7.9013400567 2.9944098623 +0.13 -7.8193634857 2.9212400287 +0.14 -7.7314001139 2.8426624957 +0.15 -7.6375822685 2.7587897114 +0.16 -7.5380521411 2.6697455952 +0.17 -7.4329617463 2.5756660223 +0.18 -7.3224728588 2.4766992817 +0.19 -7.2067569298 2.3730065 +0.2 -7.0859949741 2.264762015 +0.21 -6.9603774349 2.1521537041 +0.22 -6.8301040119 2.0353832393 +0.23 -6.6953834613 1.9146662776 +0.24 -6.5564333561 1.7902325648 +0.25 -6.4134798193 1.6623259662 +0.26 -6.2667572033 1.5312043795 +0.27 -6.1165077488 1.3971395792 +0.28 -5.962981185 1.2604169264 +0.29 -5.8064342997 1.1213349889 +0.3 -5.647130465 0.98020504493 +0.31 -5.4853391206 0.83735047259 +0.32 -5.321335212 0.69310601673 +0.33 -5.1553985987 0.54781696226 +0.34 -4.9878134137 0.4018381796 +0.35 -4.8188673807 0.2555330423 +0.36 -4.6488511141 0.10927228913 +0.37 -4.4780573629 -0.036567271748 +0.38 -4.3067802363 -0.18160413005 +0.39 -4.1353144013 -0.32545351173 +0.4 -3.9639542476 -0.46772893557 +0.41 -3.7929930387 -0.60804381665 +0.42 -3.6227220421 -0.74601313801 +0.43 -3.4534296479 -0.88125517478 +0.44 -3.2854004768 -1.0133932584 +0.45 -3.1189144859 -1.1420575692 +0.46 -2.9542460734 -1.2668869442 +0.47 -2.7916631913 -1.387530686 +0.48 -2.6314265018 -1.5036503015 +0.49 -2.47378844 -1.6149214048 +0.5 -2.3189924173 -1.7210353819 +0.51 -2.1672720468 -1.8217009988 +0.52 -2.0188503481 -1.9166460389 +0.53 -1.8739388537 -2.0056191124 +0.54 -1.732737225 -2.0883905138 +0.55 -1.5954322069 -2.1647542719 +0.56 -1.4621973571 -2.2345287551 +0.57 -1.333192259 -2.2975581981 +0.58 -1.2085621978 -2.353713346 +0.59 -1.0884375079 -2.4028926755 +0.6 -0.97293360827 -2.4450223075 +0.61 -0.86215002573 -2.4800577568 +0.62 -0.75617104038 -2.5079826487 +0.63 -0.65506469693 -2.5288104095 +0.64 -0.55888310603 -2.5425835395 +0.65 -0.46766273408 -2.5493728888 +0.66 -0.38142374876 -2.5492785866 +0.67 -0.30017057762 -2.5424287393 +0.68 -0.22389217436 -2.5289786283 +0.69 -0.15256244154 -2.5091096125 +0.7 -0.086139788103 -2.4830294406 +0.71 -0.024568342949 -2.4509697004 +0.72 0.032221809773 -2.413184916 +0.73 0.084314093363 -2.3699510915 +0.74 0.1318047363 -2.321563976 +0.75 0.17480239707 -2.2683377698 +0.76 0.2134272741 -2.2106030175 +0.77 0.24781028338 -2.1487046133 +0.78 0.27809232097 -2.0829999018 +0.79 0.30442342743 -2.0138566245 +0.8 0.32696191973 -1.9416508145 +0.81 0.34587349606 -1.8667646538 +0.82 0.36133031921 -1.7895843109 +0.83 0.37351008427 -1.7104977709 +0.84 0.38259507609 -1.6298926758 +0.85 0.38877122239 -1.5481541871 +0.86 0.39222714765 -1.4656628877 +0.87 0.39315323345 -1.3827927341 +0.88 0.39174069021 -1.2999090738 +0.89 0.3881806455 -1.21736674 +0.9 0.38266325363 -1.1355082343 +0.91 0.3753768312 -1.05466201 +0.92 0.36650702291 -0.97514086392 +0.93 0.35623600185 -0.8972404476 +0.94 0.34474170806 -0.82123790516 +0.95 0.33219712893 -0.74739064427 +0.96 0.31876980308 -0.67593583254 +0.97 0.30462077846 -0.60708800794 +0.98 0.28990423397 -0.54103907779 +0.99 0.27476701087 -0.47795782396 +1.0 0.25934807756 -0.41798919045 +1.01 0.24377810285 -0.3612539062 +1.02 0.22817911217 -0.30784832508 +1.03 0.21266470226 -0.2578456722 +1.04 0.19733797262 -0.21129138178 +1.05 0.18229298665 -0.1682077754 +1.06 0.16761409975 -0.12859294839 +1.07 0.15337584366 -0.092421151151 +1.08 0.1396438166 -0.059645387463 +1.09 0.12647159361 -0.030191037026 +1.1 0.11390353081 -0.0039630063463 +1.11 0.10197406871 0.019155281934 +1.12 0.090708799358 0.039298886262 +1.13 0.080123593815 0.056619471672 +1.14 0.070227673946 0.071275245192 +1.15 0.061024372085 0.083429790064 +1.16 0.052511962988 0.093250798533 +1.17 0.044682028129 0.10090978584 +1.18 0.037521888086 0.10657912289 +1.19 0.031015059574 0.11043138827 +1.2 0.025141213318 0.11263866157 +1.21 0.019876005667 0.11336980971 +1.22 0.015192955473 0.11279033889 +1.23 0.011063595893 0.11106306695 +1.24 0.007456655257 0.10834247812 +1.25 0.0043400745021 0.10477759314 +1.26 0.0016808247737 0.10051281117 +1.27 -0.00055510433714 0.095680812463 +1.28 -0.0024016580764 0.090406927903 +1.29 -0.0038929622779 0.084811273797 +1.3 -0.0050625477429 0.078998968474 +1.31 -0.0059432416054 0.073068149935 +1.32 -0.0065674921377 0.067110734764 +1.33 -0.0069657771649 0.061201995595 +1.34 -0.007167828516 0.055412646763 +1.35 -0.0072020326121 0.049803640184 +1.36 -0.0070943360586 0.044422493293 +1.37 -0.0068701508663 0.039313829408 +1.38 -0.0065520891543 0.034508390981 +1.39 -0.0061611134842 0.030030366818 +1.4 -0.005717227564 0.025899372872 +1.41 -0.005236852799 0.022121815245 +1.42 -0.0047362333134 0.018703666441 +1.43 -0.0042285949508 0.01564086524 +1.44 -0.0037256393417 0.012925597893 +1.45 -0.003238024767 0.010546723262 +1.46 -0.0027732187555 0.0084862043738 +1.47 -0.0023389767364 0.0067264437607 +1.48 -0.0019396582406 0.005243996278 +1.49 -0.0015795366936 0.0040158727896 +1.5 -0.0012608949043 0.0030165453042 +1.51 -0.00098475359848 0.0022205393899 +1.52 -0.00075153947202 0.0016016330403 +1.53 -0.00055956428127 0.0011349570223 +1.54 -0.000407575033 0.00079484813589 +1.55 -0.00029180979095 0.00055922039698 +1.56 -0.00020917873683 0.00040490419124 +1.57 -0.00015459824577 0.00031332706041 +1.58 -0.00012325344328 0.00026554043958 +1.59 -0.00011211107894 0.00025459512511 +1.6 -0.00011091815786 0.00025685343044 +1.61 -8.0065645976e-05 0.00018636643746 +1.62 -2.207959279e-05 5.1598480417e-05 +1.63 6.3416160474e-06 -1.4718026015e-05 +1.64 3.4369124819e-06 -7.9766051653e-06 +1.65 0.0 0.0 +1.66 0.0 0.0 +1.67 0.0 0.0 +1.68 0.0 0.0 +1.69 0.0 0.0 +1.7 0.0 0.0 +1.71 0.0 0.0 +1.72 0.0 0.0 +1.73 0.0 0.0 +1.74 0.0 0.0 +1.75 0.0 0.0 +1.76 0.0 0.0 +1.77 0.0 0.0 +1.78 0.0 0.0 +1.79 0.0 0.0 +1.8 0.0 0.0 +1.81 0.0 0.0 +1.82 0.0 0.0 +1.83 0.0 0.0 +1.84 0.0 0.0 +1.85 0.0 0.0 +1.86 0.0 0.0 +1.87 0.0 0.0 +1.88 0.0 0.0 +1.89 0.0 0.0 +1.9 0.0 0.0 +1.91 0.0 0.0 +1.92 0.0 0.0 +1.93 0.0 0.0 +1.94 0.0 0.0 +1.95 0.0 0.0 +1.96 0.0 0.0 +1.97 0.0 0.0 +1.98 0.0 0.0 +1.99 0.0 0.0 +2.0 0.0 0.0 +2.01 0.0 0.0 +2.02 0.0 0.0 +2.03 0.0 0.0 +2.04 0.0 0.0 +2.05 0.0 0.0 +2.06 0.0 0.0 +2.07 0.0 0.0 +2.08 0.0 0.0 +2.09 0.0 0.0 +2.1 0.0 0.0 +2.11 0.0 0.0 +2.12 0.0 0.0 +2.13 0.0 0.0 +2.14 0.0 0.0 +2.15 0.0 0.0 +2.16 0.0 0.0 +2.17 0.0 0.0 +2.18 0.0 0.0 +2.19 0.0 0.0 +2.2 0.0 0.0 +2.21 0.0 0.0 +2.22 0.0 0.0 +2.23 0.0 0.0 +2.24 0.0 0.0 +2.25 0.0 0.0 +2.26 0.0 0.0 +2.27 0.0 0.0 +2.28 0.0 0.0 +2.29 0.0 0.0 +2.3 0.0 0.0 +2.31 0.0 0.0 +2.32 0.0 0.0 +2.33 0.0 0.0 +2.34 0.0 0.0 +2.35 0.0 0.0 +2.36 0.0 0.0 +2.37 0.0 0.0 +2.38 0.0 0.0 +2.39 0.0 0.0 +2.4 0.0 0.0 +2.41 0.0 0.0 +2.42 0.0 0.0 +2.43 0.0 0.0 +2.44 0.0 0.0 +2.45 0.0 0.0 +2.46 0.0 0.0 +2.47 0.0 0.0 +2.48 0.0 0.0 +2.49 0.0 0.0 +2.5 0.0 0.0 +2.51 0.0 0.0 +2.52 0.0 0.0 +2.53 0.0 0.0 +2.54 0.0 0.0 +2.55 0.0 0.0 +2.56 0.0 0.0 +2.57 0.0 0.0 +2.58 0.0 0.0 +2.59 0.0 0.0 +2.6 0.0 0.0 +2.61 0.0 0.0 +2.62 0.0 0.0 +2.63 0.0 0.0 +2.64 0.0 0.0 +2.65 0.0 0.0 +2.66 0.0 0.0 +2.67 0.0 0.0 +2.68 0.0 0.0 +2.69 0.0 0.0 +2.7 0.0 0.0 +2.71 0.0 0.0 +2.72 0.0 0.0 +2.73 0.0 0.0 +2.74 0.0 0.0 +2.75 0.0 0.0 +2.76 0.0 0.0 +2.77 0.0 0.0 +2.78 0.0 0.0 +2.79 0.0 0.0 +2.8 0.0 0.0 +2.81 0.0 0.0 +2.82 0.0 0.0 +2.83 0.0 0.0 +2.84 0.0 0.0 +2.85 0.0 0.0 +2.86 0.0 0.0 +2.87 0.0 0.0 +2.88 0.0 0.0 +2.89 0.0 0.0 +2.9 0.0 0.0 +2.91 0.0 0.0 +2.92 0.0 0.0 +2.93 0.0 0.0 +2.94 0.0 0.0 +2.95 0.0 0.0 +2.96 0.0 0.0 +2.97 0.0 0.0 +2.98 0.0 0.0 +2.99 0.0 0.0 +3.0 0.0 0.0 +3.01 0.0 0.0 +3.02 0.0 0.0 +3.03 0.0 0.0 +3.04 0.0 0.0 +3.05 0.0 0.0 +3.06 0.0 0.0 +3.07 0.0 0.0 +3.08 0.0 0.0 +3.09 0.0 0.0 +3.1 0.0 0.0 +3.11 0.0 0.0 +3.12 0.0 0.0 +3.13 0.0 0.0 +3.14 0.0 0.0 +3.15 0.0 0.0 +3.16 0.0 0.0 +3.17 0.0 0.0 +3.18 0.0 0.0 +3.19 0.0 0.0 +3.2 0.0 0.0 +3.21 0.0 0.0 +3.22 0.0 0.0 +3.23 0.0 0.0 +3.24 0.0 0.0 +3.25 0.0 0.0 +3.26 0.0 0.0 +3.27 0.0 0.0 +3.28 0.0 0.0 +3.29 0.0 0.0 +3.3 0.0 0.0 +3.31 0.0 0.0 +3.32 0.0 0.0 +3.33 0.0 0.0 +3.34 0.0 0.0 +3.35 0.0 0.0 +3.36 0.0 0.0 +3.37 0.0 0.0 +3.38 0.0 0.0 +3.39 0.0 0.0 +3.4 0.0 0.0 +3.41 0.0 0.0 +3.42 0.0 0.0 +3.43 0.0 0.0 +3.44 0.0 0.0 +3.45 0.0 0.0 +3.46 0.0 0.0 +3.47 0.0 0.0 +3.48 0.0 0.0 +3.49 0.0 0.0 +3.5 0.0 0.0 +3.51 0.0 0.0 +3.52 0.0 0.0 +3.53 0.0 0.0 +3.54 0.0 0.0 +3.55 0.0 0.0 +3.56 0.0 0.0 +3.57 0.0 0.0 +3.58 0.0 0.0 +3.59 0.0 0.0 +3.6 0.0 0.0 +3.61 0.0 0.0 +3.62 0.0 0.0 +3.63 0.0 0.0 +3.64 0.0 0.0 +3.65 0.0 0.0 +3.66 0.0 0.0 +3.67 0.0 0.0 +3.68 0.0 0.0 +3.69 0.0 0.0 +3.7 0.0 0.0 +3.71 0.0 0.0 +3.72 0.0 0.0 +3.73 0.0 0.0 +3.74 0.0 0.0 +3.75 0.0 0.0 +3.76 0.0 0.0 +3.77 0.0 0.0 +3.78 0.0 0.0 +3.79 0.0 0.0 +3.8 0.0 0.0 +3.81 0.0 0.0 +3.82 0.0 0.0 +3.83 0.0 0.0 +3.84 0.0 0.0 +3.85 0.0 0.0 +3.86 0.0 0.0 +3.87 0.0 0.0 +3.88 0.0 0.0 +3.89 0.0 0.0 +3.9 0.0 0.0 +3.91 0.0 0.0 +3.92 0.0 0.0 +3.93 0.0 0.0 +3.94 0.0 0.0 +3.95 0.0 0.0 +3.96 0.0 0.0 +3.97 0.0 0.0 +3.98 0.0 0.0 +3.99 0.0 0.0 +4.0 0.0 0.0 +4.01 0.0 0.0 +4.02 0.0 0.0 +4.03 0.0 0.0 +4.04 0.0 0.0 +4.05 0.0 0.0 +4.06 0.0 0.0 +4.07 0.0 0.0 +4.08 0.0 0.0 +4.09 0.0 0.0 +4.1 0.0 0.0 +4.11 0.0 0.0 +4.12 0.0 0.0 +4.13 0.0 0.0 +4.14 0.0 0.0 +4.15 0.0 0.0 +4.16 0.0 0.0 +4.17 0.0 0.0 +4.18 0.0 0.0 +4.19 0.0 0.0 +4.2 0.0 0.0 +4.21 0.0 0.0 +4.22 0.0 0.0 +4.23 0.0 0.0 +4.24 0.0 0.0 +4.25 0.0 0.0 +4.26 0.0 0.0 +4.27 0.0 0.0 +4.28 0.0 0.0 +4.29 0.0 0.0 +4.3 0.0 0.0 +4.31 0.0 0.0 +4.32 0.0 0.0 +4.33 0.0 0.0 +4.34 0.0 0.0 +4.35 0.0 0.0 +4.36 0.0 0.0 +4.37 0.0 0.0 +4.38 0.0 0.0 +4.39 0.0 0.0 +4.4 0.0 0.0 +4.41 0.0 0.0 +4.42 0.0 0.0 +4.43 0.0 0.0 +4.44 0.0 0.0 +4.45 0.0 0.0 +4.46 0.0 0.0 +4.47 0.0 0.0 +4.48 0.0 0.0 +4.49 0.0 0.0 +4.5 0.0 0.0 +4.51 0.0 0.0 +4.52 0.0 0.0 +4.53 0.0 0.0 +4.54 0.0 0.0 +4.55 0.0 0.0 +4.56 0.0 0.0 +4.57 0.0 0.0 +4.58 0.0 0.0 +4.59 0.0 0.0 +4.6 0.0 0.0 +4.61 0.0 0.0 +4.62 0.0 0.0 +4.63 0.0 0.0 +4.64 0.0 0.0 +4.65 0.0 0.0 +4.66 0.0 0.0 +4.67 0.0 0.0 +4.68 0.0 0.0 +4.69 0.0 0.0 +4.7 0.0 0.0 +4.71 0.0 0.0 +4.72 0.0 0.0 +4.73 0.0 0.0 +4.74 0.0 0.0 +4.75 0.0 0.0 +4.76 0.0 0.0 +4.77 0.0 0.0 +4.78 0.0 0.0 +4.79 0.0 0.0 +4.8 0.0 0.0 +4.81 0.0 0.0 +4.82 0.0 0.0 +4.83 0.0 0.0 +4.84 0.0 0.0 +4.85 0.0 0.0 +4.86 0.0 0.0 +4.87 0.0 0.0 +4.88 0.0 0.0 +4.89 0.0 0.0 +4.9 0.0 0.0 +4.91 0.0 0.0 +4.92 0.0 0.0 +4.93 0.0 0.0 +4.94 0.0 0.0 +4.95 0.0 0.0 +4.96 0.0 0.0 +4.97 0.0 0.0 +4.98 0.0 0.0 +4.99 0.0 0.0 +5.0 0.0 0.0 +5.01 0.0 0.0 +5.02 0.0 0.0 +5.03 0.0 0.0 +5.04 0.0 0.0 +5.05 0.0 0.0 +5.06 0.0 0.0 +5.07 0.0 0.0 +5.08 0.0 0.0 +5.09 0.0 0.0 +5.1 0.0 0.0 +5.11 0.0 0.0 +5.12 0.0 0.0 +5.13 0.0 0.0 +5.14 0.0 0.0 +5.15 0.0 0.0 +5.16 0.0 0.0 +5.17 0.0 0.0 +5.18 0.0 0.0 +5.19 0.0 0.0 +5.2 0.0 0.0 +5.21 0.0 0.0 +5.22 0.0 0.0 +5.23 0.0 0.0 +5.24 0.0 0.0 +5.25 0.0 0.0 +5.26 0.0 0.0 +5.27 0.0 0.0 +5.28 0.0 0.0 +5.29 0.0 0.0 +5.3 0.0 0.0 +5.31 0.0 0.0 +5.32 0.0 0.0 +5.33 0.0 0.0 +5.34 0.0 0.0 +5.35 0.0 0.0 +5.36 0.0 0.0 +5.37 0.0 0.0 +5.38 0.0 0.0 +5.39 0.0 0.0 +5.4 0.0 0.0 +5.41 0.0 0.0 +5.42 0.0 0.0 +5.43 0.0 0.0 +5.44 0.0 0.0 +5.45 0.0 0.0 +5.46 0.0 0.0 +5.47 0.0 0.0 +5.48 0.0 0.0 +5.49 0.0 0.0 +5.5 0.0 0.0 +5.51 0.0 0.0 +5.52 0.0 0.0 +5.53 0.0 0.0 +5.54 0.0 0.0 +5.55 0.0 0.0 +5.56 0.0 0.0 +5.57 0.0 0.0 +5.58 0.0 0.0 +5.59 0.0 0.0 +5.6 0.0 0.0 +5.61 0.0 0.0 +5.62 0.0 0.0 +5.63 0.0 0.0 +5.64 0.0 0.0 +5.65 0.0 0.0 +5.66 0.0 0.0 +5.67 0.0 0.0 +5.68 0.0 0.0 +5.69 0.0 0.0 +5.7 0.0 0.0 +5.71 0.0 0.0 +5.72 0.0 0.0 +5.73 0.0 0.0 +5.74 0.0 0.0 +5.75 0.0 0.0 +5.76 0.0 0.0 +5.77 0.0 0.0 +5.78 0.0 0.0 +5.79 0.0 0.0 +5.8 0.0 0.0 +5.81 0.0 0.0 +5.82 0.0 0.0 +5.83 0.0 0.0 +5.84 0.0 0.0 +5.85 0.0 0.0 +5.86 0.0 0.0 +5.87 0.0 0.0 +5.88 0.0 0.0 +5.89 0.0 0.0 +5.9 0.0 0.0 +5.91 0.0 0.0 +5.92 0.0 0.0 +5.93 0.0 0.0 +5.94 0.0 0.0 +5.95 0.0 0.0 +5.96 0.0 0.0 +5.97 0.0 0.0 +5.98 0.0 0.0 +5.99 0.0 0.0 +6.0 0.0 0.0 +6.01 0.0 0.0 +# l= 1 +0.0 0.0 0.0 +0.01 0.31691115296 0.12833055002 +0.02 0.6327506292 0.25591217486 +0.03 0.94645136007 0.38199950664 +0.04 1.2569554737 0.50585428246 +0.05 1.5632188455 0.62674887265 +0.06 1.8642155902 0.74396977966 +0.07 2.1589424785 0.85682109725 +0.08 2.4464232565 0.96462791932 +0.09 2.7257128515 1.0667396871 +0.1 2.9959014448 1.1625334629 +0.11 3.2561183922 1.2514171183 +0.12 3.5055359765 1.3328324232 +0.13 3.7433729717 1.4062580228 +0.14 3.9688980035 1.4712122871 +0.15 4.1814326898 1.5272560198 +0.16 4.3803545423 1.5739950083 +0.17 4.5650996186 1.6110824023 +0.18 4.7351649045 1.6382209014 +0.19 4.8901104195 1.6551647388 +0.2 5.0295610238 1.6617214382 +0.21 5.1532079278 1.6577533385 +0.22 5.2608098792 1.643178856 +0.23 5.3521940304 1.6179734804 +0.24 5.4272564708 1.5821704817 +0.25 5.4859624266 1.5358613221 +0.26 5.5283461091 1.4791957464 +0.27 5.5545102276 1.4123815589 +0.28 5.5646251474 1.3356840555 +0.29 5.5589277078 1.2494251212 +0.3 5.5377196963 1.1539819761 +0.31 5.5013659832 1.0497855704 +0.32 5.4502923332 0.93731863025 +0.33 5.3849828862 0.81711334783 +0.34 5.3059773248 0.68974871947 +0.35 5.213867786 0.55584758189 +0.36 5.1092954099 0.41607324677 +0.37 4.9929467086 0.27112589469 +0.38 4.8655496976 0.12173868067 +0.39 4.7278696661 -0.031326572994 +0.4 4.5807049337 -0.18728376199 +0.41 4.4248823343 -0.34532720879 +0.42 4.261252569 -0.50463655721 +0.43 4.0906854444 -0.66438183174 +0.44 3.9140650229 -0.82372863041 +0.45 3.732284716 -0.98184341479 +0.46 3.5462423515 -1.137898858 +0.47 3.3568352506 -1.2910792079 +0.48 3.1649554901 -1.4405854704 +0.49 2.9714847265 -1.5856410667 +0.5 2.777289713 -1.7254967608 +0.51 2.5832178663 -1.8594355206 +0.52 2.3900927615 -1.9867774392 +0.53 2.1987094003 -2.1068848753 +0.54 2.009831375 -2.2191655283 +0.55 1.8241857462 -2.323077914 +0.56 1.6424607668 -2.4181337401 +0.57 1.4653019142 -2.5039020266 +0.58 1.293309626 -2.5800113322 +0.59 1.1270361291 -2.6461528615 +0.6 0.96698445288 -2.702081178 +0.61 0.81360483138 -2.7476175365 +0.62 0.66729559136 -2.782648464 +0.63 0.52840007985 -2.8071282798 +0.64 0.39720703887 -2.8210779659 +0.65 0.27395097135 -2.8245839931 +0.66 0.15881083858 -2.8177985492 +0.67 0.051911531399 -2.8009369687 +0.68 -0.04667512666 -2.774275508 +0.69 -0.13692897946 -2.7381486826 +0.7 -0.21888098639 -2.6929468757 +0.71 -0.29261008059 -2.6391119188 +0.72 -0.35824128404 -2.5771334414 +0.73 -0.4159431628 -2.5075445643 +0.74 -0.46592500191 -2.4309173651 +0.75 -0.50843386132 -2.3478577659 +0.76 -0.54375116939 -2.2590003427 +0.77 -0.57218933031 -2.165003103 +0.78 -0.59408818976 -2.0665419618 +0.79 -0.60981137125 -1.9643051699 +0.8 -0.61974254036 -1.8589877089 +0.81 -0.62428163146 -1.7512857073 +0.82 -0.62384107136 -1.6418909303 +0.83 -0.61884203343 -1.5314853966 +0.84 -0.60971075517 -1.4207361709 +0.85 -0.59687495066 -1.3102903826 +0.86 -0.58076034815 -1.2007705148 +0.87 -0.56178738113 -1.0927700081 +0.88 -0.54036805961 -0.98684921825 +0.89 -0.51690304588 -0.88353176365 +0.9 -0.49177895717 -0.78330129473 +0.91 -0.46536591463 -0.68659871302 +0.92 -0.43801535576 -0.59381986204 +0.93 -0.41005812411 -0.50531370759 +0.94 -0.38180284687 -0.42138101873 +0.95 -0.3535346078 -0.3422735548 +0.96 -0.32551478108 -0.26819491238 +0.97 -0.29797801258 -0.1992975509 +0.98 -0.27113326694 -0.13568562865 +0.99 -0.24516387242 -0.077416374488 +1.0 -0.22022736766 -0.024501296954 +1.01 -0.19645591928 0.023091770031 +1.02 -0.17395706242 0.065436168456 +1.03 -0.15281569619 0.10264273095 +1.04 -0.13309095544 0.13486073225 +1.05 -0.11482209959 0.16226970049 +1.06 -0.098028572297 0.18507736826 +1.07 -0.082711657304 0.20351587633 +1.08 -0.068857068133 0.21783821631 +1.09 -0.056433565659 0.22831418337 +1.1 -0.04539878345 0.23522480425 +1.11 -0.035699969459 0.23885904541 +1.12 -0.027275757447 0.23951189031 +1.13 -0.020055775351 0.23748125528 +1.14 -0.013962010699 0.23306834239 +1.15 -0.0089106101841 0.2265749932 +1.16 -0.0048126887966 0.21830334176 +1.17 -0.0015766618338 0.20854486815 +1.18 0.00088988846604 0.19758215521 +1.19 0.0026797111161 0.18568844316 +1.2 0.0038846691896 0.17312405795 +1.21 0.0045919134973 0.16012585041 +1.22 0.0048855679028 0.14691636088 +1.23 0.004847386319 0.13370571947 +1.24 0.0045491243956 0.12066933844 +1.25 0.0040581363722 0.10796761423 +1.26 0.0034376915115 0.095746225461 +1.27 0.0027393344079 0.084115610895 +1.28 0.0020090206457 0.073169616761 +1.29 0.001288444668 0.062988588425 +1.3 0.00060672379237 0.053618373655 +1.31 -1.1366147491e-05 0.045092609726 +1.32 -0.00054582135205 0.037431155887 +1.33 -0.00098759768257 0.030625723663 +1.34 -0.0013297945509 0.024661714734 +1.35 -0.001571493513 0.019508874959 +1.36 -0.0017177187464 0.015122411378 +1.37 -0.0017749428805 0.011452922158 +1.38 -0.0017547314206 0.0084389831708 +1.39 -0.0016697251667 0.0060163835385 +1.4 -0.0015344078581 0.0041157761059 +1.41 -0.0013638321981 0.0026676915805 +1.42 -0.0011732258973 0.0016011633474 +1.43 -0.00097689571102 0.00084843776735 +1.44 -0.00078726131467 0.00034630915422 +1.45 -0.00061636763112 3.2136810107e-05 +1.46 -0.00047079465995 -0.00014219650548 +1.47 -0.00035797587788 -0.00022811628966 +1.48 -0.00027793709055 -0.00025754617817 +1.49 -0.00023007362243 -0.00026174974903 +1.5 -0.00021697247884 -0.00026976292795 +1.51 -0.00020908394519 -0.00026344496824 +1.52 -0.00011615291035 -0.0001468686914 +1.53 -2.9555906996e-06 -3.7564455382e-06 +1.54 1.2960730815e-05 1.6280198529e-05 +1.55 0.0 0.0 +1.56 0.0 0.0 +1.57 0.0 0.0 +1.58 0.0 0.0 +1.59 0.0 0.0 +1.6 0.0 0.0 +1.61 0.0 0.0 +1.62 0.0 0.0 +1.63 0.0 0.0 +1.64 0.0 0.0 +1.65 0.0 0.0 +1.66 0.0 0.0 +1.67 0.0 0.0 +1.68 0.0 0.0 +1.69 0.0 0.0 +1.7 0.0 0.0 +1.71 0.0 0.0 +1.72 0.0 0.0 +1.73 0.0 0.0 +1.74 0.0 0.0 +1.75 0.0 0.0 +1.76 0.0 0.0 +1.77 0.0 0.0 +1.78 0.0 0.0 +1.79 0.0 0.0 +1.8 0.0 0.0 +1.81 0.0 0.0 +1.82 0.0 0.0 +1.83 0.0 0.0 +1.84 0.0 0.0 +1.85 0.0 0.0 +1.86 0.0 0.0 +1.87 0.0 0.0 +1.88 0.0 0.0 +1.89 0.0 0.0 +1.9 0.0 0.0 +1.91 0.0 0.0 +1.92 0.0 0.0 +1.93 0.0 0.0 +1.94 0.0 0.0 +1.95 0.0 0.0 +1.96 0.0 0.0 +1.97 0.0 0.0 +1.98 0.0 0.0 +1.99 0.0 0.0 +2.0 0.0 0.0 +2.01 0.0 0.0 +2.02 0.0 0.0 +2.03 0.0 0.0 +2.04 0.0 0.0 +2.05 0.0 0.0 +2.06 0.0 0.0 +2.07 0.0 0.0 +2.08 0.0 0.0 +2.09 0.0 0.0 +2.1 0.0 0.0 +2.11 0.0 0.0 +2.12 0.0 0.0 +2.13 0.0 0.0 +2.14 0.0 0.0 +2.15 0.0 0.0 +2.16 0.0 0.0 +2.17 0.0 0.0 +2.18 0.0 0.0 +2.19 0.0 0.0 +2.2 0.0 0.0 +2.21 0.0 0.0 +2.22 0.0 0.0 +2.23 0.0 0.0 +2.24 0.0 0.0 +2.25 0.0 0.0 +2.26 0.0 0.0 +2.27 0.0 0.0 +2.28 0.0 0.0 +2.29 0.0 0.0 +2.3 0.0 0.0 +2.31 0.0 0.0 +2.32 0.0 0.0 +2.33 0.0 0.0 +2.34 0.0 0.0 +2.35 0.0 0.0 +2.36 0.0 0.0 +2.37 0.0 0.0 +2.38 0.0 0.0 +2.39 0.0 0.0 +2.4 0.0 0.0 +2.41 0.0 0.0 +2.42 0.0 0.0 +2.43 0.0 0.0 +2.44 0.0 0.0 +2.45 0.0 0.0 +2.46 0.0 0.0 +2.47 0.0 0.0 +2.48 0.0 0.0 +2.49 0.0 0.0 +2.5 0.0 0.0 +2.51 0.0 0.0 +2.52 0.0 0.0 +2.53 0.0 0.0 +2.54 0.0 0.0 +2.55 0.0 0.0 +2.56 0.0 0.0 +2.57 0.0 0.0 +2.58 0.0 0.0 +2.59 0.0 0.0 +2.6 0.0 0.0 +2.61 0.0 0.0 +2.62 0.0 0.0 +2.63 0.0 0.0 +2.64 0.0 0.0 +2.65 0.0 0.0 +2.66 0.0 0.0 +2.67 0.0 0.0 +2.68 0.0 0.0 +2.69 0.0 0.0 +2.7 0.0 0.0 +2.71 0.0 0.0 +2.72 0.0 0.0 +2.73 0.0 0.0 +2.74 0.0 0.0 +2.75 0.0 0.0 +2.76 0.0 0.0 +2.77 0.0 0.0 +2.78 0.0 0.0 +2.79 0.0 0.0 +2.8 0.0 0.0 +2.81 0.0 0.0 +2.82 0.0 0.0 +2.83 0.0 0.0 +2.84 0.0 0.0 +2.85 0.0 0.0 +2.86 0.0 0.0 +2.87 0.0 0.0 +2.88 0.0 0.0 +2.89 0.0 0.0 +2.9 0.0 0.0 +2.91 0.0 0.0 +2.92 0.0 0.0 +2.93 0.0 0.0 +2.94 0.0 0.0 +2.95 0.0 0.0 +2.96 0.0 0.0 +2.97 0.0 0.0 +2.98 0.0 0.0 +2.99 0.0 0.0 +3.0 0.0 0.0 +3.01 0.0 0.0 +3.02 0.0 0.0 +3.03 0.0 0.0 +3.04 0.0 0.0 +3.05 0.0 0.0 +3.06 0.0 0.0 +3.07 0.0 0.0 +3.08 0.0 0.0 +3.09 0.0 0.0 +3.1 0.0 0.0 +3.11 0.0 0.0 +3.12 0.0 0.0 +3.13 0.0 0.0 +3.14 0.0 0.0 +3.15 0.0 0.0 +3.16 0.0 0.0 +3.17 0.0 0.0 +3.18 0.0 0.0 +3.19 0.0 0.0 +3.2 0.0 0.0 +3.21 0.0 0.0 +3.22 0.0 0.0 +3.23 0.0 0.0 +3.24 0.0 0.0 +3.25 0.0 0.0 +3.26 0.0 0.0 +3.27 0.0 0.0 +3.28 0.0 0.0 +3.29 0.0 0.0 +3.3 0.0 0.0 +3.31 0.0 0.0 +3.32 0.0 0.0 +3.33 0.0 0.0 +3.34 0.0 0.0 +3.35 0.0 0.0 +3.36 0.0 0.0 +3.37 0.0 0.0 +3.38 0.0 0.0 +3.39 0.0 0.0 +3.4 0.0 0.0 +3.41 0.0 0.0 +3.42 0.0 0.0 +3.43 0.0 0.0 +3.44 0.0 0.0 +3.45 0.0 0.0 +3.46 0.0 0.0 +3.47 0.0 0.0 +3.48 0.0 0.0 +3.49 0.0 0.0 +3.5 0.0 0.0 +3.51 0.0 0.0 +3.52 0.0 0.0 +3.53 0.0 0.0 +3.54 0.0 0.0 +3.55 0.0 0.0 +3.56 0.0 0.0 +3.57 0.0 0.0 +3.58 0.0 0.0 +3.59 0.0 0.0 +3.6 0.0 0.0 +3.61 0.0 0.0 +3.62 0.0 0.0 +3.63 0.0 0.0 +3.64 0.0 0.0 +3.65 0.0 0.0 +3.66 0.0 0.0 +3.67 0.0 0.0 +3.68 0.0 0.0 +3.69 0.0 0.0 +3.7 0.0 0.0 +3.71 0.0 0.0 +3.72 0.0 0.0 +3.73 0.0 0.0 +3.74 0.0 0.0 +3.75 0.0 0.0 +3.76 0.0 0.0 +3.77 0.0 0.0 +3.78 0.0 0.0 +3.79 0.0 0.0 +3.8 0.0 0.0 +3.81 0.0 0.0 +3.82 0.0 0.0 +3.83 0.0 0.0 +3.84 0.0 0.0 +3.85 0.0 0.0 +3.86 0.0 0.0 +3.87 0.0 0.0 +3.88 0.0 0.0 +3.89 0.0 0.0 +3.9 0.0 0.0 +3.91 0.0 0.0 +3.92 0.0 0.0 +3.93 0.0 0.0 +3.94 0.0 0.0 +3.95 0.0 0.0 +3.96 0.0 0.0 +3.97 0.0 0.0 +3.98 0.0 0.0 +3.99 0.0 0.0 +4.0 0.0 0.0 +4.01 0.0 0.0 +4.02 0.0 0.0 +4.03 0.0 0.0 +4.04 0.0 0.0 +4.05 0.0 0.0 +4.06 0.0 0.0 +4.07 0.0 0.0 +4.08 0.0 0.0 +4.09 0.0 0.0 +4.1 0.0 0.0 +4.11 0.0 0.0 +4.12 0.0 0.0 +4.13 0.0 0.0 +4.14 0.0 0.0 +4.15 0.0 0.0 +4.16 0.0 0.0 +4.17 0.0 0.0 +4.18 0.0 0.0 +4.19 0.0 0.0 +4.2 0.0 0.0 +4.21 0.0 0.0 +4.22 0.0 0.0 +4.23 0.0 0.0 +4.24 0.0 0.0 +4.25 0.0 0.0 +4.26 0.0 0.0 +4.27 0.0 0.0 +4.28 0.0 0.0 +4.29 0.0 0.0 +4.3 0.0 0.0 +4.31 0.0 0.0 +4.32 0.0 0.0 +4.33 0.0 0.0 +4.34 0.0 0.0 +4.35 0.0 0.0 +4.36 0.0 0.0 +4.37 0.0 0.0 +4.38 0.0 0.0 +4.39 0.0 0.0 +4.4 0.0 0.0 +4.41 0.0 0.0 +4.42 0.0 0.0 +4.43 0.0 0.0 +4.44 0.0 0.0 +4.45 0.0 0.0 +4.46 0.0 0.0 +4.47 0.0 0.0 +4.48 0.0 0.0 +4.49 0.0 0.0 +4.5 0.0 0.0 +4.51 0.0 0.0 +4.52 0.0 0.0 +4.53 0.0 0.0 +4.54 0.0 0.0 +4.55 0.0 0.0 +4.56 0.0 0.0 +4.57 0.0 0.0 +4.58 0.0 0.0 +4.59 0.0 0.0 +4.6 0.0 0.0 +4.61 0.0 0.0 +4.62 0.0 0.0 +4.63 0.0 0.0 +4.64 0.0 0.0 +4.65 0.0 0.0 +4.66 0.0 0.0 +4.67 0.0 0.0 +4.68 0.0 0.0 +4.69 0.0 0.0 +4.7 0.0 0.0 +4.71 0.0 0.0 +4.72 0.0 0.0 +4.73 0.0 0.0 +4.74 0.0 0.0 +4.75 0.0 0.0 +4.76 0.0 0.0 +4.77 0.0 0.0 +4.78 0.0 0.0 +4.79 0.0 0.0 +4.8 0.0 0.0 +4.81 0.0 0.0 +4.82 0.0 0.0 +4.83 0.0 0.0 +4.84 0.0 0.0 +4.85 0.0 0.0 +4.86 0.0 0.0 +4.87 0.0 0.0 +4.88 0.0 0.0 +4.89 0.0 0.0 +4.9 0.0 0.0 +4.91 0.0 0.0 +4.92 0.0 0.0 +4.93 0.0 0.0 +4.94 0.0 0.0 +4.95 0.0 0.0 +4.96 0.0 0.0 +4.97 0.0 0.0 +4.98 0.0 0.0 +4.99 0.0 0.0 +5.0 0.0 0.0 +5.01 0.0 0.0 +5.02 0.0 0.0 +5.03 0.0 0.0 +5.04 0.0 0.0 +5.05 0.0 0.0 +5.06 0.0 0.0 +5.07 0.0 0.0 +5.08 0.0 0.0 +5.09 0.0 0.0 +5.1 0.0 0.0 +5.11 0.0 0.0 +5.12 0.0 0.0 +5.13 0.0 0.0 +5.14 0.0 0.0 +5.15 0.0 0.0 +5.16 0.0 0.0 +5.17 0.0 0.0 +5.18 0.0 0.0 +5.19 0.0 0.0 +5.2 0.0 0.0 +5.21 0.0 0.0 +5.22 0.0 0.0 +5.23 0.0 0.0 +5.24 0.0 0.0 +5.25 0.0 0.0 +5.26 0.0 0.0 +5.27 0.0 0.0 +5.28 0.0 0.0 +5.29 0.0 0.0 +5.3 0.0 0.0 +5.31 0.0 0.0 +5.32 0.0 0.0 +5.33 0.0 0.0 +5.34 0.0 0.0 +5.35 0.0 0.0 +5.36 0.0 0.0 +5.37 0.0 0.0 +5.38 0.0 0.0 +5.39 0.0 0.0 +5.4 0.0 0.0 +5.41 0.0 0.0 +5.42 0.0 0.0 +5.43 0.0 0.0 +5.44 0.0 0.0 +5.45 0.0 0.0 +5.46 0.0 0.0 +5.47 0.0 0.0 +5.48 0.0 0.0 +5.49 0.0 0.0 +5.5 0.0 0.0 +5.51 0.0 0.0 +5.52 0.0 0.0 +5.53 0.0 0.0 +5.54 0.0 0.0 +5.55 0.0 0.0 +5.56 0.0 0.0 +5.57 0.0 0.0 +5.58 0.0 0.0 +5.59 0.0 0.0 +5.6 0.0 0.0 +5.61 0.0 0.0 +5.62 0.0 0.0 +5.63 0.0 0.0 +5.64 0.0 0.0 +5.65 0.0 0.0 +5.66 0.0 0.0 +5.67 0.0 0.0 +5.68 0.0 0.0 +5.69 0.0 0.0 +5.7 0.0 0.0 +5.71 0.0 0.0 +5.72 0.0 0.0 +5.73 0.0 0.0 +5.74 0.0 0.0 +5.75 0.0 0.0 +5.76 0.0 0.0 +5.77 0.0 0.0 +5.78 0.0 0.0 +5.79 0.0 0.0 +5.8 0.0 0.0 +5.81 0.0 0.0 +5.82 0.0 0.0 +5.83 0.0 0.0 +5.84 0.0 0.0 +5.85 0.0 0.0 +5.86 0.0 0.0 +5.87 0.0 0.0 +5.88 0.0 0.0 +5.89 0.0 0.0 +5.9 0.0 0.0 +5.91 0.0 0.0 +5.92 0.0 0.0 +5.93 0.0 0.0 +5.94 0.0 0.0 +5.95 0.0 0.0 +5.96 0.0 0.0 +5.97 0.0 0.0 +5.98 0.0 0.0 +5.99 0.0 0.0 +6.0 0.0 0.0 +6.01 0.0 0.0 +# local +0.0 -9.7953079069E+00 +0.01 -1.0400311409E+01 +0.02 -1.0888406632E+01 +0.03 -1.1200435759E+01 +0.04 -1.1277240974E+01 +0.05 -1.1209255048E+01 +0.06 -1.1082805427E+01 +0.07 -1.0947620831E+01 +0.08 -1.0824723658E+01 +0.09 -1.0719296734E+01 +0.1 -1.0629697231E+01 +0.11 -1.0552200672E+01 +0.12 -1.0483027171E+01 +0.13 -1.0419016416E+01 +0.14 -1.0357741770E+01 +0.15 -1.0297434680E+01 +0.16 -1.0236868006E+01 +0.17 -1.0175247464E+01 +0.18 -1.0112121292E+01 +0.19 -1.0047305980E+01 +0.2 -9.9808234246E+00 +0.21 -9.9128463822E+00 +0.22 -9.8436503701E+00 +0.23 -9.7735719697E+00 +0.24 -9.7029738376E+00 +0.25 -9.6322169180E+00 +0.26 -9.5616399874E+00 +0.27 -9.4915460360E+00 +0.28 -9.4221946690E+00 +0.29 -9.3537991265E+00 +0.3 -9.2865265731E+00 +0.31 -9.2205003077E+00 +0.32 -9.1558027948E+00 +0.33 -9.0924785976E+00 +0.34 -9.0305368263E+00 +0.35 -8.9699529387E+00 +0.36 -8.9106696470E+00 +0.37 -8.8525976536E+00 +0.38 -8.7956163315E+00 +0.39 -8.7395746729E+00 +0.4 -8.6842935256E+00 +0.41 -8.6295690309E+00 +0.42 -8.5751777523E+00 +0.43 -8.5208835233E+00 +0.44 -8.4664458092E+00 +0.45 -8.4116291209E+00 +0.46 -8.3562127953E+00 +0.47 -8.3000003168E+00 +0.48 -8.2428271911E+00 +0.49 -8.1845673552E+00 +0.5 -8.1251364881E+00 +0.51 -8.0644926540E+00 +0.52 -8.0026346953E+00 +0.53 -7.9395985751E+00 +0.54 -7.8754509157E+00 +0.55 -7.8102832830E+00 +0.56 -7.7442045501E+00 +0.57 -7.6773346132E+00 +0.58 -7.6097981896E+00 +0.59 -7.5417198577E+00 +0.6 -7.4732201798E+00 +0.61 -7.4044126046E+00 +0.62 -7.3354020632E+00 +0.63 -7.2662835375E+00 +0.64 -7.1971420148E+00 +0.65 -7.1280529085E+00 +0.66 -7.0590822441E+00 +0.67 -6.9902877566E+00 +0.68 -6.9217198555E+00 +0.69 -6.8534226258E+00 +0.7 -6.7854345116E+00 +0.71 -6.7177895103E+00 +0.72 -6.6505178739E+00 +0.73 -6.5836468129E+00 +0.74 -6.5172011142E+00 +0.75 -6.4512036077E+00 +0.76 -6.3856756197E+00 +0.77 -6.3206373121E+00 +0.78 -6.2561079364E+00 +0.79 -6.1921060357E+00 +0.8 -6.1286495929E+00 +0.81 -6.0657561346E+00 +0.82 -6.0034427989E+00 +0.83 -5.9417263711E+00 +0.84 -5.8806232961E+00 +0.85 -5.8201496694E+00 +0.86 -5.7603212132E+00 +0.87 -5.7011532385E+00 +0.88 -5.6426605986E+00 +0.89 -5.5848576333E+00 +0.9 -5.5277581083E+00 +0.91 -5.4713751498E+00 +0.92 -5.4157211748E+00 +0.93 -5.3608078208E+00 +0.94 -5.3066458721E+00 +0.95 -5.2532451864E+00 +0.96 -5.2006146500E+00 +0.97 -5.1487620332E+00 +0.98 -5.0976939521E+00 +0.99 -5.0474158034E+00 +1.0 -4.9979316802E+00 +1.01 -4.9492442981E+00 +1.02 -4.9013549271E+00 +1.03 -4.8542634146E+00 +1.04 -4.8079678026E+00 +1.05 -4.7624645403E+00 +1.06 -4.7177483280E+00 +1.07 -4.6738120477E+00 +1.08 -4.6306468941E+00 +1.09 -4.5882416969E+00 +1.1 -4.5465833716E+00 +1.11 -4.5056567463E+00 +1.12 -4.4654448539E+00 +1.13 -4.4259298481E+00 +1.14 -4.3870946833E+00 +1.15 -4.3489228924E+00 +1.16 -4.3113985506E+00 +1.17 -4.2745064577E+00 +1.18 -4.2382318993E+00 +1.19 -4.2025606063E+00 +1.2 -4.1674787716E+00 +1.21 -4.1329731249E+00 +1.22 -4.0990307195E+00 +1.23 -4.0656388980E+00 +1.24 -4.0327855198E+00 +1.25 -4.0004586688E+00 +1.26 -3.9686466511E+00 +1.27 -3.9373381834E+00 +1.28 -3.9065221765E+00 +1.29 -3.8761876992E+00 +1.3 -3.8463241824E+00 +1.31 -3.8169211999E+00 +1.32 -3.7879684651E+00 +1.33 -3.7594560262E+00 +1.34 -3.7313740213E+00 +1.35 -3.7037127811E+00 +1.36 -3.6764629142E+00 +1.37 -3.6496150905E+00 +1.38 -3.6231602853E+00 +1.39 -3.5970896520E+00 +1.4 -3.5713944446E+00 +1.41 -3.5460663178E+00 +1.42 -3.5210969494E+00 +1.43 -3.4964783706E+00 +1.44 -3.4722028189E+00 +1.45 -3.4482626548E+00 +1.46 -3.4246506993E+00 +1.47 -3.4013597263E+00 +1.48 -3.3783830415E+00 +1.49 -3.3557139795E+00 +1.5 -3.3333461989E+00 +1.51 -3.3112736518E+00 +1.52 -3.2894903707E+00 +1.53 -3.2679910124E+00 +1.54 -3.2467701693E+00 +1.55 -3.2258229108E+00 +1.56 -3.2051442058E+00 +1.57 -3.1847290067E+00 +1.58 -3.1645723052E+00 +1.59 -3.1446693226E+00 +1.6 -3.1250152668E+00 +1.61 -3.1056053810E+00 +1.62 -3.0864350882E+00 +1.63 -3.0674999620E+00 +1.64 -3.0487957324E+00 +1.65 -3.0303182431E+00 +1.66 -3.0120634156E+00 +1.67 -2.9940271521E+00 +1.68 -2.9762055892E+00 +1.69 -2.9585948974E+00 +1.7 -2.9411913974E+00 +1.71 -2.9239914167E+00 +1.72 -2.9069914394E+00 +1.73 -2.8901879713E+00 +1.74 -2.8735776456E+00 +1.75 -2.8571571408E+00 +1.76 -2.8409232236E+00 +1.77 -2.8248727423E+00 +1.78 -2.8090025834E+00 +1.79 -2.7933097620E+00 +1.8 -2.7777912717E+00 +1.81 -2.7624442908E+00 +1.82 -2.7472659091E+00 +1.83 -2.7322534622E+00 +1.84 -2.7174041320E+00 +1.85 -2.7027153831E+00 +1.86 -2.6881845381E+00 +1.87 -2.6738091286E+00 +1.88 -2.6595866351E+00 +1.89 -2.6455146440E+00 +1.9 -2.6315907931E+00 +1.91 -2.6178127142E+00 +1.92 -2.6041781993E+00 +1.93 -2.5906849203E+00 +1.94 -2.5773307925E+00 +1.95 -2.5641135986E+00 +1.96 -2.5510312912E+00 +1.97 -2.5380818020E+00 +1.98 -2.5252630984E+00 +1.99 -2.5125732639E+00 +2.0 -2.5000102769E+00 +2.01 -2.4875723360E+00 +2.02 -2.4752575202E+00 +2.03 -2.4630640399E+00 +2.04 -2.4509901208E+00 +2.05 -2.4390339638E+00 +2.06 -2.4271939251E+00 +2.07 -2.4154682479E+00 +2.08 -2.4038553365E+00 +2.09 -2.3923535617E+00 +2.1 -2.3809613030E+00 +2.11 -2.3696770644E+00 +2.12 -2.3584992450E+00 +2.13 -2.3474264006E+00 +2.14 -2.3364570512E+00 +2.15 -2.3255897148E+00 +2.16 -2.3148230361E+00 +2.17 -2.3041555626E+00 +2.18 -2.2935859647E+00 +2.19 -2.2831129150E+00 +2.2 -2.2727350358E+00 +2.21 -2.2624511056E+00 +2.22 -2.2522598146E+00 +2.23 -2.2421599154E+00 +2.24 -2.2321502244E+00 +2.25 -2.2222294734E+00 +2.26 -2.2123965330E+00 +2.27 -2.2026502440E+00 +2.28 -2.1929894132E+00 +2.29 -2.1834129847E+00 +2.3 -2.1739198255E+00 +2.31 -2.1645088428E+00 +2.32 -2.1551790170E+00 +2.33 -2.1459292535E+00 +2.34 -2.1367585500E+00 +2.35 -2.1276659213E+00 +2.36 -2.1186503127E+00 +2.37 -2.1097108018E+00 +2.38 -2.1008464331E+00 +2.39 -2.0920562046E+00 +2.4 -2.0833392519E+00 +2.41 -2.0746946489E+00 +2.42 -2.0661214596E+00 +2.43 -2.0576188552E+00 +2.44 -2.0491859464E+00 +2.45 -2.0408218556E+00 +2.46 -2.0325257888E+00 +2.47 -2.0242968934E+00 +2.48 -2.0161343429E+00 +2.49 -2.0080373776E+00 +2.5 -2.0000051821E+00 +2.51 -1.9920369736E+00 +2.52 -1.9841320267E+00 +2.53 -1.9762895631E+00 +2.54 -1.9685088371E+00 +2.55 -1.9607891574E+00 +2.56 -1.9531297830E+00 +2.57 -1.9455299990E+00 +2.58 -1.9379891482E+00 +2.59 -1.9305065270E+00 +2.6 -1.9230814452E+00 +2.61 -1.9157132799E+00 +2.62 -1.9084013649E+00 +2.63 -1.9011450292E+00 +2.64 -1.8939436842E+00 +2.65 -1.8867967011E+00 +2.66 -1.8797034234E+00 +2.67 -1.8726632960E+00 +2.68 -1.8656757242E+00 +2.69 -1.8587400729E+00 +2.7 -1.8518558088E+00 +2.71 -1.8450223684E+00 +2.72 -1.8382391535E+00 +2.73 -1.8315056327E+00 +2.74 -1.8248212788E+00 +2.75 -1.8181855313E+00 +2.76 -1.8115978559E+00 +2.77 -1.8050577615E+00 +2.78 -1.7985647255E+00 +2.79 -1.7921182080E+00 +2.8 -1.7857177520E+00 +2.81 -1.7793628661E+00 +2.82 -1.7730530256E+00 +2.83 -1.7667877827E+00 +2.84 -1.7605666767E+00 +2.85 -1.7543892206E+00 +2.86 -1.7482549508E+00 +2.87 -1.7421634433E+00 +2.88 -1.7361142500E+00 +2.89 -1.7301068892E+00 +2.9 -1.7241409715E+00 +2.91 -1.7182160709E+00 +2.92 -1.7123317389E+00 +2.93 -1.7064875664E+00 +2.94 -1.7006831640E+00 +2.95 -1.6949181227E+00 +2.96 -1.6891920072E+00 +2.97 -1.6835044635E+00 +2.98 -1.6778551046E+00 +2.99 -1.6722435230E+00 +3.0 -1.6666693437E+00 +3.01 -1.6611322150E+00 +3.02 -1.6556317677E+00 +3.03 -1.6501675996E+00 +3.04 -1.6447393891E+00 +3.05 -1.6393467861E+00 +3.06 -1.6339894265E+00 +3.07 -1.6286669522E+00 +3.08 -1.6233790509E+00 +3.09 -1.6181253872E+00 +3.1 -1.6129056026E+00 +3.11 -1.6077193828E+00 +3.12 -1.6025664196E+00 +3.13 -1.5974463932E+00 +3.14 -1.5923589527E+00 +3.15 -1.5873038213E+00 +3.16 -1.5822806958E+00 +3.17 -1.5772892645E+00 +3.18 -1.5723292057E+00 +3.19 -1.5674002544E+00 +3.2 -1.5625021199E+00 +3.21 -1.5576344984E+00 +3.22 -1.5527970977E+00 +3.23 -1.5479896601E+00 +3.24 -1.5432119082E+00 +3.25 -1.5384635475E+00 +3.26 -1.5337443106E+00 +3.27 -1.5290539473E+00 +3.28 -1.5243921936E+00 +3.29 -1.5197587649E+00 +3.3 -1.5151534139E+00 +3.31 -1.5105758996E+00 +3.32 -1.5060259698E+00 +3.33 -1.5015033514E+00 +3.34 -1.4970078128E+00 +3.35 -1.4925391227E+00 +3.36 -1.4880970407E+00 +3.37 -1.4836813055E+00 +3.38 -1.4792916976E+00 +3.39 -1.4749279962E+00 +3.4 -1.4705899722E+00 +3.41 -1.4662773767E+00 +3.42 -1.4619899986E+00 +3.43 -1.4577276285E+00 +3.44 -1.4534900480E+00 +3.45 -1.4492770216E+00 +3.46 -1.4450883429E+00 +3.47 -1.4409238146E+00 +3.48 -1.4367832287E+00 +3.49 -1.4326663637E+00 +3.5 -1.4285730146E+00 +3.51 -1.4245029971E+00 +3.52 -1.4204561129E+00 +3.53 -1.4164321554E+00 +3.54 -1.4124309178E+00 +3.55 -1.4084522296E+00 +3.56 -1.4044959016E+00 +3.57 -1.4005617428E+00 +3.58 -1.3966495421E+00 +3.59 -1.3927591429E+00 +3.6 -1.3888903647E+00 +3.61 -1.3850430270E+00 +3.62 -1.3812169299E+00 +3.63 -1.3774119129E+00 +3.64 -1.3736278098E+00 +3.65 -1.3698644486E+00 +3.66 -1.3661216462E+00 +3.67 -1.3623992310E+00 +3.68 -1.3586970529E+00 +3.69 -1.3550149479E+00 +3.7 -1.3513527506E+00 +3.71 -1.3477102763E+00 +3.72 -1.3440873909E+00 +3.73 -1.3404839380E+00 +3.74 -1.3368997612E+00 +3.75 -1.3333346896E+00 +3.76 -1.3297885764E+00 +3.77 -1.3262612817E+00 +3.78 -1.3227526565E+00 +3.79 -1.3192625488E+00 +3.8 -1.3157907932E+00 +3.81 -1.3123372674E+00 +3.82 -1.3089018292E+00 +3.83 -1.3054843366E+00 +3.84 -1.3020846342E+00 +3.85 -1.2987025869E+00 +3.86 -1.2953380691E+00 +3.87 -1.2919909450E+00 +3.88 -1.2886610792E+00 +3.89 -1.2853483151E+00 +3.9 -1.2820525436E+00 +3.91 -1.2787736362E+00 +3.92 -1.2755114636E+00 +3.93 -1.2722658895E+00 +3.94 -1.2690367790E+00 +3.95 -1.2658240235E+00 +3.96 -1.2626274997E+00 +3.97 -1.2594470844E+00 +3.98 -1.2562826417E+00 +3.99 -1.2531340562E+00 +4.0 -1.2500012187E+00 +4.01 -1.2468840119E+00 +4.02 -1.2437823182E+00 +4.03 -1.2406960031E+00 +4.04 -1.2376249674E+00 +4.05 -1.2345691024E+00 +4.06 -1.2315282961E+00 +4.07 -1.2285024364E+00 +4.08 -1.2254913921E+00 +4.09 -1.2224950757E+00 +4.1 -1.2195133805E+00 +4.11 -1.2165461998E+00 +4.12 -1.2135934242E+00 +4.13 -1.2106549331E+00 +4.14 -1.2077306417E+00 +4.15 -1.2048204481E+00 +4.16 -1.2019242507E+00 +4.17 -1.1990419438E+00 +4.18 -1.1961734148E+00 +4.19 -1.1933185819E+00 +4.2 -1.1904773482E+00 +4.21 -1.1876496167E+00 +4.22 -1.1848352867E+00 +4.23 -1.1820342502E+00 +4.24 -1.1792464301E+00 +4.25 -1.1764717336E+00 +4.26 -1.1737100682E+00 +4.27 -1.1709613393E+00 +4.28 -1.1682254410E+00 +4.29 -1.1655023010E+00 +4.3 -1.1627918312E+00 +4.31 -1.1600939433E+00 +4.32 -1.1574085490E+00 +4.33 -1.1547355427E+00 +4.34 -1.1520748576E+00 +4.35 -1.1494264097E+00 +4.36 -1.1467901148E+00 +4.37 -1.1441658889E+00 +4.38 -1.1415536337E+00 +4.39 -1.1389532790E+00 +4.4 -1.1363647480E+00 +4.41 -1.1337879605E+00 +4.42 -1.1312228362E+00 +4.43 -1.1286692854E+00 +4.44 -1.1261272320E+00 +4.45 -1.1235966071E+00 +4.46 -1.1210773343E+00 +4.47 -1.1185693371E+00 +4.48 -1.1160725347E+00 +4.49 -1.1135868436E+00 +4.5 -1.1111122031E+00 +4.51 -1.1086485404E+00 +4.52 -1.1061957826E+00 +4.53 -1.1037538568E+00 +4.54 -1.1013226758E+00 +4.55 -1.0989021820E+00 +4.56 -1.0964923080E+00 +4.57 -1.0940929841E+00 +4.58 -1.0917041410E+00 +4.59 -1.0893257018E+00 +4.6 -1.0869575968E+00 +4.61 -1.0845997688E+00 +4.62 -1.0822521513E+00 +4.63 -1.0799146781E+00 +4.64 -1.0775872830E+00 +4.65 -1.0752698851E+00 +4.66 -1.0729624349E+00 +4.67 -1.0706648699E+00 +4.68 -1.0683771270E+00 +4.69 -1.0660991430E+00 +4.7 -1.0638308489E+00 +4.71 -1.0615721790E+00 +4.72 -1.0593230825E+00 +4.73 -1.0570834992E+00 +4.74 -1.0548533689E+00 +4.75 -1.0526326312E+00 +4.76 -1.0504212151E+00 +4.77 -1.0482190695E+00 +4.78 -1.0460261407E+00 +4.79 -1.0438423714E+00 +4.8 -1.0416677041E+00 +4.81 -1.0395020809E+00 +4.82 -1.0373454313E+00 +4.83 -1.0351977142E+00 +4.84 -1.0330588748E+00 +4.85 -1.0309288585E+00 +4.86 -1.0288076105E+00 +4.87 -1.0266950722E+00 +4.88 -1.0245911830E+00 +4.89 -1.0224959010E+00 +4.9 -1.0204091741E+00 +4.91 -1.0183309501E+00 +4.92 -1.0162611767E+00 +4.93 -1.0141997958E+00 +4.94 -1.0121467543E+00 +4.95 -1.0101020102E+00 +4.96 -1.0080655139E+00 +4.97 -1.0060372156E+00 +4.98 -1.0040170655E+00 +4.99 -1.0020050068E+00 +5.0 -1.0000009915E+00 +5.01 -9.9800497853E-01 +5.02 -9.9601692053E-01 +5.03 -9.9403677005E-01 +5.04 -9.9206447965E-01 +5.05 -9.9009999481E-01 +5.06 -9.8814327013E-01 +5.07 -9.8619426653E-01 +5.08 -9.8425293881E-01 +5.09 -9.8231924171E-01 +5.1 -9.8039313002E-01 +5.11 -9.7847455227E-01 +5.12 -9.7656346372E-01 +5.13 -9.7465982794E-01 +5.14 -9.7276360183E-01 +5.15 -9.7087474224E-01 +5.16 -9.6899320607E-01 +5.17 -9.6711894568E-01 +5.18 -9.6525191502E-01 +5.19 -9.6339208106E-01 +5.2 -9.6153940271E-01 +5.21 -9.5969383884E-01 +5.22 -9.5785534834E-01 +5.23 -9.5602388818E-01 +5.24 -9.5419940913E-01 +5.25 -9.5238188227E-01 +5.26 -9.5057126840E-01 +5.27 -9.4876752834E-01 +5.28 -9.4697062287E-01 +5.29 -9.4518051281E-01 +5.3 -9.4339714846E-01 +5.31 -9.4162050123E-01 +5.32 -9.3985053520E-01 +5.33 -9.3808721299E-01 +5.34 -9.3633049724E-01 +5.35 -9.3458035056E-01 +5.36 -9.3283672935E-01 +5.37 -9.3109959783E-01 +5.38 -9.2936892590E-01 +5.39 -9.2764467793E-01 +5.4 -9.2592681829E-01 +5.41 -9.2421531134E-01 +5.42 -9.2251012024E-01 +5.43 -9.2081120058E-01 +5.44 -9.1911852851E-01 +5.45 -9.1743207008E-01 +5.46 -9.1575179130E-01 +5.47 -9.1407765821E-01 +5.48 -9.1240963683E-01 +5.49 -9.1074768644E-01 +5.5 -9.0909177628E-01 +5.51 -9.0744187842E-01 +5.52 -9.0579796048E-01 +5.53 -9.0415999006E-01 +5.54 -9.0252793478E-01 +5.55 -9.0090176200E-01 +5.56 -8.9928142906E-01 +5.57 -8.9766691560E-01 +5.58 -8.9605819074E-01 +5.59 -8.9445522360E-01 +5.6 -8.9285798330E-01 +5.61 -8.9126643897E-01 +5.62 -8.8968055584E-01 +5.63 -8.8810030060E-01 +5.64 -8.8652565061E-01 +5.65 -8.8495657644E-01 +5.66 -8.8339304863E-01 +5.67 -8.8183503776E-01 +5.68 -8.8028251439E-01 +5.69 -8.7873544257E-01 +5.7 -8.7719379623E-01 +5.71 -8.7565755125E-01 +5.72 -8.7412667957E-01 +5.73 -8.7260115311E-01 +5.74 -8.7108094382E-01 +5.75 -8.6956602361E-01 +5.76 -8.6805635630E-01 +5.77 -8.6655192092E-01 +5.78 -8.6505269276E-01 +5.79 -8.6355864504E-01 +5.8 -8.6206975101E-01 +5.81 -8.6058598390E-01 +5.82 -8.5910731697E-01 +5.83 -8.5763371456E-01 +5.84 -8.5616515888E-01 +5.85 -8.5470162543E-01 +5.86 -8.5324308868E-01 +5.87 -8.5178952313E-01 +5.88 -8.5034090325E-01 +5.89 -8.4889720355E-01 +5.9 -8.4745838969E-01 +5.91 -8.4602444519E-01 +5.92 -8.4459534656E-01 +5.93 -8.4317106947E-01 +5.94 -8.4175158961E-01 +5.95 -8.4033688264E-01 +5.96 -8.3892692424E-01 +5.97 -8.3752168212E-01 +5.98 -8.3612113938E-01 +5.99 -8.3472527429E-01 +6.0 -8.3333406366E-01 +6.01 -8.3194748432E-01 From bc897c6fbb7dade0d0991ef1d8691866eb3aa270 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Mon, 22 Jul 2024 09:48:03 -0400 Subject: [PATCH 09/17] Update 2-pyridone example (#259) --- examples/2-pyridone/2-pyridone.cfg | 26 +++++++------------------- examples/2-pyridone/README | 10 +++++----- examples/2-pyridone/coords.in | 25 ++++++++++++------------- 3 files changed, 24 insertions(+), 37 deletions(-) diff --git a/examples/2-pyridone/2-pyridone.cfg b/examples/2-pyridone/2-pyridone.cfg index 82457d9b..23f53496 100644 --- a/examples/2-pyridone/2-pyridone.cfg +++ b/examples/2-pyridone/2-pyridone.cfg @@ -1,5 +1,6 @@ verbosity=1 xcFunctional=PBE +FDtype=4th [Mesh] nx=128 ny=128 @@ -12,30 +13,17 @@ lx=24. ly=24. lz=12. [Potentials] -pseudopotential=pseudo.N_soft_pbe -pseudopotential=pseudo.O_pbe -pseudopotential=pseudo.C_pbe -pseudopotential=pseudo.H_pbe -[Poisson] -solver=MG -max_steps_initial=10 -max_steps=10 -bcx=periodic -bcy=periodic -bcz=periodic +pseudopotential=pseudo.N_ONCV_PBE_SG15 +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.C_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 [Run] type=QUENCH [Quench] -step_length=1.5 max_steps=200 atol=1.e-7 -num_lin_iterations=2 -ortho_freq=20 -interval_print_residual=10 [Orbitals] initial_type=Gaussian initial_width=1.5 -[ProjectedMatrices] -solver=exact -[Restart] -output_type=single_file +[Poisson] +FDtype=4th diff --git a/examples/2-pyridone/README b/examples/2-pyridone/README index 351b7012..663f5855 100644 --- a/examples/2-pyridone/README +++ b/examples/2-pyridone/README @@ -1,6 +1,6 @@ -ln -s ../../potentials/pseudo.N_soft_pbe -ln -s ../../potentials/pseudo.O_pbe -ln -s ../../potentials/pseudo.C_pbe -ln -s ../../potentials/pseudo.H_pbe +ln -s ../../potentials/pseudo.N_ONCV_PBE_SG15 +ln -s ../../potentials/pseudo.O_ONCV_PBE_SG15 +ln -s ../../potentials/pseudo.C_ONCV_PBE_SG15 +ln -s ../../potentials/pseudo.H_ONCV_PBE_SG15 -srun -ppdebug -n32 /nfs/tmp2/jeanluc/SVN/MGmol/mgmol/branches/bin/mgmol-pel -c 2-pyridone.cfg -l lrs.in -i coords.in +srun -ppdebug -n32 mgmol-opt -c 2-pyridone.cfg -l lrs.in -i coords.in diff --git a/examples/2-pyridone/coords.in b/examples/2-pyridone/coords.in index 911b5cdf..390ff902 100644 --- a/examples/2-pyridone/coords.in +++ b/examples/2-pyridone/coords.in @@ -1,13 +1,12 @@ -N1 1 10.0818 7.9149 6. 1 -O1 2 14.2787 7.9126 6. 1 -C1 3 7.9625 9.1867 6. 0 -C2 3 7.8388 11.8070 6. 0 -C3 3 10.1050 13.1610 6. 0 -C4 3 12.3830 11.8680 6. 0 -C5 3 12.2810 9.2255 6. 0 -H1 4 6.2219 8.0878 6. 1 -H2 4 6.0028 12.7142 6. 1 -H3 4 10.0754 15.2192 6. 1 -H4 4 14.2295 12.7482 6. 1 -H5 4 10.2332 6.0038 6. 1 - +N1 1 10.0818 7.9149 6. +O1 2 14.2787 7.9126 6. +C1 3 7.9625 9.1867 6. +C2 3 7.8388 11.8070 6. +C3 3 10.1050 13.1610 6. +C4 3 12.3830 11.8680 6. +C5 3 12.2810 9.2255 6. +H1 4 6.2219 8.0878 6. +H2 4 6.0028 12.7142 6. +H3 4 10.0754 15.2192 6. +H4 4 14.2295 12.7482 6. +H5 4 10.2332 6.0038 6. From 3c111219eb07558c4f6f2784b0b9ddd387742957 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Mon, 29 Jul 2024 11:45:47 -0400 Subject: [PATCH 10/17] Adjust verbosity in some functions (#260) --- src/DavidsonSolver.cc | 12 ++++++------ src/DensityMatrix.cc | 4 ++-- src/Hartree_CG.cc | 4 +++- src/MVPSolver.cc | 8 +++++--- src/ProjectedMatricesInterface.h | 2 +- 5 files changed, 17 insertions(+), 13 deletions(-) diff --git a/src/DavidsonSolver.cc b/src/DavidsonSolver.cc index c53e11c0..0dc88ef4 100644 --- a/src/DavidsonSolver.cc +++ b/src/DavidsonSolver.cc @@ -406,7 +406,7 @@ int DavidsonSolver::solve( // orthonormalization over outer iterations (MD or geometry optimization) // jlf, 01/02/2021 if (mmpi.PE0() && ct.verbose > 1) - os_ << "Orthonormalize wavefunctions ate start of DavidsonSolver" + os_ << "Orthonormalize wavefunctions at start of DavidsonSolver" << std::endl; orbitals.orthonormalizeLoewdin(false, nullptr, false); @@ -418,7 +418,7 @@ int DavidsonSolver::solve( ProjectedMatrices* proj_matN = dynamic_cast*>( orbitals.getProjMatrices()); - if (mmpi.PE0() && ct.verbose > 1) + if (mmpi.PE0() && ct.verbose > 0) { os_ << "###########################" << std::endl; os_ << "DavidsonSolver -> Iteration " << outer_it << std::endl; @@ -498,7 +498,7 @@ int DavidsonSolver::solve( // if( onpe0 )os_<<"Matrices..."<printMatrices(os_); - ts0 = evalEntropy(projmatrices, true, os_); + ts0 = evalEntropy(projmatrices, (ct.verbose > 1), os_); e0 = energy_->evaluateTotal( ts0, projmatrices, orbitals, printE, os_); @@ -547,7 +547,7 @@ int DavidsonSolver::solve( proj_mat2N_->assignBlocksH(h11, h12, h21, h22); proj_mat2N_->setHB2H(); - ts0 = evalEntropy(proj_mat2N_.get(), true, os_); + ts0 = evalEntropy(proj_mat2N_.get(), (ct.verbose > 1), os_); e0 = energy_->evaluateTotal( ts0, proj_mat2N_.get(), orbitals, printE, os_); } @@ -627,7 +627,7 @@ int DavidsonSolver::solve( // line minimization beta = minQuadPolynomial(e0, e1, de0, (ct.verbose > 2), os_); - if (mmpi.PE0() && ct.verbose > 1) + if (mmpi.PE0() && ct.verbose > 0) { os_ << std::setprecision(12); os_ << "ts1=" << ts1 << std::endl; @@ -690,7 +690,7 @@ int DavidsonSolver::solve( os_ << "Total occupations for top half states=" << std::setprecision(15) << tot << std::endl; } - if (mmpi.PE0() && ct.verbose > 0) + if (mmpi.PE0() && ct.verbose > 1) { os_ << std::setprecision(15) << "Last level occupancy = " << eval[numst_] << std::endl; diff --git a/src/DensityMatrix.cc b/src/DensityMatrix.cc index 32877d01..87b30224 100644 --- a/src/DensityMatrix.cc +++ b/src/DensityMatrix.cc @@ -118,13 +118,13 @@ void DensityMatrix::build( template void DensityMatrix::build(const int new_orbitals_index) { - //#ifdef PRINT_OPERATIONS MGmol_MPI& mmpi = *(MGmol_MPI::instance()); +#ifdef PRINT_OPERATIONS if (mmpi.PE0()) std::cout << "DensityMatrix::build() for diagonal occupation..." << std::endl; - //#endif +#endif if (!occ_uptodate_ && mmpi.instancePE0()) std::cout << "Warning: occupations not up to date to build DM!!!" << std::endl; diff --git a/src/Hartree_CG.cc b/src/Hartree_CG.cc index db156a9f..f9bf3860 100644 --- a/src/Hartree_CG.cc +++ b/src/Hartree_CG.cc @@ -90,8 +90,10 @@ void Hartree_CG::solve( double residual_reduction = poisson_solver_->getResidualReduction(); double final_residual = poisson_solver_->getFinalResidual(); + const bool large_residual + = (residual_reduction > 1.e-3 || final_residual > 1.e-3); - if (onpe0) + if (onpe0 && (large_residual || ct.verbose > 1)) (*MPIdata::sout) << setprecision(2) << scientific << "Hartree_CG: residual reduction = " << residual_reduction diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc index 63ddf690..45237e36 100644 --- a/src/MVPSolver.cc +++ b/src/MVPSolver.cc @@ -59,7 +59,8 @@ MVPSolver::MVPSolver(MPI_Comm comm, std::ostream& os, if (onpe0 && ct.verbose > 0) { os_ << "MVPSolver..." << std::endl; - if (use_old_dm_) os_ << "MVPSolver uses old DM..." << std::endl; + if (use_old_dm_ && ct.verbose > 1) + os_ << "MVPSolver uses old DM..." << std::endl; } rho_ = rho; @@ -252,8 +253,9 @@ int MVPSolver::solve(OrbitalsType& orbitals) dmInit = proj_mat_work_->dm(); } - const double ts0 = evalEntropyMVP(current_proj_mat, true, os_); - const double e0 = energy_->evaluateTotal( + const double ts0 + = evalEntropyMVP(current_proj_mat, (ct.verbose > 1), os_); + const double e0 = energy_->evaluateTotal( ts0, current_proj_mat, orbitals, printE, os_); MatrixType target("target", numst_, numst_); diff --git a/src/ProjectedMatricesInterface.h b/src/ProjectedMatricesInterface.h index 4041ba76..653b386d 100644 --- a/src/ProjectedMatricesInterface.h +++ b/src/ProjectedMatricesInterface.h @@ -71,7 +71,7 @@ class ProjectedMatricesInterface : public ChebyshevApproximationFunction Control& ct = *(Control::instance()); nel_ = ct.getNelSpin(); MGmol_MPI& mmpi = *(MGmol_MPI::instance()); - if (mmpi.PE0()) + if (mmpi.PE0() && ct.verbose > 1) std::cout << "ProjectedMatricesInterface: nel_=" << nel_ << std::endl; }; From 35ba70fb9eaef3cc3f5e089cc59ef761a9cad0b2 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Mon, 29 Jul 2024 12:47:56 -0400 Subject: [PATCH 11/17] Add new example: pinned H2O (#261) --- examples/PinnedH2O/README | 3 +++ examples/PinnedH2O/coords.in | 3 +++ examples/PinnedH2O/mgmol.cfg | 35 +++++++++++++++++++++++++++++++++++ 3 files changed, 41 insertions(+) create mode 100644 examples/PinnedH2O/README create mode 100644 examples/PinnedH2O/coords.in create mode 100644 examples/PinnedH2O/mgmol.cfg diff --git a/examples/PinnedH2O/README b/examples/PinnedH2O/README new file mode 100644 index 00000000..12ca9138 --- /dev/null +++ b/examples/PinnedH2O/README @@ -0,0 +1,3 @@ +Single water molecule MD, with Oxigen atom pinned to (0,0,0) coordinate: + +mpirun -np 4 mgmol-opt -c mgmol.cfg -i coords.in diff --git a/examples/PinnedH2O/coords.in b/examples/PinnedH2O/coords.in new file mode 100644 index 00000000..307e170c --- /dev/null +++ b/examples/PinnedH2O/coords.in @@ -0,0 +1,3 @@ +O1 1 0.00 0.00 0.00 0 +H1 2 -0.45 1.42 -1.07 1 +H2 2 -0.45 -1.48 -0.97 1 diff --git a/examples/PinnedH2O/mgmol.cfg b/examples/PinnedH2O/mgmol.cfg new file mode 100644 index 00000000..32325e47 --- /dev/null +++ b/examples/PinnedH2O/mgmol.cfg @@ -0,0 +1,35 @@ +verbosity=1 +xcFunctional=PBE +FDtype=Mehrstellen +[Mesh] +nx=64 +ny=64 +nz=64 +[Domain] +ox=-6. +oy=-6. +oz=-6. +lx=12. +ly=12. +lz=12. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=MD +[MD] +num_steps=50 +dt=40. +thermostat=ON +[Thermostat] +type=Berendsen +temperature=1000. +relax_time=800. +[Quench] +max_steps=100 +atol=1.e-8 +[Orbitals] +initial_type=Random +initial_width=2. +[Restart] +output_level=0 From bde85069f875ef3e410b7b5a24269bb914a2e7eb Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Tue, 30 Jul 2024 11:45:29 -0400 Subject: [PATCH 12/17] Print out eigenvalues out of MVP solver (#262) Previously, the wrong eigenvalues (0) were printed out because eigenvalues outside solver were not up-to-date. --- src/MVPSolver.cc | 1 + src/ProjectedMatrices.h | 6 ++++++ tests/MVP/test.py | 24 ++++++++++++++++++++++++ 3 files changed, 31 insertions(+) diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc index 45237e36..7feb564d 100644 --- a/src/MVPSolver.cc +++ b/src/MVPSolver.cc @@ -347,6 +347,7 @@ int MVPSolver::solve(OrbitalsType& orbitals) = dynamic_cast*>( orbitals.getProjMatrices()); projmatrices->setDM(*work_, orbitals.getIterativeIndex()); + projmatrices->setEigenvalues(proj_mat_work_->getEigenvalues()); } // Generate new density diff --git a/src/ProjectedMatrices.h b/src/ProjectedMatrices.h index 702637f4..6725b827 100644 --- a/src/ProjectedMatrices.h +++ b/src/ProjectedMatrices.h @@ -409,6 +409,12 @@ class ProjectedMatrices : public ProjectedMatricesInterface { dm_->setMatrix(mat, orbitals_index); } + void setEigenvalues(const std::vector& eigenvalues) + { + memcpy(eigenvalues_.data(), eigenvalues.data(), + eigenvalues.size() * sizeof(double)); + } + const std::vector& getEigenvalues() const { return eigenvalues_; } void computeGenEigenInterval(std::vector& interval, const int maxits, const double padding = 0.01); DensityMatrix& getDM() { return *dm_; } diff --git a/tests/MVP/test.py b/tests/MVP/test.py index 708fd45c..bd5d708e 100644 --- a/tests/MVP/test.py +++ b/tests/MVP/test.py @@ -70,6 +70,30 @@ print("Found HDF5 error") sys.exit(1) +flag = 0 +eigenvalues=[] +for line in lines: + if line.count(b'FERMI'): + flag = 0 + if flag==1: + words=line.split() + for w in words: + eigenvalues.append(eval(w)) + if line.count(b'Eigenvalues'): + flag = 1 + eigenvalues=[] + +print(eigenvalues) +tol = 1.e-4 +eigenvalue0 = -0.208 +if abs(eigenvalues[0]-eigenvalue0)>tol: + print("Expected eigenvalue 0 to be {}".format(eigenvalue0)) + sys.exit(1) +eigenvalue50 = 0.208 +if abs(eigenvalues[50]-eigenvalue50)>tol: + print("Expected eigenvalue 50 to be {}".format(eigenvalue50)) + sys.exit(1) + niterations = len(energies) print("MVP solver ran for {} iterations".format(niterations)) if niterations>180: From 4bbb64e9c9138f56b6d8377e9480a515d22c1145 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Wed, 31 Jul 2024 09:04:46 -0400 Subject: [PATCH 13/17] Add Li2 example with local GTH potential (#263) --- examples/Li2GTH/README | 11 + examples/Li2GTH/davidson.cfg | 46 +++++ examples/Li2GTH/li2.xyz | 4 + potentials/pseudo.Li_GTH_PBE | 332 +++++++++++++++++++++++++++++++ util/generateLiLocalGTHpseudo.py | 65 ++++++ 5 files changed, 458 insertions(+) create mode 100644 examples/Li2GTH/README create mode 100644 examples/Li2GTH/davidson.cfg create mode 100644 examples/Li2GTH/li2.xyz create mode 100644 potentials/pseudo.Li_GTH_PBE create mode 100644 util/generateLiLocalGTHpseudo.py diff --git a/examples/Li2GTH/README b/examples/Li2GTH/README new file mode 100644 index 00000000..f411a01e --- /dev/null +++ b/examples/Li2GTH/README @@ -0,0 +1,11 @@ +#Run MGmol +mpirun -np 4 mgmol-opt -c davidson.cfg -i li2.xyz > davidson.out + +#extract visit files from HDF5 file +python read_hdf5.py -bov li2.h5 Vtotal + +python read_hdf5.py -bov li2.h50 Function0002 + +Note: the eigenfunctions corresponding to the 3 occupied states are stored +in functions 2, 3 and 4 because they are the result of diagonalizing the DM +in Davidson algorithm, not H. diff --git a/examples/Li2GTH/davidson.cfg b/examples/Li2GTH/davidson.cfg new file mode 100644 index 00000000..d211b812 --- /dev/null +++ b/examples/Li2GTH/davidson.cfg @@ -0,0 +1,46 @@ +verbosity=2 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx=96 +ny=96 +nz=96 +[Domain] +ox=-9. +oy=-9. +oz=-9. +lx=18. +ly=18. +lz=18. +[Potentials] +pseudopotential=pseudo.Li_GTH_PBE +[Poisson] +solver=CG +bcx=0 +bcy=0 +bcz=0 +[Run] +type=QUENCH +[Quench] +solver=Davidson +max_steps=200 +atol=1.e-8 +preconditioner_num_levels=3 +[Orbitals] +nempty=2 +initial_type=random +temperature=300. +bcx=0 +bcy=0 +bcz=0 +[ProjectedMatrices] +solver=exact +[DensityMatrix] +nb_inner_it=1 +mixing=1. +[Restart] +output_level=3 +output_filename=li2.h5 +output_type=single_file +[Coloring] +scope=global diff --git a/examples/Li2GTH/li2.xyz b/examples/Li2GTH/li2.xyz new file mode 100644 index 00000000..ccf3a3f1 --- /dev/null +++ b/examples/Li2GTH/li2.xyz @@ -0,0 +1,4 @@ +2 +Experimental geometry (Menconi and Tozer, Chem. Phys. Lett. 2002) +Li 1.3365 0.0 0.0 +Li -1.3365 0.0 0.0 diff --git a/potentials/pseudo.Li_GTH_PBE b/potentials/pseudo.Li_GTH_PBE new file mode 100644 index 00000000..7da944e0 --- /dev/null +++ b/potentials/pseudo.Li_GTH_PBE @@ -0,0 +1,332 @@ +# Goedecker, Teter, Hutter, Phys. rev. B 54 (3), 1996 +# Parameters from M. Krack, Theor. Chem. Acc. 114, 2005 +# Short description of the species type. One line only! +LithiumGTH_PBE +# +White +#radii of balls and covalent bonds +0.4 1.0 +# Nlcc flag +0 +# Atomic number +3 +# Atomic mass +3.0 +# Number of valence electrons +3.0 +# Gaussian core charge parameter rc +1. +# Number of potentials +1 +# l-value for state which is local +0 0 +# Local potential radius +3. +# Non-local potential radius +3. +# number of points in radial grid +301 +# log mesh parameter +0. +# radial grid, reference state, and potential for l=0 +0.0 -20.06528920602149 +0.01 -20.054252432503986 +0.02 -20.021181949571293 +0.03 -19.966197051499293 +0.04 -19.889495835016746 +0.05 -19.791354114111485 +0.06 -19.672123911517346 +0.07 -19.53223153752974 +0.08 -19.372175269684018 +0.09 -19.192522649601603 +0.1 -18.993907415941983 +0.11 -18.777026094871008 +0.12 -18.54263427174855 +0.13 -18.291542569831762 +0.14 -18.02461236366824 +0.15 -17.742751256501393 +0.16 -17.446908352415903 +0.17 -17.138069355105003 +0.18 -16.81725152603506 +0.19 -16.485498535412276 +0.2 -16.14387523971814 +0.21 -15.7934624196748 +0.22 -15.435351512330925 +0.23 -15.07063937052802 +0.24 -14.700423082322942 +0.25 -14.325794882014955 +0.26 -13.947837183265412 +0.27 -13.567617763419314 +0.28 -13.186185126555326 +0.29 -12.80456407102134 +0.3 -12.423751485274781 +0.31 -12.0447123937598 +0.32 -11.668376272337989 +0.33 -11.295633650466765 +0.34 -10.927333014911554 +0.35 -10.564278027307214 +0.36 -10.207225065372501 +0.37 -9.856881095051628 +0.38 -9.513901878330472 +0.39 -9.178890518973706 +0.4 -8.852396345973204 +0.41 -8.534914132107685 +0.42 -8.226883642706916 +0.43 -7.928689507508834 +0.44 -7.64066140641021 +0.45 -7.363074557955612 +0.46 -7.096150497598659 +0.47 -6.840058131114424 +0.48 -6.594915047052536 +0.49 -6.3607890708040955 +0.5 -6.1377000417180545 +0.51 -5.925621793748402 +0.52 -5.724484319343899 +0.53 -5.534176095708496 +0.54 -5.354546552160849 +0.55 -5.1854086571031885 +0.56 -5.026541603068161 +0.57 -4.877693568441651 +0.58 -4.73858453475222 +0.59 -4.60890913886544 +0.6 -4.488339540014026 +0.61 -4.376528282322026 +0.62 -4.273111134331415 +0.63 -4.17770988800058 +0.64 -4.089935100703225 +0.65 -4.009388764900398 +0.66 -3.935666891373983 +0.67 -3.868361993183916 +0.68 -3.807065458829906 +0.69 -3.751369804448444 +0.7 -3.7008707962442955 +0.71 -3.655169435730016 +0.72 -3.6138738017150143 +0.73 -3.5766007443361536 +0.74 -3.5429774277437494 +0.75 -3.512642719340314 +0.76 -3.4852484247051203 +0.77 -3.4604603685176145 +0.78 -3.4379593229091037 +0.79 -3.4174417857190456 +0.8 -3.3986206121036857 +0.81 -3.38122550383682 +0.82 -3.3650033614510035 +0.83 -3.3497185050904883 +0.84 -3.335152770582715 +0.85 -3.3211054877824484 +0.86 -3.307393348702202 +0.87 -3.293850173314738 +0.88 -3.2803265812005225 +0.89 -3.2666895774169027 +0.9 -3.2528220610898493 +0.91 -3.2386222652767795 +0.92 -3.224003136624447 +0.93 -3.208891663253622 +0.94 -3.1932281591472704 +0.95 -3.1769655131062664 +0.96 -3.160068410071886 +0.97 -3.142512532302873 +0.98 -3.1242837475424827 +0.99 -3.10537729092335 +1.0 -3.085796946940759 +1.01 -3.065554237383791 +1.02 -3.0446676206540504 +1.03 -3.0231617074287267 +1.04 -3.001066497143526 +1.05 -2.978416639286402 +1.06 -2.9552507230094722 +1.07 -2.931610598088238 +1.08 -2.9075407297881393 +1.09 -2.8830875897419252 +1.1 -2.8582990845006977 +1.11 -2.833224022999286 +1.12 -2.8079116237754826 +1.13 -2.7824110624044724 +1.14 -2.7567710592563195 +1.15 -2.7310395073568614 +1.16 -2.705263139831839 +1.17 -2.6794872361411723 +1.18 -2.653755366065358 +1.19 -2.628109170188997 +1.2 -2.60258817543732 +1.21 -2.57722964405972 +1.22 -2.552068454319101 +1.23 -2.527137011036365 +1.24 -2.50246518405455 +1.25 -2.4780802726257742 +1.26 -2.4540069936848927 +1.27 -2.430267491955173 +1.28 -2.4068813698318525 +1.29 -2.383865735007489 +1.3 -2.361235263837057 +1.31 -2.339002278488971 +1.32 -2.31717683598915 +1.33 -2.2957668273371157 +1.34 -2.274778084954434 +1.35 -2.2542144968149556 +1.36 -2.2340781257018256 +1.37 -2.2143693321366866 +1.38 -2.1950868996304913 +1.39 -2.176228161011706 +1.4 -2.1577891246951135 +1.41 -2.1397645998619574 +1.42 -2.1221483196287054 +1.43 -2.1049330613864585 +1.44 -2.0881107635950986 +1.45 -2.07167263841506 +1.46 -2.055609279654423 +1.47 -2.03991076559945 +1.48 -2.0245667563822187 +1.49 -2.0095665856193814 +1.5 -1.9948993461310054 +1.51 -1.9805539696177794 +1.52 -1.96651930023847 +1.53 -1.95278416208739 +1.54 -1.9393374206237712 +1.55 -1.9261680381514323 +1.56 -1.9132651234880866 +1.57 -1.900617975999247 +1.58 -1.8882161242021371 +1.59 -1.876049359170528 +1.6 -1.8641077629922906 +1.61 -1.8523817325479044 +1.62 -1.8408619988905475 +1.63 -1.8295396425169552 +1.64 -1.8184061048233324 +1.65 -1.807453196042526 +1.66 -1.796673099957719 +1.67 -1.786058375684427 +1.68 -1.7756019568068189 +1.69 -1.7652971481466952 +1.7 -1.7551376204340732 +1.71 -1.7451174031375227 +1.72 -1.7352308757004553 +1.73 -1.7254727574166897 +1.74 -1.7158380961650346 +1.75 -1.7063222562085625 +1.76 -1.696920905249867 +1.77 -1.6876300009190652 +1.78 -1.6784457768568064 +1.79 -1.6693647285401818 +1.8 -1.6603835989853395 +1.81 -1.6514993644468863 +1.82 -1.6427092202208944 +1.83 -1.6340105666456095 +1.84 -1.625400995381811 +1.85 -1.6168782760432898 +1.86 -1.608440343237089 +1.87 -1.6000852840630355 +1.88 -1.5918113261126972 +1.89 -1.5836168259992183 +1.9 -1.5755002584415392 +1.91 -1.567460205919256 +1.92 -1.5594953489078422 +1.93 -1.5516044566980827 +1.94 -1.54378637879836 +1.95 -1.5360400369138487 +1.96 -1.5283644174926794 +1.97 -1.5207585648257078 +1.98 -1.513221574683613 +1.99 -1.5057525884726553 +2.0 -1.4983507878884341 +2.01 -1.4910153900454863 +2.02 -1.4837456430593685 +2.03 -1.476540822057088 +2.04 -1.4694002255912106 +2.05 -1.4623231724327952 +2.06 -1.4553089987182684 +2.07 -1.4483570554256544 +2.08 -1.441466706155941 +2.09 -1.4346373251959794 +2.1 -1.4278682958400157 +2.11 -1.4211590089477848 +2.12 -1.4145088617179962 +2.13 -1.407917256657042 +2.14 -1.4013836007237732 +2.15 -1.3949073046322673 +2.16 -1.38848778229559 +2.17 -1.3821244503946528 +2.18 -1.3758167280573612 +2.19 -1.3695640366343098 +2.2 -1.3633657995583632 +2.21 -1.3572214422764575 +2.22 -1.3511303922429727 +2.23 -1.3450920789649687 +2.24 -1.3391059340904876 +2.25 -1.3331713915319947 +2.26 -1.3272878876178418 +2.27 -1.321454861265409 +2.28 -1.3156717541702987 +2.29 -1.3099380110066239 +2.3 -1.304253079634056 +2.31 -1.2986164113078702 +2.32 -1.2930274608887506 +2.33 -1.2874856870495992 +2.34 -1.2819905524770356 +2.35 -1.2765415240656612 +2.36 -1.271138073103527 +2.37 -1.2657796754475625 +2.38 -1.260465811688003 +2.39 -1.2551959673011142 +2.4 -1.2499696327897274 +2.41 -1.244786303811294 +2.42 -1.2396454812933442 +2.43 -1.2345466715363633 +2.44 -1.229489386304243 +2.45 -1.224473142902546 +2.46 -1.219497464244928 +2.47 -1.214561878908113 +2.48 -1.2096659211758896 +2.49 -1.204809131072622 +2.5 -1.1999910543868184 +2.51 -1.1952112426853008 +2.52 -1.1904692533185557 +2.53 -1.1857646494178307 +2.54 -1.1810969998845513 +2.55 -1.1764658793726301 +2.56 -1.1718708682642154 +2.57 -1.1673115526394273 +2.58 -1.1627875242406034 +2.59 -1.158298380431558 +2.6 -1.1538437241523383 +2.61 -1.149423163869939 +2.62 -1.1450363135254051 +2.63 -1.1406827924777474 +2.64 -1.1363622254450418 +2.65 -1.1320742424430907 +2.66 -1.1278184787219725 +2.67 -1.1235945747008054 +2.68 -1.1194021759010047 +2.69 -1.1152409328783188 +2.7 -1.1111105011538749 +2.71 -1.107010541144477 +2.72 -1.1029407180923525 +2.73 -1.0989007019945478 +2.74 -1.094890167532131 +2.75 -1.0909087939993736 +2.76 -1.0869562652330373 +2.77 -1.0830322695418986 +2.78 -1.079136499636622 +2.79 -1.0752686525600812 +2.8 -1.0714284296182166 +2.81 -1.067615536311507 +2.82 -1.0638296822671227 +2.83 -1.0600705811718185 +2.84 -1.0563379507056208 +2.85 -1.0526315124763437 +2.86 -1.0489509919549826 +2.87 -1.0452961184120013 +2.88 -1.0416666248545523 +2.89 -1.0380622479646329 +2.9 -1.0344827280382114 +2.91 -1.0309278089253149 +2.92 -1.0273972379711012 +2.93 -1.0238907659579113 +2.94 -1.0204081470483046 +2.95 -1.0169491387290766 +2.96 -1.013513501756257 +2.97 -1.0101010001010782 +2.98 -1.0067114008969125 +2.99 -1.0033444743871645 +3.0 -0.9999999938741118 diff --git a/util/generateLiLocalGTHpseudo.py b/util/generateLiLocalGTHpseudo.py new file mode 100644 index 00000000..1b553370 --- /dev/null +++ b/util/generateLiLocalGTHpseudo.py @@ -0,0 +1,65 @@ +# Generate local potential according to +# Goedecker, Teter, Hutter, Phys. rev. B 54 (3), 1996 +# Parameters from M. Krack, Theor. Chem. Acc. 114, 2005 +from math import exp, erf, sqrt, pi + +#coefficients for H +rloc = 0.4 +c1 = -14.081155 +c2 = 9.626220 +c3 = -1.783616 +c4 = 0.085152 +zion = 3. +anumber = 3 +name = "LithiumGTH_PBE" +mass = 3. + +def radialfunction(r): + alpha = (r/rloc)**2 + val = exp(-0.5*alpha)*(c1+c2*alpha+c3*alpha*alpha+c4*alpha*alpha*alpha) + if r>1.e-8: + val = val - zion*erf(r/(sqrt(2.)*rloc))/r + else: + #print("special case for r = {}".format(r)) + val = val -zion*sqrt(2.)/(sqrt(pi)*rloc) + + return val + +npts = 301 + +#header +print("# Short description of the species type. One line only!") +print(name) +print("#") +print("White") +print("#radii of balls and covalent bonds") +print("0.4 1.0") +print("# Nlcc flag") +print("0") +print("# Atomic number") +print(anumber) +print("# Atomic mass") +print(mass) +print("# Number of valence electrons") +print(zion) +print("# Gaussian core charge parameter rc") +print("1.") +print("# Number of potentials") +print("1") +print("# l-value for state which is local") +print("0 0") +print("# Local potential radius") +print("3.") +print("# Non-local potential radius") +print("3.") +print("# number of points in radial grid") +print(npts) +print("# log mesh parameter") +print("0.") +print("# radial grid, reference state, and potential for l=0") + +#potential +for i in range(npts): + r = round(0.01*i,2) + f = radialfunction(r) + print("{} {}".format(r,f)) From 0ead573485e01ef40833106876864dcdd322ab37 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Wed, 31 Jul 2024 15:26:06 -0400 Subject: [PATCH 14/17] Fix LBFGS termination when converged (#264) --- src/lbfgsrlx.cc | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/lbfgsrlx.cc b/src/lbfgsrlx.cc index b2aff795..6698eef0 100644 --- a/src/lbfgsrlx.cc +++ b/src/lbfgsrlx.cc @@ -77,6 +77,7 @@ void MGmol::lbfgsrlx(OrbitalsType** orbitals, Ions& ions) lbfgs.computeForces(); + // check for convergence int flag_convF = lbfgs.checkTolForces(ct.tol_forces); int conv = 0; @@ -87,10 +88,8 @@ void MGmol::lbfgsrlx(OrbitalsType** orbitals, Ions& ions) os_ << endl << endl << "LBFGS: convergence in forces has been achieved. " - "stopping ..." << endl; } - conv = 1; } else { @@ -116,10 +115,17 @@ void MGmol::lbfgsrlx(OrbitalsType** orbitals, Ions& ions) // Write down positions and displacements ions.printPositions(os_); + // early termination if convergence achieved + if (flag_convF) + { + if (onpe0) os_ << "LBFGS: stopping ..." << endl; + break; + } + + // reset iterations if last step failed if (conv != 0) { if (onpe0) os_ << "LBFGS Geometry optimization stopped" << endl; - // break; if (onpe0) os_ << "LBFGS Geometry optimization reset" << endl; lbfgs.reset(ct.dt); } From b7e75a24ede116397a26edd7ce18388cd9024267 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Wed, 31 Jul 2024 15:55:00 -0400 Subject: [PATCH 15/17] Remove unused code to extrapolate rho (#265) --- src/md.cc | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/src/md.cc b/src/md.cc index 5e016941..f7fe934d 100644 --- a/src/md.cc +++ b/src/md.cc @@ -576,18 +576,8 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) ct.steps = md_iteration_; -#if EXTRAPOLATE_RHO - if (onpe0) os_ << "Extrapolate rho..." << std::endl; - rho_->axpyRhoc(-1., rhoc_); - rho_->extrapolate(); -#endif - moveVnuc(ions); -#if EXTRAPOLATE_RHO - rho_->axpyRhoc(1., rhoc_); -#endif - if (!small_move) { orbitals_extrapol_->clearOldOrbitals(); @@ -669,8 +659,7 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) } template -void MGmol::loadRestartFile( - const std::string filename) +void MGmol::loadRestartFile(const std::string filename) { MGmol_MPI& mmpi(*(MGmol_MPI::instance())); Control& ct = *(Control::instance()); @@ -688,7 +677,9 @@ void MGmol::loadRestartFile( if (ierr < 0) { if (onpe0) - (*MPIdata::serr) << "loadRestartFile: failed to read the restart file." << std::endl; + (*MPIdata::serr) + << "loadRestartFile: failed to read the restart file." + << std::endl; global_exit(0); } From 7bcd787a3a9fd5344c73adedb69cc9f91a2b613f Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Thu, 1 Aug 2024 14:01:04 -0400 Subject: [PATCH 16/17] Fix and test EnergyAndForces interface (#266) * Atomic potentials were not updated when atomic positions were changed * Added test to make sure energies and forces are the same after positions move by one mesh spacing --- src/MGmol.cc | 2 + src/MGmol.h | 5 +- tests/CMakeLists.txt | 11 ++ tests/EnergyAndForces/coords.in | 2 + tests/EnergyAndForces/lrs.in | 6 + tests/EnergyAndForces/mgmol.cfg | 33 ++++ tests/EnergyAndForces/test.py | 87 ++++++++ tests/EnergyAndForces/testEnergyAndForces.cc | 198 +++++++++++++++++++ 8 files changed, 343 insertions(+), 1 deletion(-) create mode 100644 tests/EnergyAndForces/coords.in create mode 100644 tests/EnergyAndForces/lrs.in create mode 100644 tests/EnergyAndForces/mgmol.cfg create mode 100755 tests/EnergyAndForces/test.py create mode 100644 tests/EnergyAndForces/testEnergyAndForces.cc diff --git a/src/MGmol.cc b/src/MGmol.cc index e1f65d21..682a0e05 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -1428,6 +1428,8 @@ double MGmol::evaluateEnergyAndForces( ions_->setPositions(tau, atnumbers); + moveVnuc(*ions_); + double eks = 0.; quench(current_orbitals_, *ions_, ct.max_electronic_steps, 20, eks); diff --git a/src/MGmol.h b/src/MGmol.h index 12096a5b..8a3ae1bb 100644 --- a/src/MGmol.h +++ b/src/MGmol.h @@ -180,7 +180,10 @@ class MGmol : public MGmolInterface /* access functions */ OrbitalsType* getOrbitals() { return current_orbitals_; } - std::shared_ptr> getHamiltonian() { return hamiltonian_; } + std::shared_ptr> getHamiltonian() + { + return hamiltonian_; + } void run() override; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index fe7b571b..3f221912 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -246,6 +246,8 @@ add_executable(testDensityMatrix ${CMAKE_SOURCE_DIR}/src/tools/mgmol_mpi_tools.cc ${CMAKE_SOURCE_DIR}/src/tools/random.cc ${CMAKE_SOURCE_DIR}/tests/ut_magma_main.cc) +add_executable(testEnergyAndForces + ${CMAKE_SOURCE_DIR}/tests/EnergyAndForces/testEnergyAndForces.cc) if(${MAGMA_FOUND}) add_executable(testOpenmpOffload @@ -334,6 +336,14 @@ add_test(NAME testGramMatrix add_test(NAME testDensityMatrix COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/testDensityMatrix) +add_test(NAME testEnergyAndForces + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/EnergyAndForces/test.py + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/testEnergyAndForces + ${CMAKE_CURRENT_SOURCE_DIR}/EnergyAndForces/mgmol.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/EnergyAndForces/coords.in + ${CMAKE_CURRENT_SOURCE_DIR}/EnergyAndForces/lrs.in + ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) if(${MAGMA_FOUND}) add_test(NAME testOpenmpOffload @@ -518,6 +528,7 @@ target_link_libraries(testBlacsContext PRIVATE ${SCALAPACK_LIBRARIES} ${BLAS_LIBRARIES} MPI::MPI_CXX) target_link_libraries(testSuperSampling PRIVATE MPI::MPI_CXX) target_link_libraries(testDirectionalReduce PRIVATE MPI::MPI_CXX) +target_link_libraries(testEnergyAndForces PRIVATE mgmol_src) if(${MAGMA_FOUND}) target_link_libraries(testDistVector PRIVATE ${SCALAPACK_LIBRARIES} diff --git a/tests/EnergyAndForces/coords.in b/tests/EnergyAndForces/coords.in new file mode 100644 index 00000000..cda898b9 --- /dev/null +++ b/tests/EnergyAndForces/coords.in @@ -0,0 +1,2 @@ +N1 1 0. 0. -1.0345 +N2 1 0. 0. 1.0345 diff --git a/tests/EnergyAndForces/lrs.in b/tests/EnergyAndForces/lrs.in new file mode 100644 index 00000000..9136a500 --- /dev/null +++ b/tests/EnergyAndForces/lrs.in @@ -0,0 +1,6 @@ +0.00 0.7 1.4 +0.00 -0.7 -1.4 +0.00 0.7 -1.4 +0.00 -0.7 1.4 +0.46 0.0 0.0 +-0.46 0.0 0.0 diff --git a/tests/EnergyAndForces/mgmol.cfg b/tests/EnergyAndForces/mgmol.cfg new file mode 100644 index 00000000..cf23889e --- /dev/null +++ b/tests/EnergyAndForces/mgmol.cfg @@ -0,0 +1,33 @@ +verbosity=1 +xcFunctional=LDA +FDtype=Mehrstellen +[Mesh] +nx=64 +ny=64 +nz=64 +[Domain] +ox=-6. +oy=-6. +oz=-6. +lx=12. +ly=12. +lz=12. +[Potentials] +pseudopotential=pseudo.N_ONCVPSP_LDA +[Run] +type=QUENCH +[Quench] +solver=PSD +max_steps=100 +atol=1.e-9 +step_length=2. +ortho_freq=10 +[Orbitals] +initial_type=Gaussian +initial_width=1.5 +temperature=10. +nempty=1 +[Restart] +output_level=0 +[DensityMatrix] +mixing=0.5 diff --git a/tests/EnergyAndForces/test.py b/tests/EnergyAndForces/test.py new file mode 100755 index 00000000..e066bfc6 --- /dev/null +++ b/tests/EnergyAndForces/test.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python +import sys +import os +import subprocess +import string + +print("Test EnergyAndForces...") + +nargs=len(sys.argv) + +mpicmd = sys.argv[1]+" "+sys.argv[2]+" "+sys.argv[3] +for i in range(4,nargs-6): + mpicmd = mpicmd + " "+sys.argv[i] +print("MPI run command: {}".format(mpicmd)) + +exe = sys.argv[nargs-5] +inp = sys.argv[nargs-4] +coords = sys.argv[nargs-3] +print("coordinates file: %s"%coords) +lrs = sys.argv[-2] + +#create links to potentials files +dst = 'pseudo.N_ONCVPSP_LDA' +src = sys.argv[-1] + '/' + dst + +if not os.path.exists(dst): + print("Create link to %s"%dst) + os.symlink(src, dst) + +#run +command = "{} {} -c {} -i {} -l {}".format(mpicmd,exe,inp,coords,lrs) +print("Run command: {}".format(command)) +output = subprocess.check_output(command,shell=True) +lines=output.split(b'\n') + +#analyse output +energies=[] +for line in lines: + if line.count(b'%%'): + print(line) + words=line.split() + words=words[5].split(b',')[0] + energy = words.decode() + if line.count(b'achieved'): + energies.append(energy) + +flag=0 +forces=[] +for line in lines: + if flag>0: + print(line) + words=line.split(b' ') + forces.append(words[1].decode()) + forces.append(words[2].decode()) + forces.append(words[3].decode()) + flag=flag-1 + if line.count(b'Forces:'): + flag=2 + + +print("Check energies...") +print( energies ) +if len(energies)<2: + print("Expected two converged energies") + sys.exit(1) + +tol = 1.e-6 +diff=eval(energies[1])-eval(energies[0]) +print(diff) +if abs(diff)>tol: + print("Energies differ: {} vs {} !!!".format(energies[0],energies[1])) + sys.exit(1) + +print("Check forces...") +print(forces) +flag=0 +for i in range(6): + diff=eval(forces[i+6])-eval(forces[i]) + print(diff) + if abs(diff)>1.e-3: + print("Forces difference larger than tol") + flag=1 +if flag>0: + sys.exit(1) + +print("Test SUCCESSFUL!") +sys.exit(0) diff --git a/tests/EnergyAndForces/testEnergyAndForces.cc b/tests/EnergyAndForces/testEnergyAndForces.cc new file mode 100644 index 00000000..c8f7ab95 --- /dev/null +++ b/tests/EnergyAndForces/testEnergyAndForces.cc @@ -0,0 +1,198 @@ +// Copyright (c) 2017, Lawrence Livermore National Security, LLC and +// UT-Battelle, LLC. +// Produced at the Lawrence Livermore National Laboratory and the Oak Ridge +// National Laboratory. +// LLNL-CODE-743438 +// All rights reserved. +// 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 "Control.h" +#include "ExtendedGridOrbitals.h" +#include "LocGridOrbitals.h" +#include "MGmol.h" +#include "MGmol_MPI.h" +#include "MPIdata.h" +#include "mgmol_run.h" + +#include +#include +#include +#include + +#include +namespace po = boost::program_options; + +int main(int argc, char** argv) +{ + int mpirc = MPI_Init(&argc, &argv); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Initialization failed!!!" << std::endl; + MPI_Abort(MPI_COMM_WORLD, 0); + } + + MPI_Comm comm = MPI_COMM_WORLD; + + /* + * Initialize general things, like magma, openmp, IO, ... + */ + mgmol_init(comm); + + /* + * read runtime parameters + */ + std::string input_filename(""); + std::string lrs_filename; + std::string constraints_filename(""); + + float total_spin = 0.; + bool with_spin = false; + + po::variables_map vm; + + // read from PE0 only + if (MPIdata::onpe0) + { + read_config(argc, argv, vm, input_filename, lrs_filename, + constraints_filename, total_spin, with_spin); + } + + MGmol_MPI::setup(comm, std::cout, with_spin); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + MPI_Comm global_comm = mmpi.commGlobal(); + + /* + * Setup control struct with run time parameters + */ + Control::setup(global_comm, with_spin, total_spin); + Control& ct = *(Control::instance()); + + ct.setOptions(vm); + + int ret = ct.checkOptions(); + if (ret < 0) return ret; + + mmpi.bcastGlobal(input_filename); + mmpi.bcastGlobal(lrs_filename); + + // Enter main scope + { + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "Construct MGmol object..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + + MGmolInterface* mgmol; + if (ct.isLocMode()) + 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); + + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "MGmol setup..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + mgmol->setup(); + + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "Setup done..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + + // here we just use the atomic positions read in and used + // to initialize MGmol + std::vector positions; + mgmol->getAtomicPositions(positions); + std::vector anumbers; + mgmol->getAtomicNumbers(anumbers); + if (MPIdata::onpe0) + { + std::cout << "Positions:" << std::endl; + std::vector::iterator ita = anumbers.begin(); + for (std::vector::iterator it = positions.begin(); + it != positions.end(); it += 3) + { + std::cout << *ita; + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + ita++; + } + } + + // compute energy and forces using all MPI tasks + // expect positions to be replicated on all MPI tasks + std::vector forces; + double eks + = mgmol->evaluateEnergyAndForces(positions, anumbers, forces); + + // print out results + if (MPIdata::onpe0) + { + std::cout << "Eks: " << eks << std::endl; + std::cout << "Forces:" << std::endl; + for (std::vector::iterator it = forces.begin(); + it != forces.end(); it += 3) + { + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + } + } + + // move atoms one mesh spacing in all directions + Mesh* mymesh = Mesh::instance(); + const pb::Grid& mygrid = mymesh->grid(); + const double hspacing[3] + = { mygrid.hgrid(0), mygrid.hgrid(1), mygrid.hgrid(2) }; + int i = 0; + for (auto& pos : positions) + { + pos += hspacing[i % 3]; + i++; + } + + // evaluate energy and forces again + eks = mgmol->evaluateEnergyAndForces(positions, anumbers, forces); + + // print out results + if (MPIdata::onpe0) + { + std::cout << "Eks: " << eks << std::endl; + std::cout << "Forces:" << std::endl; + for (std::vector::iterator it = forces.begin(); + it != forces.end(); it += 3) + { + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + } + } + + delete mgmol; + + } // close main scope + + mgmol_finalize(); + + mpirc = MPI_Finalize(); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Finalize failed!!!" << std::endl; + } + + time_t tt; + time(&tt); + if (onpe0) std::cout << " Run ended at " << ctime(&tt) << std::endl; + + return 0; +} From e0ad1a9604ccce6baaa57b4b5d6a34c1986f4d15 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Thu, 8 Aug 2024 08:17:20 -0400 Subject: [PATCH 17/17] Add new functionality to compute energy and forces (#267) * use specified initial conditions for wavefunctions --- src/Control.cc | 36 --- src/Control.h | 1 - src/IonicAlgorithm.cc | 2 +- src/LBFGS.cc | 2 +- src/MGmol.cc | 50 +++-- src/MGmol.h | 18 +- src/MGmolInterface.h | 7 + src/Orbitals.h | 2 + src/md.cc | 20 +- src/quench.cc | 24 +- tests/CMakeLists.txt | 10 + tests/WFEnergyAndForces/mgmol.cfg | 28 +++ tests/WFEnergyAndForces/sih4.xyz | 8 + tests/WFEnergyAndForces/test.py | 91 ++++++++ .../testWFEnergyAndForces.cc | 210 ++++++++++++++++++ 15 files changed, 431 insertions(+), 78 deletions(-) create mode 100644 tests/WFEnergyAndForces/mgmol.cfg create mode 100644 tests/WFEnergyAndForces/sih4.xyz create mode 100755 tests/WFEnergyAndForces/test.py create mode 100644 tests/WFEnergyAndForces/testWFEnergyAndForces.cc diff --git a/src/Control.cc b/src/Control.cc index 7d8dc33f..af3d8ed9 100644 --- a/src/Control.cc +++ b/src/Control.cc @@ -997,42 +997,6 @@ void Control::readRestartInfo(std::ifstream* tfile) } } -void Control::readRestartOutputInfo(std::ifstream* tfile) -{ - const std::string zero = "0"; - const std::string one = "1"; - if (tfile != nullptr) - { - // Read in the output restart filename - std::string filename; - (*tfile) >> filename; - // int dpcs_chkpoint=0; - if (zero.compare(filename) == 0) // no restart dump - out_restart_info = 0; - else - { - if (one.compare(filename) == 0) - { // automatic naming of dump - filename = "snapshot"; - out_restart_file_naming_strategy = 1; - } - (*tfile) >> out_restart_info; - (*tfile) >> out_restart_file_type; - //(*tfile)>>dpcs_chkpoint; - // timeout_.set(dpcs_chkpoint); - } - - out_restart_file.assign(run_directory_); - out_restart_file.append("/"); - out_restart_file.append(filename); - (*MPIdata::sout) << "Output restart file: " << out_restart_file - << " with info level " << out_restart_info - << std::endl; - //(*MPIdata::sout)<<"Time for DPCS checkpoint: - //"<::init(HDFrestart* h5f_file) template int IonicAlgorithm::quenchElectrons(const int itmax, double& etot) { - int ret = mgmol_strategy_.quench(*orbitals_, ions_, itmax, 0, etot); + int ret = mgmol_strategy_.quench(**orbitals_, ions_, itmax, 0, etot); return ret; } diff --git a/src/LBFGS.cc b/src/LBFGS.cc index 0837a0e1..9154b7b3 100644 --- a/src/LBFGS.cc +++ b/src/LBFGS.cc @@ -97,7 +97,7 @@ int LBFGS::quenchElectrons(const int itmax, double& etot) { etot_i_[0] = etot_i_[1]; etot_i_[1] = etot_i_[2]; - int ret = mgmol_strategy_.quench(*orbitals_, ions_, itmax, 0, etot_i_[2]); + int ret = mgmol_strategy_.quench(**orbitals_, ions_, itmax, 0, etot_i_[2]); etot = etot_i_[2]; return ret; diff --git a/src/MGmol.cc b/src/MGmol.cc index 682a0e05..d5daabe6 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -521,7 +521,8 @@ void MGmol::run() switch (ct.AtomsDynamic()) { case AtomsDynamicType::Quench: - quench(current_orbitals_, *ions_, ct.max_electronic_steps, 20, eks); + quench( + *current_orbitals_, *ions_, ct.max_electronic_steps, 20, eks); // Forces for the last states force(*current_orbitals_, *ions_); @@ -1097,25 +1098,21 @@ void MGmol::setup() } template -void MGmol::cleanup() +void MGmol::dumpRestart() { - closing_tm_.start(); - - Mesh* mymesh = Mesh::instance(); - const pb::PEenv& myPEenv = mymesh->peenv(); - Control& ct = *(Control::instance()); - - printTimers(); + Control& ct = *(Control::instance()); // Save data to restart file - if (ct.out_restart_info > 0 && !ct.AtomsMove()) + if (ct.out_restart_info > 0) { + Mesh* mymesh = Mesh::instance(); const pb::Grid& mygrid = mymesh->grid(); unsigned gdim[3] = { mygrid.gdim(0), mygrid.gdim(1), mygrid.gdim(2) }; + const pb::PEenv& myPEenv = mymesh->peenv(); // create restart file std::string filename(std::string(ct.out_restart_file)); - filename += "0"; + if (ct.out_restart_file_naming_strategy) filename += "0"; HDFrestart h5restartfile( filename, myPEenv, gdim, ct.out_restart_file_type); @@ -1125,6 +1122,22 @@ void MGmol::cleanup() if (ierr < 0) os_ << "WARNING: writing restart data failed!!!" << std::endl; } +} + +template +void MGmol::cleanup() +{ + closing_tm_.start(); + + Control& ct = *(Control::instance()); + + printTimers(); + + // Save data to restart file + if (!ct.AtomsMove()) + { + dumpRestart(); + } MPI_Barrier(comm_); closing_tm_.stop(); @@ -1421,6 +1434,14 @@ template double MGmol::evaluateEnergyAndForces( const std::vector& tau, std::vector& atnumbers, std::vector& forces) +{ + return evaluateEnergyAndForces(current_orbitals_, tau, atnumbers, forces); +} + +template +double MGmol::evaluateEnergyAndForces(Orbitals* orbitals, + const std::vector& tau, std::vector& atnumbers, + std::vector& forces) { assert(tau.size() == 3 * atnumbers.size()); @@ -1430,10 +1451,11 @@ double MGmol::evaluateEnergyAndForces( moveVnuc(*ions_); - double eks = 0.; - quench(current_orbitals_, *ions_, ct.max_electronic_steps, 20, eks); + double eks = 0.; + OrbitalsType* dorbitals = dynamic_cast(orbitals); + quench(*dorbitals, *ions_, ct.max_electronic_steps, 20, eks); - force(*current_orbitals_, *ions_); + force(*dorbitals, *ions_); ions_->getForces(forces); diff --git a/src/MGmol.h b/src/MGmol.h index 8a3ae1bb..625c7b3f 100644 --- a/src/MGmol.h +++ b/src/MGmol.h @@ -42,7 +42,6 @@ class IonicAlgorithm; #include "AOMMprojector.h" #include "ClusterOrbitals.h" #include "DMStrategy.h" -#include "Energy.h" #include "ExtendedGridOrbitals.h" #include "Forces.h" #include "Ions.h" @@ -131,7 +130,7 @@ class MGmol : public MGmolInterface KBPsiMatrixSparse* kbpsi, dist_matrix::DistMatrix& hij); void computeHnlPhiAndAdd2HPhi(Ions& ions, OrbitalsType& phi, OrbitalsType& hphi, const KBPsiMatrixSparse* const kbpsi); - int dumprestartFile(OrbitalsType** orbitals, Ions& ions, + int dumpMDrestartFile(OrbitalsType** orbitals, Ions& ions, Rho& rho, const bool write_extrapolated_wf, const short count); @@ -190,6 +189,9 @@ class MGmol : public MGmolInterface double evaluateEnergyAndForces(const std::vector& tau, std::vector& atnumbers, std::vector& forces); + double evaluateEnergyAndForces(Orbitals*, const std::vector& tau, + std::vector& atnumbers, std::vector& forces); + /* * get internal atomic positions */ @@ -247,7 +249,7 @@ class MGmol : public MGmolInterface void update_pot(const pb::GridFunc& vh_init, const Ions& ions); void update_pot(const Ions& ions); - int quench(OrbitalsType* orbitals, Ions& ions, const int max_steps, + int quench(OrbitalsType& orbitals, Ions& ions, const int max_steps, const int iprint, double& last_eks); int outerSolve(OrbitalsType& orbitals, OrbitalsType& work_orbitals, Ions& ions, const int max_steps, const int iprint, double& last_eks); @@ -320,7 +322,17 @@ class MGmol : public MGmolInterface forces_->force(orbitals, ions); } + /* + * simply dump current state + */ + void dumpRestart(); + void loadRestartFile(const std::string filename); + + std::shared_ptr getProjectedMatrices() + { + return proj_matrices_; + } }; // Instantiate static variables here to avoid clang warnings template diff --git a/src/MGmolInterface.h b/src/MGmolInterface.h index 046459d6..af7be67b 100644 --- a/src/MGmolInterface.h +++ b/src/MGmolInterface.h @@ -28,8 +28,15 @@ class MGmolInterface virtual double evaluateEnergyAndForces(const std::vector& tau, std::vector& atnumbers, std::vector& forces) = 0; + virtual double evaluateEnergyAndForces(Orbitals*, + const std::vector& tau, std::vector& atnumbers, + std::vector& forces) + = 0; virtual void getAtomicPositions(std::vector& tau) = 0; virtual void getAtomicNumbers(std::vector& an) = 0; + virtual std::shared_ptr getProjectedMatrices() + = 0; + virtual void dumpRestart() = 0; }; #endif diff --git a/src/Orbitals.h b/src/Orbitals.h index d58f46bb..bc524aff 100644 --- a/src/Orbitals.h +++ b/src/Orbitals.h @@ -26,6 +26,8 @@ class Orbitals Orbitals() { iterative_index_ = -10; } + virtual ~Orbitals(){}; + Orbitals(const Orbitals& A, const bool copy_data) { if (copy_data) diff --git a/src/md.cc b/src/md.cc index f7fe934d..7c6e2aa8 100644 --- a/src/md.cc +++ b/src/md.cc @@ -226,7 +226,7 @@ void checkMaxForces(const std::vector& fion, } template -int MGmol::dumprestartFile(OrbitalsType** orbitals, Ions& ions, +int MGmol::dumpMDrestartFile(OrbitalsType** orbitals, Ions& ions, Rho& rho, const bool write_extrapolated_wf, const short count) { MGmol_MPI& mmpi(*(MGmol_MPI::instance())); @@ -255,7 +255,7 @@ int MGmol::dumprestartFile(OrbitalsType** orbitals, Ions& ions, { if (onpe0) (*MPIdata::serr) - << "dumprestartFile: cannot write ...previous_orbitals..." + << "dumpMDrestartFile: cannot write ...previous_orbitals..." << std::endl; return ierr; } @@ -269,7 +269,7 @@ int MGmol::dumprestartFile(OrbitalsType** orbitals, Ions& ions, if (ierr < 0) { if (onpe0) - (*MPIdata::serr) << "dumprestartFile: cannot write " + (*MPIdata::serr) << "dumpMDrestartFile: cannot write " "...ExtrapolatedFunction..." << std::endl; return ierr; @@ -282,7 +282,7 @@ int MGmol::dumprestartFile(OrbitalsType** orbitals, Ions& ions, { if (onpe0) (*MPIdata::serr) - << "dumprestartFile: cannot close file..." << std::endl; + << "dumpMDrestartFile: cannot close file..." << std::endl; return ierr; } @@ -398,7 +398,7 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) if (ct.restart_info < 3) { double eks = 0.; - quench(*orbitals, ions, ct.max_electronic_steps, 20, eks); + quench(**orbitals, ions, ct.max_electronic_steps, 20, eks); } ct.max_changes_pot = 0; @@ -426,7 +426,7 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) bool last_move_is_small = true; do { - retval = quench(*orbitals, ions, ct.max_electronic_steps, 0, eks); + retval = quench(**orbitals, ions, ct.max_electronic_steps, 0, eks); // update localization regions if (ct.adaptiveLRs()) @@ -614,12 +614,12 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) while (ierr < 0 && count < DUMP_MAX_NUM_TRY) { dump_tm_.start(); - ierr = dumprestartFile( + ierr = dumpMDrestartFile( orbitals, ions, *rho_, extrapolated_flag, count); dump_tm_.stop(); if (onpe0 && ierr < 0 && count < (DUMP_MAX_NUM_TRY - 1)) std::cout - << "dumprestartFile() failed... try again..." + << "dumpMDrestartFile() failed... try again..." << std::endl; if (ierr < 0) sleep(1.); count++; @@ -640,12 +640,12 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) while (ierr < 0 && count < DUMP_MAX_NUM_TRY) { dump_tm_.start(); - ierr = dumprestartFile( + ierr = dumpMDrestartFile( orbitals, ions, *rho_, extrapolated_flag, count); dump_tm_.stop(); if (onpe0 && ierr < 0 && count < (DUMP_MAX_NUM_TRY - 1)) - std::cout << "dumprestartFile() failed... try again..." + std::cout << "dumpMDrestartFile() failed... try again..." << std::endl; if (ierr < 0) sleep(1.); count++; diff --git a/src/quench.cc b/src/quench.cc index 614aa636..59a16d50 100644 --- a/src/quench.cc +++ b/src/quench.cc @@ -519,7 +519,7 @@ int MGmol::outerSolve(OrbitalsType& orbitals, } template -int MGmol::quench(OrbitalsType* orbitals, Ions& ions, +int MGmol::quench(OrbitalsType& orbitals, Ions& ions, const int max_inner_steps, const int iprint, double& last_eks) { assert(max_inner_steps > -1); @@ -543,33 +543,33 @@ int MGmol::quench(OrbitalsType* orbitals, Ions& ions, } // get actual indexes of stored functions - const std::vector>& gids(orbitals->getOverlappingGids()); + const std::vector>& gids(orbitals.getOverlappingGids()); g_kbpsi_->setup(*ions_); electrostat_->setup(ct.vh_its); rho_->setup(ct.getOrthoType(), gids); - OrbitalsType work_orbitals("Work", *orbitals); + OrbitalsType work_orbitals("Work", orbitals); - orbitals->setDataWithGhosts(); - orbitals->trade_boundaries(); + orbitals.setDataWithGhosts(); + orbitals.trade_boundaries(); - disentangleOrbitals(*orbitals, work_orbitals, ions, max_steps); + disentangleOrbitals(orbitals, work_orbitals, ions, max_steps); // setup "kernel" functions for AOMM algorithm if (ct.use_kernel_functions) { - applyAOMMprojection(*orbitals); + applyAOMMprojection(orbitals); } orbitals_precond_.reset(new OrbitalsPreconditioning()); orbitals_precond_->setup( - *orbitals, ct.getMGlevels(), ct.lap_type, currentMasks_.get(), lrs_); + orbitals, ct.getMGlevels(), ct.lap_type, currentMasks_.get(), lrs_); // solve electronic structure problem // (inner iterations) retval = outerSolve( - *orbitals, work_orbitals, ions, max_steps, iprint, last_eks); + orbitals, work_orbitals, ions, max_steps, iprint, last_eks); if (retval == -1) return -1; if (ct.use_kernel_functions) @@ -602,7 +602,7 @@ int MGmol::quench(OrbitalsType* orbitals, Ions& ions, } } last_eks - = energy_->evaluateTotal(ts, proj_matrices_.get(), *orbitals, 2, os_); + = energy_->evaluateTotal(ts, proj_matrices_.get(), orbitals, 2, os_); if (ct.computeCondGramMD()) { @@ -615,7 +615,7 @@ int MGmol::quench(OrbitalsType* orbitals, Ions& ions, if (ct.isLocMode() || ct.isSpreadFunctionalActive()) { // build matrices necessary to compute spreads and centers - spreadf_->computePositionMatrix(*orbitals, work_orbitals); + spreadf_->computePositionMatrix(orbitals, work_orbitals); if (ct.verbose > 0 || !ct.AtomsMove()) { @@ -628,7 +628,7 @@ int MGmol::quench(OrbitalsType* orbitals, Ions& ions, { if (ct.wannier_transform_type >= 1) { - wftransform(orbitals, &work_orbitals, ions); + wftransform(&orbitals, &work_orbitals, ions); } } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 3f221912..11125081 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -248,6 +248,8 @@ add_executable(testDensityMatrix ${CMAKE_SOURCE_DIR}/tests/ut_magma_main.cc) add_executable(testEnergyAndForces ${CMAKE_SOURCE_DIR}/tests/EnergyAndForces/testEnergyAndForces.cc) +add_executable(testWFEnergyAndForces + ${CMAKE_SOURCE_DIR}/tests/WFEnergyAndForces/testWFEnergyAndForces.cc) if(${MAGMA_FOUND}) add_executable(testOpenmpOffload @@ -344,6 +346,13 @@ add_test(NAME testEnergyAndForces ${CMAKE_CURRENT_SOURCE_DIR}/EnergyAndForces/coords.in ${CMAKE_CURRENT_SOURCE_DIR}/EnergyAndForces/lrs.in ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) +add_test(NAME testWFEnergyAndForces + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/WFEnergyAndForces/test.py + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/testWFEnergyAndForces + ${CMAKE_CURRENT_SOURCE_DIR}/WFEnergyAndForces/mgmol.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/WFEnergyAndForces/sih4.xyz + ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) if(${MAGMA_FOUND}) add_test(NAME testOpenmpOffload @@ -529,6 +538,7 @@ target_link_libraries(testBlacsContext PRIVATE ${SCALAPACK_LIBRARIES} target_link_libraries(testSuperSampling PRIVATE MPI::MPI_CXX) target_link_libraries(testDirectionalReduce PRIVATE MPI::MPI_CXX) target_link_libraries(testEnergyAndForces PRIVATE mgmol_src) +target_link_libraries(testWFEnergyAndForces PRIVATE mgmol_src) if(${MAGMA_FOUND}) target_link_libraries(testDistVector PRIVATE ${SCALAPACK_LIBRARIES} diff --git a/tests/WFEnergyAndForces/mgmol.cfg b/tests/WFEnergyAndForces/mgmol.cfg new file mode 100644 index 00000000..d543b626 --- /dev/null +++ b/tests/WFEnergyAndForces/mgmol.cfg @@ -0,0 +1,28 @@ +verbosity=2 +xcFunctional=LDA +[Mesh] +nx=40 +ny=40 +nz=40 +[Domain] +ox=-6.75 +oy=-6.75 +oz=-6.75 +lx=13.5 +ly=13.5 +lz=13.5 +[Potentials] +pseudopotential=pseudo.Si +pseudopotential=pseudo.H +[Run] +type=QUENCH +[Quench] +max_steps=50 +atol=1.e-9 +num_lin_iterations=2 +[Orbitals] +initial_type=Gaussian +initial_width=2. +[Restart] +output_level=3 +output_filename=WF diff --git a/tests/WFEnergyAndForces/sih4.xyz b/tests/WFEnergyAndForces/sih4.xyz new file mode 100644 index 00000000..b3f921e3 --- /dev/null +++ b/tests/WFEnergyAndForces/sih4.xyz @@ -0,0 +1,8 @@ +5 +SiH4 molecule (coordinates in Angstrom) +Si 0.0 0.0 0.0 +H 0.885 0.885 0.885 +H -0.885 -0.885 0.885 +H -0.885 0.885 -0.885 +H 0.885 -0.885 -0.885 + diff --git a/tests/WFEnergyAndForces/test.py b/tests/WFEnergyAndForces/test.py new file mode 100755 index 00000000..45420ddf --- /dev/null +++ b/tests/WFEnergyAndForces/test.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python +import sys +import os +import subprocess +import string + +print("Test WFEnergyAndForces...") + +nargs=len(sys.argv) + +mpicmd = sys.argv[1]+" "+sys.argv[2]+" "+sys.argv[3] +for i in range(4,nargs-4): + mpicmd = mpicmd + " "+sys.argv[i] +print("MPI run command: {}".format(mpicmd)) + +exe = sys.argv[nargs-4] +inp = sys.argv[nargs-3] +coords = sys.argv[nargs-2] +print("coordinates file: %s"%coords) + +#create links to potentials files +dst1 = 'pseudo.Si' +dst2 = 'pseudo.H' +src1 = sys.argv[nargs-1] + '/' + dst1 +src2 = sys.argv[nargs-1] + '/' + dst2 + +if not os.path.exists(dst1): + print("Create link to %s"%dst1) + os.symlink(src1, dst1) +if not os.path.exists(dst2): + print("Create link to %s"%dst2) + os.symlink(src2, dst2) + +#run +command = "{} {} -c {} -i {}".format(mpicmd,exe,inp,coords) +print("Run command: {}".format(command)) +output = subprocess.check_output(command,shell=True) +lines=output.split(b'\n') + +#analyse output +energies=[] +for line in lines: + if line.count(b'%%'): + print(line) + words=line.split() + words=words[5].split(b',')[0] + energy = words.decode() + if line.count(b'achieved'): + energies.append(energy) + +flag=0 +forces=[] +for line in lines: + if flag>0: + print(line) + words=line.split(b' ') + forces.append(words[1].decode()) + forces.append(words[2].decode()) + forces.append(words[3].decode()) + flag=flag-1 + if line.count(b'Forces:'): + flag=2 + + +print("Check energies...") +print( energies ) +if len(energies)<2: + print("Expected two converged energies") + sys.exit(1) + +tol = 1.e-6 +diff=eval(energies[1])-eval(energies[0]) +print(diff) +if abs(diff)>tol: + print("Energies differ: {} vs {} !!!".format(energies[0],energies[1])) + sys.exit(1) + +print("Check forces...") +print(forces) +flag=0 +for i in range(6): + diff=eval(forces[i+6])-eval(forces[i]) + print(diff) + if abs(diff)>1.e-3: + print("Forces difference larger than tol") + flag=1 +if flag>0: + sys.exit(1) + +print("Test SUCCESSFUL!") +sys.exit(0) diff --git a/tests/WFEnergyAndForces/testWFEnergyAndForces.cc b/tests/WFEnergyAndForces/testWFEnergyAndForces.cc new file mode 100644 index 00000000..c9e28d48 --- /dev/null +++ b/tests/WFEnergyAndForces/testWFEnergyAndForces.cc @@ -0,0 +1,210 @@ +// Copyright (c) 2017, Lawrence Livermore National Security, LLC and +// UT-Battelle, LLC. +// Produced at the Lawrence Livermore National Laboratory and the Oak Ridge +// National Laboratory. +// LLNL-CODE-743438 +// All rights reserved. +// 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 "Control.h" +#include "ExtendedGridOrbitals.h" +#include "LocGridOrbitals.h" +#include "MGmol.h" +#include "MGmol_MPI.h" +#include "MPIdata.h" +#include "mgmol_run.h" + +#include +#include +#include +#include + +#include +namespace po = boost::program_options; + +int main(int argc, char** argv) +{ + int mpirc = MPI_Init(&argc, &argv); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Initialization failed!!!" << std::endl; + MPI_Abort(MPI_COMM_WORLD, 0); + } + + MPI_Comm comm = MPI_COMM_WORLD; + + /* + * Initialize general things, like magma, openmp, IO, ... + */ + mgmol_init(comm); + + /* + * read runtime parameters + */ + std::string input_filename(""); + std::string lrs_filename; + std::string constraints_filename(""); + + float total_spin = 0.; + bool with_spin = false; + + po::variables_map vm; + + // read from PE0 only + if (MPIdata::onpe0) + { + read_config(argc, argv, vm, input_filename, lrs_filename, + constraints_filename, total_spin, with_spin); + } + + MGmol_MPI::setup(comm, std::cout, with_spin); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + MPI_Comm global_comm = mmpi.commGlobal(); + + /* + * Setup control struct with run time parameters + */ + Control::setup(global_comm, with_spin, total_spin); + Control& ct = *(Control::instance()); + + ct.setOptions(vm); + + int ret = ct.checkOptions(); + if (ret < 0) return ret; + + mmpi.bcastGlobal(input_filename); + mmpi.bcastGlobal(lrs_filename); + + // Enter main scope + { + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "Construct MGmol object..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + + MGmolInterface* mgmol; + if (ct.isLocMode()) + 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); + + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "MGmol setup..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + mgmol->setup(); + + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "Setup done..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + + // here we just use the atomic positions read in and used + // to initialize MGmol + std::vector positions; + mgmol->getAtomicPositions(positions); + std::vector anumbers; + mgmol->getAtomicNumbers(anumbers); + if (MPIdata::onpe0) + { + std::cout << "Positions:" << std::endl; + std::vector::iterator ita = anumbers.begin(); + for (std::vector::iterator it = positions.begin(); + it != positions.end(); it += 3) + { + std::cout << *ita; + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + ita++; + } + } + + // compute energy and forces using all MPI tasks + // expect positions to be replicated on all MPI tasks + std::vector forces; + double eks + = mgmol->evaluateEnergyAndForces(positions, anumbers, forces); + mgmol->dumpRestart(); + + // print out results + if (MPIdata::onpe0) + { + std::cout << "Eks: " << eks << std::endl; + std::cout << "Forces:" << std::endl; + for (std::vector::iterator it = forces.begin(); + it != forces.end(); it += 3) + { + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + } + } + + // compute energy and forces again using wavefunctions + // from previous call + Mesh* mymesh = Mesh::instance(); + const pb::Grid& mygrid = mymesh->grid(); + + std::shared_ptr projmatrices + = mgmol->getProjectedMatrices(); + + ExtendedGridOrbitals orbitals("new_orbitals", mygrid, mymesh->subdivx(), + ct.numst, ct.bcWF, projmatrices.get(), nullptr, nullptr, nullptr, + nullptr); + + const pb::PEenv& myPEenv = mymesh->peenv(); + HDFrestart h5file("WF", myPEenv, ct.out_restart_file_type); + orbitals.read_hdf5(h5file); + + // + // evaluate energy and forces again + // + + // convergence should be really quick since we start with an initial + // guess which is the solution + ct.max_electronic_steps = 5; + eks = mgmol->evaluateEnergyAndForces( + &orbitals, positions, anumbers, forces); + + // print out results + if (MPIdata::onpe0) + { + std::cout << "Eks: " << eks << std::endl; + std::cout << "Forces:" << std::endl; + for (std::vector::iterator it = forces.begin(); + it != forces.end(); it += 3) + { + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + } + } + + delete mgmol; + + } // close main scope + + mgmol_finalize(); + + mpirc = MPI_Finalize(); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Finalize failed!!!" << std::endl; + } + + time_t tt; + time(&tt); + if (onpe0) std::cout << " Run ended at " << ctime(&tt) << std::endl; + + return 0; +}