Skip to content

Commit

Permalink
use wavelet represenatation for xc potential in gga case
Browse files Browse the repository at this point in the history
  • Loading branch information
@moritzgubler
moritzgubler committed Aug 7, 2024
1 parent cf2fcc6 commit 6c75108
Show file tree
Hide file tree
Showing 3 changed files with 19 additions and 13 deletions.
2 changes: 1 addition & 1 deletion src/surface_forces/SurfaceForce.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -366,7 +366,7 @@ Eigen::MatrixXd surface_forces(mrchem::Molecule &mol, mrchem::OrbitalVector &Phi
VectorXd weights = integrator.getWeights();
MatrixXd normals = integrator.getNormals();

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

std::vector<Matrix3d> kstress = kineticStress(mol, Phi, nablaPhi, hessRho, prec, gridPos);
Expand Down
28 changes: 17 additions & 11 deletions src/surface_forces/xcStress.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ std::vector<Matrix3d> xcLDASpinStress(unique_ptr<mrdft::MRDFT> &mrdft_p, MatrixX
* @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){
std::vector<Matrix3d> xcGGAStress(unique_ptr<mrdft::MRDFT> &mrdft_p, mrcpp::FunctionTreeVector<3> &xc_pots, MatrixXd &rhoGrid, MatrixXd &nablaRhoGrid, Eigen::MatrixXd &gridPos){
int nGrid = rhoGrid.rows();
Eigen::MatrixXd inp(rhoGrid.rows(), 4);
inp.col(0) = rhoGrid.col(0);
Expand All @@ -74,15 +74,18 @@ std::vector<Matrix3d> xcGGAStress(unique_ptr<mrdft::MRDFT> &mrdft_p, MatrixXd &r
inp.col(3) = nablaRhoGrid.col(2);
Eigen::MatrixXd xcOUT = mrdft_p->functional().evaluate_transposed(inp);
std::vector<Matrix3d> out(nGrid);
std::array<double, 3> pos;
for (int i = 0; i < rhoGrid.rows(); i++) {
out[i] = Matrix3d::Zero();

pos[0] = gridPos(i, 0);
pos[1] = gridPos(i, 1);
pos[2] = gridPos(i, 2);
for (int j = 0; j < 3; j++) {
out[i](j, j) = xcOUT(i, 0) - xcOUT(i, 1) * rhoGrid(i);
out[i](j, j) = xcOUT(i, 0) - rhoGrid(i) * std::get<1>(xc_pots[1])->evalf(pos);
}
for (int j1 = 0; j1 < 3; j1++) {
for (int j2 = 0; j2 < 3; j2++) {
out[i](j1, j2) = out[i](j1, j2) - xcOUT(i, 2 + j1) * nablaRhoGrid(i, j2) * rhoGrid(i);
out[i](j1, j2) = out[i](j1, j2) - xcOUT(i, 2 + j1) * nablaRhoGrid(i, j2);
}
}
}
Expand All @@ -99,7 +102,7 @@ std::vector<Matrix3d> xcGGAStress(unique_ptr<mrdft::MRDFT> &mrdft_p, MatrixXd &r
* @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){
std::vector<Matrix3d> xcGGASpinStress(unique_ptr<mrdft::MRDFT> &mrdft_p, mrcpp::FunctionTreeVector<3> &xc_pots, MatrixXd &rhoGridAlpha, MatrixXd &rhoGridBeta, MatrixXd &nablaRhoGridAlpha, MatrixXd &nablaRhoGridBeta, Eigen::MatrixXd &gridPos){
int nGrid = rhoGridAlpha.rows();
Eigen::MatrixXd inp(rhoGridAlpha.rows(), 8);
std::vector<Matrix3d> out = std::vector<Eigen::Matrix3d>(nGrid);
Expand All @@ -112,16 +115,19 @@ std::vector<Matrix3d> xcGGASpinStress(unique_ptr<mrdft::MRDFT> &mrdft_p, MatrixX
inp.col(6) = nablaRhoGridBeta.col(1);
inp.col(7) = nablaRhoGridBeta.col(2);
Eigen::MatrixXd xc = mrdft_p->functional().evaluate_transposed(inp);
std::array<double, 3> pos;
for (int i = 0; i < rhoGridAlpha.rows(); i++) {
out[i] = Matrix3d::Zero();
pos[0] = gridPos(i, 0);
pos[1] = gridPos(i, 1);
pos[2] = gridPos(i, 2);
for (int j = 0; j < 3; j++) {
out[i](j, j) = xc(i, 0) - xc(i, 1) * rhoGridAlpha(i) - xc(i, 2) * rhoGridBeta(i);
out[i](j, j) = xc(i, 0) - std::get<1>(xc_pots[1])->evalf(pos) * rhoGridAlpha(i) - std::get<1>(xc_pots[2])->evalf(pos) * rhoGridBeta(i);
}
for (int j1 = 0; j1 < 3; j1++) {
for (int j2 = 0; j2 < 3; j2++) {
out[i](j1, j2) = out[i](j1, j2)
- xc(i, 3 + j1) * nablaRhoGridAlpha(i, j2) * rhoGridAlpha(i)
- xc(i, 6 + j1) * nablaRhoGridBeta(i, j2) * rhoGridBeta(i);
- xc(i, 3 + j1) * nablaRhoGridAlpha(i, j2) - xc(i, 6 + j1) * nablaRhoGridBeta(i, j2);
}
}
}
Expand All @@ -138,7 +144,7 @@ std::vector<Matrix3d> xcGGASpinStress(unique_ptr<mrdft::MRDFT> &mrdft_p, MatrixX
* @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,
std::vector<Eigen::Matrix3d> getXCStress(unique_ptr<mrdft::MRDFT> &mrdft_p, mrcpp::FunctionTreeVector<3> &xc_pots, std::shared_ptr<OrbitalVector> phi,
std::shared_ptr<NablaOperator> nabla, MatrixXd &gridPos, bool isOpenShell, double prec){

bool isGGA = mrdft_p->functional().isGGA();
Expand Down Expand Up @@ -187,7 +193,7 @@ std::vector<Eigen::Matrix3d> getXCStress(unique_ptr<mrdft::MRDFT> &mrdft_p, std:
nablaRhoGridBeta(i, 2) = nablaRhoBeta[2].real().evalf(pos);
}

xcStress = xcGGASpinStress(mrdft_p, rhoGridAlpha, rhoGridBeta, nablaRhoGridAlpha, nablaRhoGridBeta);
xcStress = xcGGASpinStress(mrdft_p, xc_pots, rhoGridAlpha, rhoGridBeta, nablaRhoGridAlpha, nablaRhoGridBeta, gridPos);
} else {
xcStress = xcLDASpinStress(mrdft_p, rhoGridAlpha, rhoGridBeta);
}
Expand Down Expand Up @@ -216,7 +222,7 @@ std::vector<Eigen::Matrix3d> getXCStress(unique_ptr<mrdft::MRDFT> &mrdft_p, std:
nablaRhoGrid(i, 1) = nablaRho[1].real().evalf(pos);
nablaRhoGrid(i, 2) = nablaRho[2].real().evalf(pos);
}
xcStress = xcGGAStress(mrdft_p, rhoGrid, nablaRhoGrid);
xcStress = xcGGAStress(mrdft_p, xc_pots, rhoGrid, nablaRhoGrid, gridPos);
} else {
xcStress = xcLDAStress(mrdft_p, rhoGrid);
}
Expand Down
2 changes: 1 addition & 1 deletion src/surface_forces/xcStress.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,6 @@ std::vector<Eigen::Matrix3d> xcLDAStress(std::unique_ptr<mrdft::MRDFT> &mrdft_p,
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);
std::vector<Eigen::Matrix3d> getXCStress(std::unique_ptr<mrdft::MRDFT> &mrdft_p, mrcpp::FunctionTreeVector<3> &xc_pots, 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 6c75108

Please sign in to comment.