Skip to content

Commit

Permalink
Merge branch 'ROMFPMD' into test_ROM_force
Browse files Browse the repository at this point in the history
  • Loading branch information
@siuwuncheung
siuwuncheung authored Oct 18, 2024
2 parents a27efba + b32d686 commit f0319f5
Show file tree
Hide file tree
Showing 32 changed files with 704 additions and 636 deletions.
1 change: 1 addition & 0 deletions src/Control.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,7 @@ Control::Control()
dm_approx_ndigits = 1;
dm_approx_power_maxits = 100;
wf_extrapolation_ = 1;
verbose = 0;

// undefined values
dm_algo_ = -1;
Expand Down
7 changes: 6 additions & 1 deletion src/DensityMatrix.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ const double factor_kernel4dot = 10.;
template <class MatrixType>
DensityMatrix<MatrixType>::DensityMatrix(const int ndim)
{
assert(ndim >= 0);
assert(ndim > 0);

dim_ = ndim;

Expand All @@ -45,6 +45,7 @@ DensityMatrix<MatrixType>::DensityMatrix(const int ndim)
kernel4dot_ = new MatrixType("K4dot", ndim, ndim);
work_ = new MatrixType("work", ndim, ndim);
occupation_.resize(dim_);
setDummyOcc();
}

template <class MatrixType>
Expand Down Expand Up @@ -109,6 +110,7 @@ void DensityMatrix<MatrixType>::build(
const std::vector<double>& occ, const int new_orbitals_index)
{
assert(dm_ != nullptr);
assert(!occ.empty());

setOccupations(occ);

Expand Down Expand Up @@ -149,6 +151,8 @@ template <class MatrixType>
void DensityMatrix<MatrixType>::setUniform(
const double nel, const int new_orbitals_index)
{
assert(!occupation_.empty());

MGmol_MPI& mmpi = *(MGmol_MPI::instance());
const double occ = (double)((double)nel / (double)dim_);
if (mmpi.instancePE0())
Expand Down Expand Up @@ -314,6 +318,7 @@ void DensityMatrix<MatrixType>::computeOccupations(const MatrixType& ls)
template <class MatrixType>
void DensityMatrix<MatrixType>::setOccupations(const std::vector<double>& occ)
{
assert(!occ.empty());
#ifdef PRINT_OPERATIONS
MGmol_MPI& mmpi = *(MGmol_MPI::instance());
if (mmpi.instancePE0())
Expand Down
13 changes: 11 additions & 2 deletions src/DensityMatrix.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,17 +84,24 @@ class DensityMatrix
*dm_ = mat;
orbitals_index_ = orbitals_index;

occupation_.clear();
setDummyOcc();

occ_uptodate_ = false;
uniform_occ_ = false;
stripped_ = false;
}

// set occupations to meaningless values to catch uninitialized use
void setDummyOcc()
{
for (auto& occ : occupation_)
occ = -1.;
}

void initMatrix(const double* const val)
{
dm_->init(val, dim_);
occupation_.clear();
setDummyOcc();

occ_uptodate_ = false;
uniform_occ_ = false;
Expand All @@ -105,8 +112,10 @@ class DensityMatrix

void getOccupations(std::vector<double>& occ) const
{
assert(!occupation_.empty());
assert(occ_uptodate_);
assert((int)occ.size() == dim_);

memcpy(&occ[0], &occupation_[0], dim_ * sizeof(double));
}

Expand Down
200 changes: 8 additions & 192 deletions src/Electrostatic.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#include "Control.h"
#include "ExtendedGridOrbitals.h"
#include "GridFactory.h"
#include "GridFunc.h"
#include "Hartree.h"
#include "Hartree_CG.h"
#include "Ions.h"
Expand All @@ -23,15 +24,6 @@
#include "ShiftedHartree.h"
#include "mputils.h"

#include "GridFunc.h"
#include "Laph2.h"
#include "Laph4.h"
#include "Laph4M.h"
#include "Laph4MP.h"
#include "Laph6.h"
#include "Laph8.h"
#include "ShiftedLaph4M.h"

Timer Electrostatic::solve_tm_("Electrostatic::solve");

Electrostatic::Electrostatic(PoissonFDtype lap_type, const short bcPoisson[3],
Expand All @@ -49,109 +41,9 @@ Electrostatic::Electrostatic(PoissonFDtype lap_type, const short bcPoisson[3],
Mesh* mymesh = Mesh::instance();
const pb::Grid& myGrid = mymesh->grid();

Control& ct = *(Control::instance());
if (ct.MGPoissonSolver()) // use MG for Poisson Solver
{
if (screening_const > 0.)
{
switch (lap_type)
{
case PoissonFDtype::h4M:
poisson_solver_
= new ShiftedHartree<pb::ShiftedLaph4M<POTDTYPE>>(
myGrid, bc_, screening_const);
break;
default:
(*MPIdata::sout)
<< "Electrostatic, shifted, Undefined option: "
<< static_cast<int>(lap_type) << std::endl;
}
}
else
{
switch (lap_type)
{
case PoissonFDtype::h4M:
poisson_solver_
= new Hartree<pb::Laph4M<POTDTYPE>>(myGrid, bc_);
break;
case PoissonFDtype::h2:
poisson_solver_
= new Hartree<pb::Laph2<POTDTYPE>>(myGrid, bc_);
break;
case PoissonFDtype::h4:
poisson_solver_
= new Hartree<pb::Laph4<POTDTYPE>>(myGrid, bc_);
break;
case PoissonFDtype::h6:
poisson_solver_
= new Hartree<pb::Laph6<POTDTYPE>>(myGrid, bc_);
break;
case PoissonFDtype::h8:
poisson_solver_
= new Hartree<pb::Laph8<POTDTYPE>>(myGrid, bc_);
break;
case PoissonFDtype::h4MP:
poisson_solver_
= new Hartree<pb::Laph4MP<POTDTYPE>>(myGrid, bc_);
break;
default:
(*MPIdata::sout) << "Electrostatic, Undefined option: "
<< static_cast<int>(lap_type) << std::endl;
}
}
}
else // use PCG for Poisson Solver
{
if (screening_const > 0.)
{
switch (lap_type)
{
case PoissonFDtype::h4M:
poisson_solver_
= new ShiftedHartree<pb::ShiftedLaph4M<POTDTYPE>>(
myGrid, bc_, screening_const);
break;
default:
(*MPIdata::sout)
<< "PCG Electrostatic, shifted, Undefined option: "
<< static_cast<int>(lap_type) << std::endl;
}
}
else
{
switch (lap_type)
{
case PoissonFDtype::h4M:
poisson_solver_
= new Hartree_CG<pb::Laph4M<POTDTYPE>>(myGrid, bc_);
break;
case PoissonFDtype::h2:
poisson_solver_
= new Hartree_CG<pb::Laph2<POTDTYPE>>(myGrid, bc_);
break;
case PoissonFDtype::h4:
poisson_solver_
= new Hartree_CG<pb::Laph4<POTDTYPE>>(myGrid, bc_);
break;
case PoissonFDtype::h6:
poisson_solver_
= new Hartree_CG<pb::Laph6<POTDTYPE>>(myGrid, bc_);
break;
case PoissonFDtype::h8:
poisson_solver_
= new Hartree_CG<pb::Laph8<POTDTYPE>>(myGrid, bc_);
break;
case PoissonFDtype::h4MP:
poisson_solver_
= new Hartree_CG<pb::Laph4MP<POTDTYPE>>(myGrid, bc_);
break;
default:
(*MPIdata::sout) << "PCG Electrostatic, Undefined option: "
<< static_cast<int>(lap_type) << std::endl;
}
}
}
// create Poisson solver
poisson_solver_ = PoissonSolverFactory::create(
myGrid, lap_type, bcPoisson, screening_const);

grhoc_ = nullptr;
diel_flag_ = false;
Expand Down Expand Up @@ -244,73 +136,8 @@ void Electrostatic::setupPB(
ngpts, origin, cell, static_cast<int>(laptype_), true, myPEenv);
if (poisson_solver_ != nullptr) delete poisson_solver_;

Control& ct = *(Control::instance());
if (ct.MGPoissonSolver()) // use MG for Poisson Solver
{
switch (laptype_)
{
case PoissonFDtype::h4M:
poisson_solver_ = new PBdiel<pb::PBh4M<POTDTYPE>>(
*pbGrid_, bc_, e0, rho0, drho0);
break;
case PoissonFDtype::h2:
poisson_solver_ = new PBdiel<pb::PBh2<POTDTYPE>>(
*pbGrid_, bc_, e0, rho0, drho0);
break;
case PoissonFDtype::h4:
poisson_solver_ = new PBdiel<pb::PBh4<POTDTYPE>>(
*pbGrid_, bc_, e0, rho0, drho0);
break;
case PoissonFDtype::h6:
poisson_solver_ = new PBdiel<pb::PBh6<POTDTYPE>>(
*pbGrid_, bc_, e0, rho0, drho0);
break;
case PoissonFDtype::h8:
poisson_solver_ = new PBdiel<pb::PBh8<POTDTYPE>>(
*pbGrid_, bc_, e0, rho0, drho0);
break;
case PoissonFDtype::h4MP:
poisson_solver_ = new PBdiel<pb::PBh4MP<POTDTYPE>>(
*pbGrid_, bc_, e0, rho0, drho0);
break;
default:
(*MPIdata::sout)
<< "Electrostatic, Undefined option" << std::endl;
}
}
else // use PCG for Poisson Solver
{
switch (laptype_)
{
case PoissonFDtype::h4M:
poisson_solver_ = new PBdiel_CG<pb::PBh4M<POTDTYPE>>(
*pbGrid_, bc_, e0, rho0, drho0);
break;
case PoissonFDtype::h2:
poisson_solver_ = new PBdiel_CG<pb::PBh2<POTDTYPE>>(
*pbGrid_, bc_, e0, rho0, drho0);
break;
case PoissonFDtype::h4:
poisson_solver_ = new PBdiel_CG<pb::PBh4<POTDTYPE>>(
*pbGrid_, bc_, e0, rho0, drho0);
break;
case PoissonFDtype::h6:
poisson_solver_ = new PBdiel_CG<pb::PBh6<POTDTYPE>>(
*pbGrid_, bc_, e0, rho0, drho0);
break;
case PoissonFDtype::h8:
poisson_solver_ = new PBdiel_CG<pb::PBh8<POTDTYPE>>(
*pbGrid_, bc_, e0, rho0, drho0);
break;
case PoissonFDtype::h4MP:
poisson_solver_ = new PBdiel_CG<pb::PBh4MP<POTDTYPE>>(
*pbGrid_, bc_, e0, rho0, drho0);
break;
default:
(*MPIdata::sout)
<< "Electrostatic, Undefined option" << std::endl;
}
}
poisson_solver_ = PoissonSolverFactory::createDiel(
*pbGrid_, laptype_, bc_, e0, rho0, drho0);

if (grhoc_ != nullptr)
{
Expand All @@ -330,6 +157,7 @@ void Electrostatic::setupPB(
poisson_solver_->set_vh(gf_vh);
}

// This function is only useful for Hartree problem with dielectric continuum
void Electrostatic::fillFuncAroundIons(const Ions& ions)
{
assert(grhod_ != nullptr);
Expand All @@ -352,7 +180,6 @@ void Electrostatic::fillFuncAroundIons(const Ions& ions)
std::vector<Ion*>::const_iterator ion = rc_ions.begin();
while (ion != rc_ions.end())
{

double rc = (*ion)->getRC();
// Special case: silicon
if ((*ion)->isMass28()) rc = 2.0;
Expand All @@ -373,43 +200,32 @@ void Electrostatic::fillFuncAroundIons(const Ions& ions)
#endif
for (unsigned int ix = 0; ix < pbGrid_->dim(0); ix++)
{

xc[1] = pbGrid_->start(1);
const int ix1 = (ix + shift) * incx;

for (unsigned int iy = 0; iy < pbGrid_->dim(1); iy++)
{

xc[2] = pbGrid_->start(2);
const int iy1 = ix1 + (iy + shift) * incy;

for (unsigned int iz = 0; iz < pbGrid_->dim(2); iz++)
{

const double r = (*ion)->minimage(xc, lattice, bc_);

if (r < rc)
{
const double alpha = 0.2 * (1. + cos(r * pi_rc));

const int iz1 = iy1 + iz + shift;
const int iz1 = iy1 + iz + shift;
vv[iz1] += alpha;
}

xc[2] += pbGrid_->hgrid(2);
}

xc[1] += pbGrid_->hgrid(1);

} // end for iy

xc[0] += pbGrid_->hgrid(0);

} // end for ix
}

ion++;

} // end loop on list of ions

return;
Expand Down
1 change: 1 addition & 0 deletions src/Electrostatic.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#include "Control.h"
#include "GridFunc.h"
#include "Poisson.h"
#include "PoissonSolverFactory.h"
#include "Rho.h"
#include "Timer.h"

Expand Down
Loading

0 comments on commit f0319f5

Please sign in to comment.