diff --git a/src/DensityMatrix.cc b/src/DensityMatrix.cc index 87b30224..6e6ea652 100644 --- a/src/DensityMatrix.cc +++ b/src/DensityMatrix.cc @@ -149,14 +149,11 @@ template void DensityMatrix::setUniform( const double nel, const int new_orbitals_index) { -#ifdef PRINT_OPERATIONS - MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + const double occ = (double)((double)nel / (double)dim_); if (mmpi.instancePE0()) - std::cout << "template " - "DensityMatrix::setUniform()" + std::cout << "DensityMatrix::setUniform(), occupation = " << occ << std::endl; -#endif - const double occ = (double)((double)nel / (double)dim_); assert(occ < 1.01); for (int i = 0; i < dim_; i++) occupation_[i] = occ; diff --git a/src/MGmol.cc b/src/MGmol.cc index d5daabe6..be4f28aa 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -1432,7 +1432,7 @@ void MGmol::getAtomicNumbers(std::vector& an) template double MGmol::evaluateEnergyAndForces( - const std::vector& tau, std::vector& atnumbers, + const std::vector& tau, const std::vector& atnumbers, std::vector& forces) { return evaluateEnergyAndForces(current_orbitals_, tau, atnumbers, forces); @@ -1440,7 +1440,7 @@ double MGmol::evaluateEnergyAndForces( template double MGmol::evaluateEnergyAndForces(Orbitals* orbitals, - const std::vector& tau, std::vector& atnumbers, + const std::vector& tau, const std::vector& atnumbers, std::vector& forces) { assert(tau.size() == 3 * atnumbers.size()); @@ -1462,6 +1462,46 @@ double MGmol::evaluateEnergyAndForces(Orbitals* orbitals, return eks; } +template +double MGmol::evaluateDMandEnergyAndForces(Orbitals* orbitals, + const std::vector& tau, const std::vector& atnumbers, + std::vector& forces) +{ + OrbitalsType* dorbitals = dynamic_cast(orbitals); + + ions_->setPositions(tau, atnumbers); + + moveVnuc(*ions_); + + // initialize electronic density + rho_->update(*dorbitals); + + // initialize potential + update_pot(*ions_); + + // initialize projected matrices + updateHmatrix(*dorbitals, *ions_); + proj_matrices_->updateThetaAndHB(); + + // compute DM + std::shared_ptr> dm_strategy( + DMStrategyFactory>::create(comm_, os_, *ions_, + rho_.get(), energy_.get(), electrostat_.get(), this, + proj_matrices_.get(), dorbitals)); + + dm_strategy->update(*dorbitals); + + // evaluate energy and forces + double ts = 0.; + double eks + = energy_->evaluateTotal(ts, proj_matrices_.get(), *dorbitals, 2, os_); + + force(*dorbitals, *ions_); + + return eks; +} + template class MGmol; template class MGmol; template int MGmol::initial(); diff --git a/src/MGmol.h b/src/MGmol.h index 625c7b3f..aee2b660 100644 --- a/src/MGmol.h +++ b/src/MGmol.h @@ -186,11 +186,30 @@ class MGmol : public MGmolInterface void run() override; + /* + * Evaluate the energy and forces for an atomic configuration + * specified by tau (input) + */ double evaluateEnergyAndForces(const std::vector& tau, - std::vector& atnumbers, std::vector& forces); + const std::vector& atnumbers, std::vector& forces); - double evaluateEnergyAndForces(Orbitals*, const std::vector& tau, - std::vector& atnumbers, std::vector& forces); + /* + * Evaluate the energy and forces for an atomic configuration + * specified by tau (input), using the input "orbitals" as initial + * guess for the wavefunctions + */ + double evaluateEnergyAndForces(Orbitals* orbitals, + const std::vector& tau, const std::vector& atnumbers, + std::vector& forces); + + /* + * Evaluate the energy and forces for an atomic configuration + * specified by tau (input), using the input "orbitals" as wavefunctions + * (fixed) + */ + double evaluateDMandEnergyAndForces(Orbitals* orbitals, + const std::vector& tau, const std::vector& atnumbers, + std::vector& forces); /* * get internal atomic positions diff --git a/src/MGmolInterface.h b/src/MGmolInterface.h index af7be67b..9a9bf8a6 100644 --- a/src/MGmolInterface.h +++ b/src/MGmolInterface.h @@ -25,13 +25,19 @@ 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) + const std::vector& atnumbers, std::vector& forces) = 0; virtual double evaluateEnergyAndForces(Orbitals*, - const std::vector& tau, std::vector& atnumbers, + const std::vector& tau, const std::vector& atnumbers, + std::vector& forces) + = 0; + virtual double evaluateDMandEnergyAndForces(Orbitals*, + const std::vector& tau, const 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() diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc index 7feb564d..6115f933 100644 --- a/src/MVPSolver.cc +++ b/src/MVPSolver.cc @@ -208,6 +208,7 @@ int MVPSolver::solve(OrbitalsType& orbitals) kbpsi.computeHvnlMatrix(&kbpsi, ions_, h11_nl); + const double tol_de0 = 1.e-12; for (int inner_it = 0; inner_it < n_inner_steps_; inner_it++) { if (onpe0 && ct.verbose > 1) @@ -270,6 +271,14 @@ int MVPSolver::solve(OrbitalsType& orbitals) double de0 = evaluateDerivative(dmInit, delta_dm, ts0); + if (std::abs(de0) < tol_de0 && inner_it > 0) + { + if (onpe0 && ct.verbose > 0) + std::cout << "MVP: de0 = " << de0 + << ", convergence achieved" << std::endl; + break; + } + // // evaluate free energy at beta=1 // @@ -306,6 +315,7 @@ int MVPSolver::solve(OrbitalsType& orbitals) // line minimization const double beta = minQuadPolynomial(e0, e1, de0, (ct.verbose > 2), os_); + assert(!std::isnan(beta)); if (onpe0 && ct.verbose > 0) { @@ -348,6 +358,8 @@ int MVPSolver::solve(OrbitalsType& orbitals) orbitals.getProjMatrices()); projmatrices->setDM(*work_, orbitals.getIterativeIndex()); projmatrices->setEigenvalues(proj_mat_work_->getEigenvalues()); + projmatrices->assignH(proj_mat_work_->getH()); + projmatrices->setHB2H(); } // Generate new density diff --git a/src/ProjectedMatrices.cc b/src/ProjectedMatrices.cc index bdfe7786..e0e679e9 100644 --- a/src/ProjectedMatrices.cc +++ b/src/ProjectedMatrices.cc @@ -650,8 +650,9 @@ double ProjectedMatrices::computeEntropy() else { if (mmpi.PE0() && ct.verbose > 1) - (*MPIdata::sout) - << "occupations uptodate, skip computation..." << std::endl; + (*MPIdata::sout) << "computeEntropy: occupations uptodate, " + "skip computation..." + << std::endl; } entropy = computeEntropy(width_); } diff --git a/src/Rho.cc b/src/Rho.cc index ca0f20ea..45a41f3e 100644 --- a/src/Rho.cc +++ b/src/Rho.cc @@ -90,8 +90,8 @@ void Rho::update(OrbitalsType& current_orbitals) update_tm_.start(); - if (verbosity_level_ > 2 && onpe0) - (*MPIdata::sout) << "Rho::update()" << std::endl; + if (verbosity_level_ > 1 && onpe0) + (*MPIdata::sout) << "Rho::update()..." << std::endl; const int new_iterative_index = ((1 + current_orbitals.getIterativeIndex()) % 100) @@ -99,7 +99,7 @@ void Rho::update(OrbitalsType& current_orbitals) if (iterative_index_ == new_iterative_index) { - if (onpe0 && verbosity_level_ > 2) + if (onpe0 && verbosity_level_ > 1) (*MPIdata::sout) << "Rho already up to date, iterative_index_=" << iterative_index_ << std::endl; return; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 11125081..cc205167 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -250,6 +250,8 @@ add_executable(testEnergyAndForces ${CMAKE_SOURCE_DIR}/tests/EnergyAndForces/testEnergyAndForces.cc) add_executable(testWFEnergyAndForces ${CMAKE_SOURCE_DIR}/tests/WFEnergyAndForces/testWFEnergyAndForces.cc) +add_executable(testDMandEnergyAndForces + ${CMAKE_SOURCE_DIR}/tests/DMandEnergyAndForces/testDMandEnergyAndForces.cc) if(${MAGMA_FOUND}) add_executable(testOpenmpOffload @@ -353,6 +355,14 @@ add_test(NAME testWFEnergyAndForces ${CMAKE_CURRENT_SOURCE_DIR}/WFEnergyAndForces/mgmol.cfg ${CMAKE_CURRENT_SOURCE_DIR}/WFEnergyAndForces/sih4.xyz ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) +add_test(NAME testDMandEnergyAndForces + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/DMandEnergyAndForces/test.py + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/testDMandEnergyAndForces + ${CMAKE_CURRENT_SOURCE_DIR}/DMandEnergyAndForces/mgmol.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/DMandEnergyAndForces/coords.in + ${CMAKE_CURRENT_SOURCE_DIR}/DMandEnergyAndForces/lrs.in + ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) if(${MAGMA_FOUND}) add_test(NAME testOpenmpOffload @@ -539,6 +549,7 @@ 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) +target_link_libraries(testDMandEnergyAndForces PRIVATE mgmol_src) if(${MAGMA_FOUND}) target_link_libraries(testDistVector PRIVATE ${SCALAPACK_LIBRARIES} diff --git a/tests/DMandEnergyAndForces/coords.in b/tests/DMandEnergyAndForces/coords.in new file mode 100644 index 00000000..cda898b9 --- /dev/null +++ b/tests/DMandEnergyAndForces/coords.in @@ -0,0 +1,2 @@ +N1 1 0. 0. -1.0345 +N2 1 0. 0. 1.0345 diff --git a/tests/DMandEnergyAndForces/lrs.in b/tests/DMandEnergyAndForces/lrs.in new file mode 100644 index 00000000..9136a500 --- /dev/null +++ b/tests/DMandEnergyAndForces/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/DMandEnergyAndForces/mgmol.cfg b/tests/DMandEnergyAndForces/mgmol.cfg new file mode 100644 index 00000000..5278f40b --- /dev/null +++ b/tests/DMandEnergyAndForces/mgmol.cfg @@ -0,0 +1,35 @@ +verbosity=2 +xcFunctional=LDA +FDtype=4th +[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-8 +step_length=2. +ortho_freq=10 +[Orbitals] +initial_type=Gaussian +initial_width=1.5 +temperature=10. +nempty=2 +[Restart] +output_level=3 +output_filename=WF +[DensityMatrix] +solver=MVP +nb_inner_it=1 diff --git a/tests/DMandEnergyAndForces/test.py b/tests/DMandEnergyAndForces/test.py new file mode 100755 index 00000000..c36f0767 --- /dev/null +++ b/tests/DMandEnergyAndForces/test.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python +import sys +import os +import subprocess +import string + +print("Test DMandEnergyAndForces...") + +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) + break + +for line in lines: + if line.count(b'MVP') and line.count(b'iteration'): + print(line) + if line.count(b'Eks2'): + print(line) + words=line.split() + word=words[2] + energy = word.decode() + energies.append(energy) + break + +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("Test SUCCESSFUL!") +sys.exit(0) diff --git a/tests/DMandEnergyAndForces/testDMandEnergyAndForces.cc b/tests/DMandEnergyAndForces/testDMandEnergyAndForces.cc new file mode 100644 index 00000000..1cf1cf88 --- /dev/null +++ b/tests/DMandEnergyAndForces/testDMandEnergyAndForces.cc @@ -0,0 +1,211 @@ +// 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 << "Eks1 : " << eks << std::endl; + std::cout << "Forces2 :" << 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, with wavefunctions + // frozen to solution of previous problem + // + + // reset initial DM to test iterative solve for it + projmatrices->setDMuniform(ct.getNelSpin(), 0); + ct.dm_inner_steps = 50; + eks = mgmol->evaluateDMandEnergyAndForces( + &orbitals, positions, anumbers, forces); + + // print out results + if (MPIdata::onpe0) + { + std::cout << "Eks2 : " << eks << std::endl; + std::cout << "Forces2 :" << 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; +}