From 6c751081d6e99ebe034a4a8136d7f0d35166c474 Mon Sep 17 00:00:00 2001 From: moritzgubler Date: Wed, 7 Aug 2024 10:14:45 +0200 Subject: [PATCH] use wavelet represenatation for xc potential in gga case --- src/surface_forces/SurfaceForce.cpp | 2 +- src/surface_forces/xcStress.cpp | 28 +++++++++++++++++----------- src/surface_forces/xcStress.h | 2 +- 3 files changed, 19 insertions(+), 13 deletions(-) diff --git a/src/surface_forces/SurfaceForce.cpp b/src/surface_forces/SurfaceForce.cpp index 5fab1d150..821f18c9a 100644 --- a/src/surface_forces/SurfaceForce.cpp +++ b/src/surface_forces/SurfaceForce.cpp @@ -366,7 +366,7 @@ Eigen::MatrixXd surface_forces(mrchem::Molecule &mol, mrchem::OrbitalVector &Phi VectorXd weights = integrator.getWeights(); MatrixXd normals = integrator.getNormals(); - std::vector xcStress = getXCStress(mrdft_p, std::make_shared(Phi), + std::vector xcStress = getXCStress(mrdft_p, *xc_pot_vector, std::make_shared(Phi), std::make_shared(nabla), gridPos, xc_spin, prec); std::vector kstress = kineticStress(mol, Phi, nablaPhi, hessRho, prec, gridPos); diff --git a/src/surface_forces/xcStress.cpp b/src/surface_forces/xcStress.cpp index b6fc16a4a..ee9422773 100644 --- a/src/surface_forces/xcStress.cpp +++ b/src/surface_forces/xcStress.cpp @@ -65,7 +65,7 @@ std::vector xcLDASpinStress(unique_ptr &mrdft_p, MatrixX * @param nablaRhoGrid MatrixXd with gradient of density values, shape (nGrid, 3) * @return std::vector vector of 3x3 matrices with stress tensor at each grid point */ -std::vector xcGGAStress(unique_ptr &mrdft_p, MatrixXd &rhoGrid, MatrixXd &nablaRhoGrid){ +std::vector xcGGAStress(unique_ptr &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); @@ -74,15 +74,18 @@ std::vector xcGGAStress(unique_ptr &mrdft_p, MatrixXd &r inp.col(3) = nablaRhoGrid.col(2); Eigen::MatrixXd xcOUT = mrdft_p->functional().evaluate_transposed(inp); std::vector out(nGrid); + std::array 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); } } } @@ -99,7 +102,7 @@ std::vector xcGGAStress(unique_ptr &mrdft_p, MatrixXd &r * @param nablaRhoGridBeta MatrixXd with gradient of beta density values, shape (nGrid, 3) * @return std::vector vector of 3x3 matrices with stress tensor at each grid point */ -std::vector xcGGASpinStress(unique_ptr &mrdft_p, MatrixXd &rhoGridAlpha, MatrixXd &rhoGridBeta, MatrixXd &nablaRhoGridAlpha, MatrixXd &nablaRhoGridBeta){ +std::vector xcGGASpinStress(unique_ptr &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 out = std::vector(nGrid); @@ -112,16 +115,19 @@ std::vector xcGGASpinStress(unique_ptr &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 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); } } } @@ -138,7 +144,7 @@ std::vector xcGGASpinStress(unique_ptr &mrdft_p, MatrixX * @param isOpenShell bool, true if open shell calculation * @param prec precision to use in density representation */ -std::vector getXCStress(unique_ptr &mrdft_p, std::shared_ptr phi, +std::vector getXCStress(unique_ptr &mrdft_p, mrcpp::FunctionTreeVector<3> &xc_pots, std::shared_ptr phi, std::shared_ptr nabla, MatrixXd &gridPos, bool isOpenShell, double prec){ bool isGGA = mrdft_p->functional().isGGA(); @@ -187,7 +193,7 @@ std::vector getXCStress(unique_ptr &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); } @@ -216,7 +222,7 @@ std::vector getXCStress(unique_ptr &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); } diff --git a/src/surface_forces/xcStress.h b/src/surface_forces/xcStress.h index c498b47b2..5eb19d892 100644 --- a/src/surface_forces/xcStress.h +++ b/src/surface_forces/xcStress.h @@ -14,6 +14,6 @@ std::vector xcLDAStress(std::unique_ptr &mrdft_p, std::vector xcLDASpinStress(std::unique_ptr &mrdft_p, Eigen::MatrixXd &rhoGridAlpha, Eigen::MatrixXd &rhoGridBeta); std::vector xcGGAStress(std::unique_ptr &mrdft_p, Eigen::MatrixXd &rhoGrid, Eigen::MatrixXd &nablaRhoGrid); std::vector xcGGASpinStress(std::unique_ptr &mrdft_p, Eigen::MatrixXd &rhoGridAlpha, Eigen::MatrixXd &rhoGridBeta, Eigen::MatrixXd &nablaRhoGridAlpha, Eigen::MatrixXd &nablaRhoGridBeta); -std::vector getXCStress(std::unique_ptr &mrdft_p, std::shared_ptr phi, std::shared_ptr nabla, Eigen::MatrixXd &gridPos, bool isOpenShell, double prec); +std::vector getXCStress(std::unique_ptr &mrdft_p, mrcpp::FunctionTreeVector<3> &xc_pots, std::shared_ptr phi, std::shared_ptr nabla, Eigen::MatrixXd &gridPos, bool isOpenShell, double prec); } // namespace surface_force