Skip to content

Commit

Permalink
gga working
Browse files Browse the repository at this point in the history
  • Loading branch information
@moritzgubler
moritzgubler committed Jul 31, 2024
1 parent 62d9774 commit 2a1da0c
Show file tree
Hide file tree
Showing 3 changed files with 38 additions and 61 deletions.
47 changes: 13 additions & 34 deletions src/surface_forces/SurfaceForce.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,34 +159,6 @@ std::vector<Eigen::Matrix3d> maxwellStress(const Molecule &mol, mrchem::OrbitalV
return stress;
}

// /**
// * @brief Calculates the exchange-correlation stress tensor for the given molecule.
// */
// std::vector<Matrix3d> xcStress(const Eigen::MatrixXd xcGrid, bool isGGA){
// int nGrid = xcGrid.cols();
// std::cout << "xcGrid: " << xcGrid.rows() << " " << xcGrid.cols() << std::endl;

// std::vector<Matrix3d> stress(nGrid);

// if (!isGGA) {
// for (int i = 0; i < nGrid; i++) {
// for (int i1 = 0; i1 < 3; i1++) {
// for (int i2 = 0; i2 < 3; i2++) {
// stress[i](i1, i2) = 0.0;
// }
// }
// for (int i1 = 0; i1 < 3; i1++) {
// stress[i](i1, i1) = xcGrid(0, i) - xcGrid(1, i);
// }
// // std::cout << "stress: " << xcGrid(0, 1) << " " << xcGrid(1, 0) << std::endl;
// }
// } else {
// MSG_ABORT("GGA not implemented");
// }

// return stress;
// }

/**
* @brief Calculates the kinetic stress tensor for the given molecule. See the function description for the formula.
*/
Expand Down Expand Up @@ -412,10 +384,10 @@ Eigen::MatrixXd surface_forces(mrchem::Molecule &mol, mrchem::OrbitalVector &Phi
nLebPoints = 194;
}
else if (leb_prec == "medium"){
nLebPoints = 350;
nLebPoints = 434;
}
else if (leb_prec == "high") {
nLebPoints = 590;
nLebPoints = 770;
}
else {
MSG_ABORT("Invalid lebedev precision");
Expand Down Expand Up @@ -457,7 +429,6 @@ Eigen::MatrixXd surface_forces(mrchem::Molecule &mol, mrchem::OrbitalVector &Phi
Eigen::MatrixXd rhoGrid(integrator.n, 1);
Eigen::MatrixXd rhoGridAlpha(integrator.n, 1);
Eigen::MatrixXd rhoGridBeta(integrator.n, 1);
std::cout << "density start" << std::endl;
for (int i = 0; i < integrator.n; i++) {
pos[0] = gridPos(i, 0);
pos[1] = gridPos(i, 1);
Expand All @@ -476,9 +447,17 @@ Eigen::MatrixXd surface_forces(mrchem::Molecule &mol, mrchem::OrbitalVector &Phi
}
} else{
if ( !xc_spin) {
// mrcpp::ComplexFunction nablaRho = nabla(rho.real());
// Eigen::MatrixXd nablaRhoGrid(integrator.n, 3);
// xcStress = xcGGA(mrdft_p, rhoGrid, nablaRhoGrid);
mrchem::OrbitalVector nablaRho = nabla(rho);
Eigen::MatrixXd nablaRhoGrid(integrator.n, 3);
for (int i = 0; i < integrator.n; i++) {
pos[0] = gridPos(i, 0);
pos[1] = gridPos(i, 1);
pos[2] = gridPos(i, 2);
nablaRhoGrid(i, 0) = nablaRho[0].real().evalf(pos);
nablaRhoGrid(i, 1) = nablaRho[1].real().evalf(pos);
nablaRhoGrid(i, 2) = nablaRho[2].real().evalf(pos);
}
xcStress = xcGGA(mrdft_p, rhoGrid, nablaRhoGrid);
} else {
// xcStress = xcGGASpin(mrdft_p, rhoGridAlpha, rhoGridBeta, nablaRhoGridAlpha, nablaRhoGridBeta);
}
Expand Down
48 changes: 22 additions & 26 deletions src/surface_forces/xcHelper.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,33 +9,24 @@ using namespace Eigen;
std::vector<Eigen::Matrix3d> xcLDA(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGrid){
int nGrid = rhoGrid.rows();
std::vector<Eigen::Matrix3d> out(nGrid);
// Eigen::MatrixXd out = Eigen::MatrixXd::Zero(rhoGrid.rows(), 5);
Eigen::MatrixXd xcOUT = mrdft_p->functional().evaluate_transposed(rhoGrid);
// out.col(0) = xcOUT.col(0);
for (int i = 0; i < nGrid; i++) {
// out(i, 1) = xcOUT(i, 1) * rhoGrid(i);
out[i] = Matrix3d::Zero();
for (int j = 0; j < 3; j++) {
out[i](j, j) = xcOUT(i, 0) - xcOUT(i, 1) * rhoGrid(i);
}
}

return out;
}

std::vector<Matrix3d> xcLDASpin(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGridAlpha, Eigen::MatrixXd &rhoGridBeta){
int nGrid = rhoGridAlpha.rows();
Eigen::MatrixXd inp(rhoGridAlpha.rows(), 2);
// Eigen::MatrixXd out = Eigen::MatrixXd::Zero(rhoGridAlpha.rows(), 5);
std::vector<Matrix3d> out = std::vector<Eigen::Matrix3d>(nGrid);
// Eigen::MatrixXd outAlpha = Eigen::MatrixXd::Zero(rhoGridAlpha.rows(), 5);
// Eigen::MatrixXd outBeta = Eigen::MatrixXd::Zero(rhoGridAlpha.rows(), 5);
inp.col(0) = rhoGridAlpha.col(0);
inp.col(1) = rhoGridBeta.col(0);
Eigen::MatrixXd xc = mrdft_p->functional().evaluate_transposed(inp);
// out.col(0) = xc.col(0);
for (int i = 0; i < rhoGridAlpha.rows(); i++) {
// out(i, 1) = xc(i, 1) * rhoGridAlpha(i) + xc(i, 2) * rhoGridBeta(i);
out[i] = Matrix3d::Zero();
for (int j = 0; j < 3; j++) {
out[i](j, j) = xc(i, 0) - xc(i, 1) * rhoGridAlpha(i) - xc(i, 2) * rhoGridBeta(i);
Expand All @@ -44,22 +35,27 @@ std::vector<Matrix3d> xcLDASpin(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::M
return out;
}

// Eigen::MatrixXd xcGGA(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGrid, Eigen::MatrixXd &nablaRhoGrid){
// Eigen::MatrixXd out = Eigen::MatrixXd::Zero(rhoGrid.rows(), 5);
// Eigen::MatrixXd inp(rhoGrid.rows(), 4);
// inp.col(0) = rhoGrid.col(0);
// inp.col(1) = nablaRhoGrid.col(0);
// inp.col(2) = nablaRhoGrid.col(1);
// inp.col(3) = nablaRhoGrid.col(2);
// Eigen::MatrixXd xcOUT = mrdft_p->functional().evaluate_transposed(inp);
// out.col(0) = xcOUT.col(0);
// for (int i = 0; i < rhoGrid.rows(); i++) {
// out(i, 1) = xcOUT(i, 1) * rhoGrid(i);
// out(i, 2) = xcOUT(i, 2) * nablaRhoGrid(i, 0);
// out(i, 3) = xcOUT(i, 3) * nablaRhoGrid(i, 1);
// out(i, 4) = xcOUT(i, 4) * nablaRhoGrid(i, 2);
// }
std::vector<Matrix3d> xcGGA(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGrid, Eigen::MatrixXd &nablaRhoGrid){
int nGrid = rhoGrid.rows();
Eigen::MatrixXd inp(rhoGrid.rows(), 4);
inp.col(0) = rhoGrid.col(0);
inp.col(1) = nablaRhoGrid.col(0);
inp.col(2) = nablaRhoGrid.col(1);
inp.col(3) = nablaRhoGrid.col(2);
Eigen::MatrixXd xcOUT = mrdft_p->functional().evaluate_transposed(inp);
std::vector<Matrix3d> out(nGrid);
for (int i = 0; i < rhoGrid.rows(); i++) {
out[i] = Matrix3d::Zero();

// return out;
// }
for (int j = 0; j < 3; j++) {
out[i](j, j) = xcOUT(i, 0) - xcOUT(i, 1) * rhoGrid(i);
}
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);
}
}
}
return out;
}

4 changes: 3 additions & 1 deletion src/surface_forces/xcHelper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,4 +8,6 @@
#include <vector>

std::vector<Eigen::Matrix3d> xcLDA(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGrid);
std::vector<Eigen::Matrix3d> xcLDASpin(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGridAlpha, Eigen::MatrixXd &rhoGridBeta);
std::vector<Eigen::Matrix3d> xcLDASpin(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGridAlpha, Eigen::MatrixXd &rhoGridBeta);
std::vector<Eigen::Matrix3d> xcGGA(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGrid, Eigen::MatrixXd &nablaRhoGrid);
std::vector<Eigen::Matrix3d> xcGGASpin(std::unique_ptr<mrdft::MRDFT> &mrdft_p, Eigen::MatrixXd &rhoGridAlpha, Eigen::MatrixXd &rhoGridBeta, Eigen::MatrixXd &nablaRhoGridAlpha, Eigen::MatrixXd &nablaRhoGridBeta);

0 comments on commit 2a1da0c

Please sign in to comment.