Skip to content

Add new functionality to compute energy and forces #267

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jeanlucf22 merged 1 commit into from
Aug 8, 2024
Merged

Add new functionality to compute energy and forces #267

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 0 additions & 36 deletions src/Control.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
//"<<dpcs_chkpoint<<"[s]"<<endl;
}
}

int Control::setPreconditionerParameters(const short type, const float factor,
const bool project_out, const short nlevels, const float fgrid_hmax)
{
Expand Down
1 change: 0 additions & 1 deletion src/Control.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -284,7 +284,6 @@ class Control
void adjust();
int checkState();
void readRestartInfo(std::ifstream* tfile);
void readRestartOutputInfo(std::ifstream* tfile);
int readThermostatInfo(std::ifstream* tfile);
void printThermostatInfo(std::ostream& os) const;
void setNumst(const short myspin, const int nval);
Expand Down
2 changes: 1 addition & 1 deletion src/IonicAlgorithm.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ void IonicAlgorithm<OrbitalsType>::init(HDFrestart* h5f_file)
template <class OrbitalsType>
int IonicAlgorithm<OrbitalsType>::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;
}
Expand Down
2 changes: 1 addition & 1 deletion src/LBFGS.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ int LBFGS<OrbitalsType>::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;
Expand Down
50 changes: 36 additions & 14 deletions src/MGmol.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -521,7 +521,8 @@ void MGmol<OrbitalsType>::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_);
Expand Down Expand Up @@ -1097,25 +1098,21 @@ void MGmol<OrbitalsType>::setup()
}

template <class OrbitalsType>
void MGmol<OrbitalsType>::cleanup()
void MGmol<OrbitalsType>::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);

Expand All @@ -1125,6 +1122,22 @@ void MGmol<OrbitalsType>::cleanup()
if (ierr < 0)
os_ << "WARNING: writing restart data failed!!!" << std::endl;
}
}

template <class OrbitalsType>
void MGmol<OrbitalsType>::cleanup()
{
closing_tm_.start();

Control& ct = *(Control::instance());

printTimers();

// Save data to restart file
if (!ct.AtomsMove())
{
dumpRestart();
}

MPI_Barrier(comm_);
closing_tm_.stop();
Expand Down Expand Up @@ -1421,6 +1434,14 @@ template <class OrbitalsType>
double MGmol<OrbitalsType>::evaluateEnergyAndForces(
const std::vector<double>& tau, std::vector<short>& atnumbers,
std::vector<double>& forces)
{
return evaluateEnergyAndForces(current_orbitals_, tau, atnumbers, forces);
}

template <class OrbitalsType>
double MGmol<OrbitalsType>::evaluateEnergyAndForces(Orbitals* orbitals,
const std::vector<double>& tau, std::vector<short>& atnumbers,
std::vector<double>& forces)
{
assert(tau.size() == 3 * atnumbers.size());

Expand All @@ -1430,10 +1451,11 @@ double MGmol<OrbitalsType>::evaluateEnergyAndForces(

moveVnuc(*ions_);

double eks = 0.;
quench(current_orbitals_, *ions_, ct.max_electronic_steps, 20, eks);
double eks = 0.;
OrbitalsType* dorbitals = dynamic_cast<OrbitalsType*>(orbitals);
quench(*dorbitals, *ions_, ct.max_electronic_steps, 20, eks);

force(*current_orbitals_, *ions_);
force(*dorbitals, *ions_);

ions_->getForces(forces);

Expand Down
18 changes: 15 additions & 3 deletions src/MGmol.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -131,7 +130,7 @@ class MGmol : public MGmolInterface
KBPsiMatrixSparse* kbpsi, dist_matrix::DistMatrix<DISTMATDTYPE>& 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<OrbitalsType>& rho, const bool write_extrapolated_wf,
const short count);

Expand Down Expand Up @@ -190,6 +189,9 @@ class MGmol : public MGmolInterface
double evaluateEnergyAndForces(const std::vector<double>& tau,
std::vector<short>& atnumbers, std::vector<double>& forces);

double evaluateEnergyAndForces(Orbitals*, const std::vector<double>& tau,
std::vector<short>& atnumbers, std::vector<double>& forces);

/*
* get internal atomic positions
*/
Expand Down Expand Up @@ -247,7 +249,7 @@ class MGmol : public MGmolInterface

void update_pot(const pb::GridFunc<POTDTYPE>& 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);
Expand Down Expand Up @@ -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<ProjectedMatricesInterface> getProjectedMatrices()
{
return proj_matrices_;
}
};
// Instantiate static variables here to avoid clang warnings
template <class OrbitalsType>
Expand Down
7 changes: 7 additions & 0 deletions src/MGmolInterface.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,15 @@ class MGmolInterface
virtual double evaluateEnergyAndForces(const std::vector<double>& tau,
std::vector<short>& atnumbers, std::vector<double>& forces)
= 0;
virtual double evaluateEnergyAndForces(Orbitals*,
const std::vector<double>& tau, std::vector<short>& atnumbers,
std::vector<double>& forces)
= 0;
virtual void getAtomicPositions(std::vector<double>& tau) = 0;
virtual void getAtomicNumbers(std::vector<short>& an) = 0;
virtual std::shared_ptr<ProjectedMatricesInterface> getProjectedMatrices()
= 0;
virtual void dumpRestart() = 0;
};

#endif
2 changes: 2 additions & 0 deletions src/Orbitals.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,8 @@ class Orbitals

Orbitals() { iterative_index_ = -10; }

virtual ~Orbitals(){};

Orbitals(const Orbitals& A, const bool copy_data)
{
if (copy_data)
Expand Down
20 changes: 10 additions & 10 deletions src/md.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -226,7 +226,7 @@ void checkMaxForces(const std::vector<double>& fion,
}

template <class OrbitalsType>
int MGmol<OrbitalsType>::dumprestartFile(OrbitalsType** orbitals, Ions& ions,
int MGmol<OrbitalsType>::dumpMDrestartFile(OrbitalsType** orbitals, Ions& ions,
Rho<OrbitalsType>& rho, const bool write_extrapolated_wf, const short count)
{
MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
Expand Down Expand Up @@ -255,7 +255,7 @@ int MGmol<OrbitalsType>::dumprestartFile(OrbitalsType** orbitals, Ions& ions,
{
if (onpe0)
(*MPIdata::serr)
<< "dumprestartFile: cannot write ...previous_orbitals..."
<< "dumpMDrestartFile: cannot write ...previous_orbitals..."
<< std::endl;
return ierr;
}
Expand All @@ -269,7 +269,7 @@ int MGmol<OrbitalsType>::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;
Expand All @@ -282,7 +282,7 @@ int MGmol<OrbitalsType>::dumprestartFile(OrbitalsType** orbitals, Ions& ions,
{
if (onpe0)
(*MPIdata::serr)
<< "dumprestartFile: cannot close file..." << std::endl;
<< "dumpMDrestartFile: cannot close file..." << std::endl;
return ierr;
}

Expand Down Expand Up @@ -398,7 +398,7 @@ void MGmol<OrbitalsType>::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;
Expand Down Expand Up @@ -426,7 +426,7 @@ void MGmol<OrbitalsType>::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())
Expand Down Expand Up @@ -614,12 +614,12 @@ void MGmol<OrbitalsType>::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++;
Expand All @@ -640,12 +640,12 @@ void MGmol<OrbitalsType>::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++;
Expand Down
Loading