From e0ad1a9604ccce6baaa57b4b5d6a34c1986f4d15 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Thu, 8 Aug 2024 08:17:20 -0400 Subject: [PATCH] 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; +}