|
| 1 | +#include "force_stress_lcao.h" |
| 2 | + |
| 3 | +#include "source_base/parallel_reduce.h" |
| 4 | +#include "source_lcao/module_dftu/dftu_lcao.h" //Quxin add for DFT+U on 20201029 |
| 5 | +#include "source_lcao/module_dftu/dftu_force.h" |
| 6 | +#include "source_io/module_output/output_log.h" |
| 7 | +#include "source_io/module_parameter/parameter.h" |
| 8 | +// new |
| 9 | +#include "source_base/timer.h" |
| 10 | +#include "source_base/tool_quit.h" |
| 11 | +#include "source_cell/module_neighbor/sltk_grid_driver.h" |
| 12 | +#include "source_estate/elecstate_lcao.h" |
| 13 | +#include "source_estate/module_pot/h_tddft_pw.h" // Taoni add 2025-02-20 |
| 14 | +#include "source_estate/module_pot/efield.h" // liuyu add 2022-05-18 |
| 15 | +#include "source_estate/module_pot/gatefield.h" // liuyu add 2022-09-13 |
| 16 | +#include "source_hamilt/module_surchem/surchem.h" //sunml add 2022-08-10 |
| 17 | +#include "source_hamilt/module_vdw/vdw.h" |
| 18 | +#include "source_io/module_parameter/parameter.h" |
| 19 | +#ifdef __MLALGO |
| 20 | +#include "source_lcao/module_deepks/lcao_deepks.h" //caoyu add for deepks 2021-06-03 |
| 21 | +#include "source_lcao/module_deepks/lcao_deepks_io.h" // mohan add 2024-07-22 |
| 22 | +#include "source_lcao/module_deepks/deepks_force.h" |
| 23 | +#endif |
| 24 | +#include "source_lcao/module_dftu/dftu_lcao_op.h" |
| 25 | +#include "source_lcao/module_operator_lcao/dspin_lcao.h" |
| 26 | +#include "source_lcao/module_operator_lcao/nonlocal.h" |
| 27 | +#include "source_lcao/module_operator_lcao/ekinetic.h" |
| 28 | +#include "source_lcao/module_operator_lcao/overlap.h" |
| 29 | +#include "source_lcao/module_operator_lcao/td_pot_hybrid.h" |
| 30 | +#include "source_lcao/pulay_fs.h" |
| 31 | +#include "source_lcao/module_rt/force_rt_overlap.h" |
| 32 | + |
| 33 | + |
| 34 | +// vlocal, hartree, ewald, core correction, exchange-correlation terms in stress |
| 35 | +template <typename T> |
| 36 | +void Force_Stress_LCAO<T>::calStressPwPart(UnitCell& ucell, |
| 37 | + ModuleBase::matrix& sigmadvl, |
| 38 | + ModuleBase::matrix& sigmahar, |
| 39 | + ModuleBase::matrix& sigmaewa, |
| 40 | + ModuleBase::matrix& sigmacc, |
| 41 | + ModuleBase::matrix& sigmaxc, |
| 42 | + const double& etxc, |
| 43 | + const Charge* const chr, |
| 44 | + ModulePW::PW_Basis* rhopw, |
| 45 | + const pseudopot_cell_vl& locpp, |
| 46 | + const Structure_Factor& sf) |
| 47 | +{ |
| 48 | + ModuleBase::TITLE("Force_Stress_LCAO", "calStressPwPart"); |
| 49 | + |
| 50 | + // local pseudopotential stress: |
| 51 | + sc_pw.stress_loc(ucell, sigmadvl, rhopw, locpp.vloc, &sf, 0, chr); |
| 52 | + |
| 53 | + // hartree term |
| 54 | + sc_pw.stress_har(ucell, sigmahar, rhopw, 0, chr); |
| 55 | + |
| 56 | + // ewald stress: use plane wave only. |
| 57 | + sc_pw.stress_ewa(ucell, sigmaewa, rhopw, 0); // remain problem |
| 58 | + |
| 59 | + // stress due to core correlation. |
| 60 | + sc_pw.stress_cc(sigmacc, rhopw, ucell, &sf, 0, locpp.numeric, chr); |
| 61 | + |
| 62 | + // stress due to self-consistent charge. |
| 63 | + for (int i = 0; i < 3; i++) |
| 64 | + { |
| 65 | + sigmaxc(i, i) = -etxc / ucell.omega; |
| 66 | + } |
| 67 | + // Exchange-correlation for PBE |
| 68 | + sc_pw.stress_gga(ucell, sigmaxc, rhopw, chr); |
| 69 | + |
| 70 | + return; |
| 71 | +} |
| 72 | + |
| 73 | +#include "source_base/mathzone.h" |
| 74 | +// do symmetry for total force |
| 75 | +template <typename T> |
| 76 | +void Force_Stress_LCAO<T>::forceSymmetry(const UnitCell& ucell, ModuleBase::matrix& fcs, ModuleSymmetry::Symmetry* symm) |
| 77 | +{ |
| 78 | + double d1, d2, d3; |
| 79 | + for (int iat = 0; iat < ucell.nat; iat++) |
| 80 | + { |
| 81 | + ModuleBase::Mathzone::Cartesian_to_Direct(fcs(iat, 0), fcs(iat, 1), fcs(iat, 2), |
| 82 | + ucell.a1.x, ucell.a1.y, ucell.a1.z, ucell.a2.x, ucell.a2.y, ucell.a2.z, |
| 83 | + ucell.a3.x, ucell.a3.y, ucell.a3.z, d1, d2, d3); |
| 84 | + |
| 85 | + fcs(iat, 0) = d1; |
| 86 | + fcs(iat, 1) = d2; |
| 87 | + fcs(iat, 2) = d3; |
| 88 | + } |
| 89 | + symm->symmetrize_vec3_nat(fcs.c); |
| 90 | + for (int iat = 0; iat < ucell.nat; iat++) |
| 91 | + { |
| 92 | + ModuleBase::Mathzone::Direct_to_Cartesian(fcs(iat, 0), fcs(iat, 1), fcs(iat, 2), |
| 93 | + ucell.a1.x, ucell.a1.y, ucell.a1.z, ucell.a2.x, ucell.a2.y, ucell.a2.z, |
| 94 | + ucell.a3.x, ucell.a3.y, ucell.a3.z, d1, d2, d3); |
| 95 | + |
| 96 | + fcs(iat, 0) = d1; |
| 97 | + fcs(iat, 1) = d2; |
| 98 | + fcs(iat, 2) = d3; |
| 99 | + } |
| 100 | + return; |
| 101 | +} |
| 102 | + |
| 103 | +template class Force_Stress_LCAO<double>; |
| 104 | +template class Force_Stress_LCAO<std::complex<double>>; |
0 commit comments