Skip to content

Commit

Permalink
add more documentation and namespaces
Browse files Browse the repository at this point in the history
  • Loading branch information
@moritzgubler
moritzgubler committed Jul 31, 2024
1 parent 2bfcf2f commit 13b89d8
Show file tree
Hide file tree
Showing 3 changed files with 53 additions and 30 deletions.
23 changes: 1 addition & 22 deletions src/surface_forces/SurfaceForce.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -304,11 +304,7 @@ Eigen::MatrixXd surface_forces(mrchem::Molecule &mol, mrchem::OrbitalVector &Phi

// setup density
mrchem::Density rho(false);
mrchem::Density rhoA(false);
mrchem::Density rhoB(false);
mrchem::density::compute(prec, rho, Phi, DensityType::Total);

Check warning on line 307 in src/surface_forces/SurfaceForce.cpp

View check run for this annotation

Codecov / codecov/patch

src/surface_forces/SurfaceForce.cpp#L306-L307 

Added lines #L306 - L307 were not covered by tests
mrchem::density::compute(prec, rhoA, Phi, DensityType::Alpha);
mrchem::density::compute(prec, rhoB, Phi, DensityType::Beta);


// setup operators and potentials
Expand Down Expand Up @@ -356,13 +352,7 @@ Eigen::MatrixXd surface_forces(mrchem::Molecule &mol, mrchem::OrbitalVector &Phi
}
std::unique_ptr<mrdft::MRDFT> mrdft_p = xc_factory.build();

Check warning on line 353 in src/surface_forces/SurfaceForce.cpp

View check run for this annotation

Codecov / codecov/patch

src/surface_forces/SurfaceForce.cpp#L353 

Added line #L353 was not covered by tests

bool isGGA = mrdft_p->functional().isGGA();
bool isHybrid = mrdft_p->functional().isHybrid();



int numAtoms = mol.getNNuclei();

Check warning on line 355 in src/surface_forces/SurfaceForce.cpp

View check run for this annotation

Codecov / codecov/patch

src/surface_forces/SurfaceForce.cpp#L355 

Added line #L355 was not covered by tests
int numOrbitals = Phi.size();

Eigen::MatrixXd forces = Eigen::MatrixXd::Zero(numAtoms, 3);
Vector3d center;
Expand Down Expand Up @@ -426,20 +416,9 @@ Eigen::MatrixXd surface_forces(mrchem::Molecule &mol, mrchem::OrbitalVector &Phi
MatrixXd gridPos = integrator.getPoints();
VectorXd weights = integrator.getWeights();
MatrixXd normals = integrator.getNormals();

Check warning on line 418 in src/surface_forces/SurfaceForce.cpp

View check run for this annotation

Codecov / codecov/patch

src/surface_forces/SurfaceForce.cpp#L414-L418 

Added lines #L414 - L418 were not covered by tests
Eigen::MatrixXd rhoGrid(integrator.n, 1);
Eigen::MatrixXd rhoGridAlpha(integrator.n, 1);
Eigen::MatrixXd rhoGridBeta(integrator.n, 1);
for (int i = 0; i < integrator.n; i++) {
pos[0] = gridPos(i, 0);
pos[1] = gridPos(i, 1);
pos[2] = gridPos(i, 2);
rhoGrid(i, 0) = rho.real().evalf(pos);
rhoGridAlpha(i, 0) = rhoA.real().evalf(pos);
rhoGridBeta(i, 0) = rhoB.real().evalf(pos);
}

std::vector<Matrix3d> xcStress = getXCStress(mrdft_p, std::make_shared<mrchem::OrbitalVector>(Phi),
std::make_shared<mrchem::NablaOperator>(nabla), gridPos, xc_spin, prec);
std::make_shared<mrchem::NablaOperator>(nabla), gridPos, xc_spin, prec);

Check warning on line 421 in src/surface_forces/SurfaceForce.cpp

View check run for this annotation

Codecov / codecov/patch

src/surface_forces/SurfaceForce.cpp#L420-L421 

Added lines #L420 - L421 were not covered by tests

std::vector<Matrix3d> kstress = kineticStress(mol, Phi, nablaPhi, hessRho, prec, gridPos);
std::vector<Matrix3d> mstress = maxwellStress(mol, negEfield, gridPos, prec);
Expand Down
57 changes: 49 additions & 8 deletions src/surface_forces/xcStress.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,15 @@ using namespace Eigen;
using namespace mrchem;
using namespace std;

namespace surface_force{

/**
* @brief Compute the exchange-correlation stress tensor for LDA functional
*
* @param mrdft_p MRDFT object
* @param rhoGrid MatrixXd with density values, shape (nGrid, 1)
* @return std::vector<Eigen::Matrix3d> vector of 3x3 matrices with stress tensor at each grid point
*/
std::vector<Eigen::Matrix3d> xcLDAStress(unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGrid){
int nGrid = rhoGrid.rows();
std::vector<Eigen::Matrix3d> out(nGrid);
Expand All @@ -24,6 +33,14 @@ std::vector<Eigen::Matrix3d> xcLDAStress(unique_ptr<mrdft::MRDFT> &mrdft_p, Eige
return out;

Check warning on line 33 in src/surface_forces/xcStress.cpp

View check run for this annotation

Codecov / codecov/patch

src/surface_forces/xcStress.cpp#L33 

Added line #L33 was not covered by tests
}

/**
* @brief Compute the exchange-correlation stress tensor for LDA functional for open shell systems
*
* @param mrdft_p MRDFT object
* @param rhoGridAlpha MatrixXd with alpha density values, shape (nGrid, 1)
* @param rhoGridBeta MatrixXd with beta density values, shape (nGrid, 1)
* @return std::vector<Eigen::Matrix3d> vector of 3x3 matrices with stress tensor at each grid point
*/
std::vector<Matrix3d> xcLDASpinStress(unique_ptr<mrdft::MRDFT> &mrdft_p, MatrixXd &rhoGridAlpha, MatrixXd &rhoGridBeta){
int nGrid = rhoGridAlpha.rows();
Eigen::MatrixXd inp(rhoGridAlpha.rows(), 2);
Expand All @@ -40,6 +57,14 @@ std::vector<Matrix3d> xcLDASpinStress(unique_ptr<mrdft::MRDFT> &mrdft_p, MatrixX
return out;

Check warning on line 57 in src/surface_forces/xcStress.cpp

View check run for this annotation

Codecov / codecov/patch

src/surface_forces/xcStress.cpp#L57 

Added line #L57 was not covered by tests
}

/**
* @brief Compute the exchange-correlation stress tensor for GGA functional
*
* @param mrdft_p MRDFT object
* @param rhoGrid MatrixXd with density values, shape (nGrid, 1)
* @param nablaRhoGrid MatrixXd with gradient of density values, shape (nGrid, 3)
* @return std::vector<Eigen::Matrix3d> vector of 3x3 matrices with stress tensor at each grid point
*/
std::vector<Matrix3d> xcGGAStress(unique_ptr<mrdft::MRDFT> &mrdft_p, MatrixXd &rhoGrid, MatrixXd &nablaRhoGrid){
int nGrid = rhoGrid.rows();
Eigen::MatrixXd inp(rhoGrid.rows(), 4);
Expand All @@ -64,6 +89,16 @@ std::vector<Matrix3d> xcGGAStress(unique_ptr<mrdft::MRDFT> &mrdft_p, MatrixXd &r
return out;

Check warning on line 89 in src/surface_forces/xcStress.cpp

View check run for this annotation

Codecov / codecov/patch

src/surface_forces/xcStress.cpp#L89 

Added line #L89 was not covered by tests
}

/**
* @brief Compute the exchange-correlation stress tensor for GGA functional for open shell systems
*
* @param mrdft_p MRDFT object
* @param rhoGridAlpha MatrixXd with alpha density values, shape (nGrid, 1)
* @param rhoGridBeta MatrixXd with beta density values, shape (nGrid, 1)
* @param nablaRhoGridAlpha MatrixXd with gradient of alpha density values, shape (nGrid, 3)
* @param nablaRhoGridBeta MatrixXd with gradient of beta density values, shape (nGrid, 3)
* @return std::vector<Eigen::Matrix3d> vector of 3x3 matrices with stress tensor at each grid point
*/
std::vector<Matrix3d> xcGGASpinStress(unique_ptr<mrdft::MRDFT> &mrdft_p, MatrixXd &rhoGridAlpha, MatrixXd &rhoGridBeta, MatrixXd &nablaRhoGridAlpha, MatrixXd &nablaRhoGridBeta){
int nGrid = rhoGridAlpha.rows();
Eigen::MatrixXd inp(rhoGridAlpha.rows(), 8);
Expand Down Expand Up @@ -93,7 +128,18 @@ std::vector<Matrix3d> xcGGASpinStress(unique_ptr<mrdft::MRDFT> &mrdft_p, MatrixX
return out;

Check warning on line 128 in src/surface_forces/xcStress.cpp

View check run for this annotation

Codecov / codecov/patch

src/surface_forces/xcStress.cpp#L128 

Added line #L128 was not covered by tests
}

std::vector<Eigen::Matrix3d> getXCStress(unique_ptr<mrdft::MRDFT> &mrdft_p, std::shared_ptr<OrbitalVector> phi, std::shared_ptr<NablaOperator> nabla, MatrixXd &gridPos, bool isOpenShell, double prec){
/**
* @brief Compute the exchange-correlation stress tensor on a grid
*
* @param mrdft_p MRDFT object
* @param phi OrbitalVector
* @param nabla NablaOperator (must be set up prior to calling this function)
* @param gridPos MatrixXd with grid positions, shape (nGrid, 3)
* @param isOpenShell bool, true if open shell calculation
* @param prec precision to use in density representation
*/
std::vector<Eigen::Matrix3d> getXCStress(unique_ptr<mrdft::MRDFT> &mrdft_p, std::shared_ptr<OrbitalVector> phi,

Check warning on line 141 in src/surface_forces/xcStress.cpp

View check run for this annotation

Codecov / codecov/patch

src/surface_forces/xcStress.cpp#L141 

Added line #L141 was not covered by tests
std::shared_ptr<NablaOperator> nabla, MatrixXd &gridPos, bool isOpenShell, double prec){

bool isGGA = mrdft_p->functional().isGGA();
bool isHybrid = mrdft_p->functional().isHybrid();
Expand Down Expand Up @@ -123,7 +169,6 @@ std::vector<Eigen::Matrix3d> getXCStress(unique_ptr<mrdft::MRDFT> &mrdft_p, std:
}

if (isGGA) {

mrchem::NablaOperator nablaOP = *nabla;
OrbitalVector nablaRhoAlpha = nablaOP(rhoA);
OrbitalVector nablaRhoBeta = nablaOP(rhoB);
Expand All @@ -147,7 +192,6 @@ std::vector<Eigen::Matrix3d> getXCStress(unique_ptr<mrdft::MRDFT> &mrdft_p, std:
xcStress = xcLDASpinStress(mrdft_p, rhoGridAlpha, rhoGridBeta);

Check warning on line 192 in src/surface_forces/xcStress.cpp

View check run for this annotation

Codecov / codecov/patch

src/surface_forces/xcStress.cpp#L190-L192 

Added lines #L190 - L192 were not covered by tests
}


} else { // closed shell
MatrixXd rhoGrid(nGrid, 1);
mrchem::Density rho(false);
Expand Down Expand Up @@ -176,10 +220,7 @@ std::vector<Eigen::Matrix3d> getXCStress(unique_ptr<mrdft::MRDFT> &mrdft_p, std:
} else {
xcStress = xcLDAStress(mrdft_p, rhoGrid);

Check warning on line 221 in src/surface_forces/xcStress.cpp

View check run for this annotation

Codecov / codecov/patch

src/surface_forces/xcStress.cpp#L219-L221 

Added lines #L219 - L221 were not covered by tests
}

}

return xcStress;

Check warning on line 224 in src/surface_forces/xcStress.cpp

View check run for this annotation

Codecov / codecov/patch

src/surface_forces/xcStress.cpp#L224 

Added line #L224 was not covered by tests


}
}
} // namespace surface_force
3 changes: 3 additions & 0 deletions src/surface_forces/xcStress.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,9 +8,12 @@
#include <vector>
#include "qmoperators/one_electron/NablaOperator.h"

namespace surface_force {

std::vector<Eigen::Matrix3d> xcLDAStress(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGrid);
std::vector<Eigen::Matrix3d> xcLDASpinStress(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGridAlpha, Eigen::MatrixXd &rhoGridBeta);
std::vector<Eigen::Matrix3d> xcGGAStress(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGrid, Eigen::MatrixXd &nablaRhoGrid);
std::vector<Eigen::Matrix3d> xcGGASpinStress(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGridAlpha, Eigen::MatrixXd &rhoGridBeta, Eigen::MatrixXd &nablaRhoGridAlpha, Eigen::MatrixXd &nablaRhoGridBeta);
std::vector<Eigen::Matrix3d> getXCStress(std::unique_ptr<mrdft::MRDFT> &mrdft_p, std::shared_ptr<mrchem::OrbitalVector> phi, std::shared_ptr<mrchem::NablaOperator> nabla, Eigen::MatrixXd &gridPos, bool isOpenShell, double prec);

} // namespace surface_force

0 comments on commit 13b89d8

Please sign in to comment.