From e5520064c4e0d25c4e1b515c51687e6c0aef2253 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Wed, 15 Oct 2025 16:08:29 +0800 Subject: [PATCH 1/3] fix relax_method parameter bug --- source/source_io/read_input_item_relax.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/source/source_io/read_input_item_relax.cpp b/source/source_io/read_input_item_relax.cpp index 59c6ab187a..f8c3453893 100644 --- a/source/source_io/read_input_item_relax.cpp +++ b/source/source_io/read_input_item_relax.cpp @@ -31,8 +31,10 @@ void ReadInput::item_relax() ModuleBase::WARNING_QUIT("ReadInput", warningstr); } }; + sync_stringvec(input.relax_method, para.input.relax_method.size(), ""); this->add_item(item); + // Input_Item item("relax_method"); // item.annotation = "cg; bfgs; sd; cg; cg_bfgs;"; // read_sync_string(input.relax_method); From cd77ad9a4bd4c223089e71164d7755a2e56ea818 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Fri, 17 Oct 2025 17:52:24 +0800 Subject: [PATCH 2/3] remove vector in bfgs and remove const number --- .../module_neighbor/test/prepare_unitcell.h | 2 +- source/source_cell/read_atom_species.cpp | 2 +- source/source_cell/read_atoms.cpp | 10 +-- source/source_cell/test/prepare_unitcell.h | 2 +- source/source_cell/update_cell.cpp | 2 +- .../charge_mixing_preconditioner.cpp | 4 +- .../module_dm/test/prepare_unitcell.h | 2 +- .../source_estate/test/charge_mixing_test.cpp | 12 +-- source/source_estate/test/prepare_unitcell.h | 2 +- .../module_surchem/test/setcell.h | 2 +- source/source_io/output_log.cpp | 2 +- source/source_io/test/prepare_unitcell.h | 2 +- .../source_io/test_serial/prepare_unitcell.h | 2 +- source/source_io/write_dipole.cpp | 6 +- source/source_lcao/FORCE_STRESS.cpp | 2 +- .../test/lambda_loop_helper_test.cpp | 2 +- .../module_deltaspin/test/prepare_unitcell.h | 2 +- .../module_hcontainer/test/prepare_unitcell.h | 2 +- source/source_md/test/setcell.h | 2 +- source/source_relax/bfgs.cpp | 89 +++++++++---------- source/source_relax/bfgs.h | 19 ++-- source/source_relax/ions_move_basic.cpp | 6 +- source/source_relax/ions_move_cg.cpp | 2 +- source/source_relax/lbfgs.cpp | 41 +++++---- source/source_relax/lbfgs.h | 34 +++---- source/source_relax/matrix_methods.cpp | 22 ++--- source/source_relax/matrix_methods.h | 6 +- 27 files changed, 141 insertions(+), 140 deletions(-) diff --git a/source/source_cell/module_neighbor/test/prepare_unitcell.h b/source/source_cell/module_neighbor/test/prepare_unitcell.h index 7af6d01e90..11bd94ba7b 100644 --- a/source/source_cell/module_neighbor/test/prepare_unitcell.h +++ b/source/source_cell/module_neighbor/test/prepare_unitcell.h @@ -98,7 +98,7 @@ class UcellTestPrepare } //lattice info ucell->lat0 = this->lat0; - ucell->lat0_angstrom = ucell->lat0 * 0.529177; + ucell->lat0_angstrom = ucell->lat0 * ModuleBase::BOHR_TO_A; ucell->tpiba = ModuleBase::TWO_PI/ucell->lat0; ucell->tpiba2 = ucell->tpiba * ucell->tpiba; ucell->latvec.e11 = this->latvec[0]; diff --git a/source/source_cell/read_atom_species.cpp b/source/source_cell/read_atom_species.cpp index 354d40cf8f..b8cf4c5744 100644 --- a/source/source_cell/read_atom_species.cpp +++ b/source/source_cell/read_atom_species.cpp @@ -150,7 +150,7 @@ bool read_lattice_constant(std::ifstream& ifa, { ModuleBase::WARNING_QUIT("read_atom_species","Lattice constant <= 0.0"); } - lat0_angstrom = lat0 * 0.529177; + lat0_angstrom = lat0 * ModuleBase::BOHR_TO_A; ModuleBase::GlobalFunc::OUT(ofs_running,"Lattice constant (Bohr)",lat0); ModuleBase::GlobalFunc::OUT(ofs_running,"Lattice constant (Angstrom)",lat0_angstrom); lat.tpiba = ModuleBase::TWO_PI / lat0; diff --git a/source/source_cell/read_atoms.cpp b/source/source_cell/read_atoms.cpp index 5162e6134e..1c1816c392 100644 --- a/source/source_cell/read_atoms.cpp +++ b/source/source_cell/read_atoms.cpp @@ -414,7 +414,7 @@ bool unitcell::read_atom_positions(UnitCell& ucell, } else if(Coordinate=="Cartesian_angstrom") { - ucell.atoms[it].tau[ia] = v / 0.529177 / ucell.lat0; + ucell.atoms[it].tau[ia] = v / ModuleBase::BOHR_TO_A / ucell.lat0; } else if(Coordinate=="Cartesian_angstrom_center_xy") { @@ -422,7 +422,7 @@ bool unitcell::read_atom_positions(UnitCell& ucell, ucell.latcenter.x = (ucell.latvec.e11 + ucell.latvec.e21 + ucell.latvec.e31)/2.0; ucell.latcenter.y = (ucell.latvec.e12 + ucell.latvec.e22 + ucell.latvec.e32)/2.0; ucell.latcenter.z = 0.0; - ucell.atoms[it].tau[ia] = v / 0.529177 / ucell.lat0 + ucell.latcenter; + ucell.atoms[it].tau[ia] = v / ModuleBase::BOHR_TO_A / ucell.lat0 + ucell.latcenter; } else if(Coordinate=="Cartesian_angstrom_center_xz") { @@ -430,7 +430,7 @@ bool unitcell::read_atom_positions(UnitCell& ucell, ucell.latcenter.x = (ucell.latvec.e11 + ucell.latvec.e21 + ucell.latvec.e31)/2.0; ucell.latcenter.y = 0.0; ucell.latcenter.z = (ucell.latvec.e13 + ucell.latvec.e23 + ucell.latvec.e33)/2.0; - ucell.atoms[it].tau[ia] = v / 0.529177 / ucell.lat0 + ucell.latcenter; + ucell.atoms[it].tau[ia] = v / ModuleBase::BOHR_TO_A / ucell.lat0 + ucell.latcenter; } else if(Coordinate=="Cartesian_angstrom_center_yz") { @@ -438,7 +438,7 @@ bool unitcell::read_atom_positions(UnitCell& ucell, ucell.latcenter.x = 0.0; ucell.latcenter.y = (ucell.latvec.e12 + ucell.latvec.e22 + ucell.latvec.e32)/2.0; ucell.latcenter.z = (ucell.latvec.e13 + ucell.latvec.e23 + ucell.latvec.e33)/2.0; - ucell.atoms[it].tau[ia] = v / 0.529177 / ucell.lat0 + ucell.latcenter; + ucell.atoms[it].tau[ia] = v / ModuleBase::BOHR_TO_A / ucell.lat0 + ucell.latcenter; } else if(Coordinate=="Cartesian_angstrom_center_xyz") { @@ -446,7 +446,7 @@ bool unitcell::read_atom_positions(UnitCell& ucell, ucell.latcenter.x = (ucell.latvec.e11 + ucell.latvec.e21 + ucell.latvec.e31)/2.0; ucell.latcenter.y = (ucell.latvec.e12 + ucell.latvec.e22 + ucell.latvec.e32)/2.0; ucell.latcenter.z = (ucell.latvec.e13 + ucell.latvec.e23 + ucell.latvec.e33)/2.0; - ucell.atoms[it].tau[ia] = v / 0.529177 / ucell.lat0 + ucell.latcenter; + ucell.atoms[it].tau[ia] = v / ModuleBase::BOHR_TO_A / ucell.lat0 + ucell.latcenter; } else if(Coordinate=="Cartesian_au") { diff --git a/source/source_cell/test/prepare_unitcell.h b/source/source_cell/test/prepare_unitcell.h index 32142daad9..7391464408 100644 --- a/source/source_cell/test/prepare_unitcell.h +++ b/source/source_cell/test/prepare_unitcell.h @@ -97,7 +97,7 @@ class UcellTestPrepare } //lattice info ucell->lat0 = this->lat0; - ucell->lat0_angstrom = ucell->lat0 * 0.529177; + ucell->lat0_angstrom = ucell->lat0 * ModuleBase::BOHR_TO_A; ucell->tpiba = ModuleBase::TWO_PI/ucell->lat0; ucell->tpiba2 = ucell->tpiba * ucell->tpiba; ucell->latvec.e11 = this->latvec[0]; diff --git a/source/source_cell/update_cell.cpp b/source/source_cell/update_cell.cpp index 187f962fd0..6a5c487581 100644 --- a/source/source_cell/update_cell.cpp +++ b/source/source_cell/update_cell.cpp @@ -319,7 +319,7 @@ void setup_cell_after_vc(UnitCell& ucell, std::ofstream& log) ucell.omega * pow(ModuleBase::BOHR_TO_A, 3)); } - ucell.lat0_angstrom = ucell.lat0 * 0.529177; + ucell.lat0_angstrom = ucell.lat0 * ModuleBase::BOHR_TO_A; ucell.tpiba = ModuleBase::TWO_PI / ucell.lat0; ucell.tpiba2 = ucell.tpiba * ucell.tpiba; diff --git a/source/source_estate/module_charge/charge_mixing_preconditioner.cpp b/source/source_estate/module_charge/charge_mixing_preconditioner.cpp index 00b5d42038..333cfb59e4 100644 --- a/source/source_estate/module_charge/charge_mixing_preconditioner.cpp +++ b/source/source_estate/module_charge/charge_mixing_preconditioner.cpp @@ -56,7 +56,7 @@ void Charge_Mixing::Kerker_screen_recip(std::complex* drhog) amin = this->mixing_beta; } - gg0 = std::pow(fac * 0.529177 / *this->tpiba, 2); + gg0 = std::pow(fac * ModuleBase::BOHR_TO_A / *this->tpiba, 2); const double gg0_amin = this->mixing_gg0_min / amin; @@ -138,7 +138,7 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) amin = this->mixing_beta; } - gg0 = std::pow(fac * 0.529177 / *this->tpiba, 2); + gg0 = std::pow(fac * ModuleBase::BOHR_TO_A / *this->tpiba, 2); const int is_idx = is * this->rhopw->npw; const double gg0_amin = this->mixing_gg0_min / amin; diff --git a/source/source_estate/module_dm/test/prepare_unitcell.h b/source/source_estate/module_dm/test/prepare_unitcell.h index e760a8da0e..0cbbd28905 100644 --- a/source/source_estate/module_dm/test/prepare_unitcell.h +++ b/source/source_estate/module_dm/test/prepare_unitcell.h @@ -96,7 +96,7 @@ class UcellTestPrepare } // lattice info ucell.lat0 = this->lat0; - ucell.lat0_angstrom = ucell.lat0 * 0.529177; + ucell.lat0_angstrom = ucell.lat0 * ModuleBase::BOHR_TO_A; ucell.tpiba = ModuleBase::TWO_PI / ucell.lat0; ucell.tpiba2 = ucell.tpiba * ucell.tpiba; ucell.latvec.e11 = this->latvec[0]; diff --git a/source/source_estate/test/charge_mixing_test.cpp b/source/source_estate/test/charge_mixing_test.cpp index 0fc96bd29e..a0e03b9911 100644 --- a/source/source_estate/test/charge_mixing_test.cpp +++ b/source/source_estate/test/charge_mixing_test.cpp @@ -619,7 +619,7 @@ TEST_F(ChargeMixingTest, KerkerScreenRecipTest) // kerker CMtest.mixing_gg0 = 1.0; CMtest.Kerker_screen_recip(drhog); - double gg0 = std::pow(0.529177, 2); + double gg0 = std::pow(ModuleBase::BOHR_TO_A, 2); for (int i = 0; i < pw_basis.npw; ++i) { double gg = this->pw_basis.gg[i]; @@ -650,7 +650,7 @@ TEST_F(ChargeMixingTest, KerkerScreenRecipTest) // mixing_gg0 = 1.0, mixing_gg0_mag = 0.0 CMtest.mixing_gg0 = 1.0; CMtest.Kerker_screen_recip(drhog); - gg0 = std::pow(0.529177, 2); + gg0 = std::pow(ModuleBase::BOHR_TO_A, 2); for (int i = 0; i < pw_basis.npw; ++i) { double gg = this->pw_basis.gg[i]; @@ -683,7 +683,7 @@ TEST_F(ChargeMixingTest, KerkerScreenRecipTest) // mixing_gg0 = 1.0, mixing_gg0_mag = 0.0 CMtest.mixing_gg0 = 1.0; CMtest.Kerker_screen_recip(drhog); - gg0 = std::pow(0.529177, 2); + gg0 = std::pow(ModuleBase::BOHR_TO_A, 2); for (int i = 0; i < pw_basis.npw; ++i) { double gg = this->pw_basis.gg[i]; @@ -701,8 +701,8 @@ TEST_F(ChargeMixingTest, KerkerScreenRecipTest) CMtest.mixing_gg0 = 1.0; CMtest.mixing_gg0_mag = 2.0; CMtest.Kerker_screen_recip(drhog); - double gg1 = std::pow(1.0 * 0.529177, 2); - double gg2 = std::pow(2.0 * 0.529177, 2); + double gg1 = std::pow(1.0 * ModuleBase::BOHR_TO_A, 2); + double gg2 = std::pow(2.0 * ModuleBase::BOHR_TO_A, 2); for (int i = 0; i < pw_basis.npw; ++i) { double gg = this->pw_basis.gg[i]; @@ -782,7 +782,7 @@ TEST_F(ChargeMixingTest, KerkerScreenRealTest) CMtest.mixing_gg0 = 1.0; PARAM.input.mixing_gg0_mag = 0.0; CMtest.Kerker_screen_recip(drhog); - const double gg0 = std::pow(0.529177, 2); + const double gg0 = std::pow(ModuleBase::BOHR_TO_A, 2); for (int i = 0; i < pw_basis.npw; ++i) { std::complex ration = drhog[i] / drhog[i+pw_basis.npw]; diff --git a/source/source_estate/test/prepare_unitcell.h b/source/source_estate/test/prepare_unitcell.h index 49635e89f3..eec2fa9a52 100644 --- a/source/source_estate/test/prepare_unitcell.h +++ b/source/source_estate/test/prepare_unitcell.h @@ -82,7 +82,7 @@ class UcellTestPrepare } //lattice info ucell->lat0 = this->lat0; - ucell->lat0_angstrom = ucell->lat0 * 0.529177; + ucell->lat0_angstrom = ucell->lat0 * ModuleBase::BOHR_TO_A; ucell->tpiba = ModuleBase::TWO_PI/ucell->lat0; ucell->tpiba2 = ucell->tpiba * ucell->tpiba; ucell->latvec.e11 = this->latvec[0]; diff --git a/source/source_hamilt/module_surchem/test/setcell.h b/source/source_hamilt/module_surchem/test/setcell.h index 6e132f7dd1..ebd38bcc09 100644 --- a/source/source_hamilt/module_surchem/test/setcell.h +++ b/source/source_hamilt/module_surchem/test/setcell.h @@ -50,7 +50,7 @@ class Setcell ucell.atoms[1].ncpp.psd = "O"; ucell.lat0 = 1; - ucell.lat0_angstrom = ucell.lat0 * 0.529177; + ucell.lat0_angstrom = ucell.lat0 * ModuleBase::BOHR_TO_A; ucell.tpiba = ModuleBase::TWO_PI / ucell.lat0; ucell.tpiba2 = ucell.tpiba * ucell.tpiba; diff --git a/source/source_io/output_log.cpp b/source/source_io/output_log.cpp index 4f73802e21..7a4471b0a6 100644 --- a/source/source_io/output_log.cpp +++ b/source/source_io/output_log.cpp @@ -224,7 +224,7 @@ void print_force(std::ofstream& ofs_running, double fac = 1.0; if (!ry) { - fac = ModuleBase::Ry_to_eV / 0.529177; + fac = ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A; } std::vector atom_label; diff --git a/source/source_io/test/prepare_unitcell.h b/source/source_io/test/prepare_unitcell.h index e05792a96d..af3cffb383 100644 --- a/source/source_io/test/prepare_unitcell.h +++ b/source/source_io/test/prepare_unitcell.h @@ -98,7 +98,7 @@ class UcellTestPrepare } //lattice info ucell->lat0 = this->lat0; - ucell->lat0_angstrom = ucell->lat0 * 0.529177; + ucell->lat0_angstrom = ucell->lat0 * ModuleBase::BOHR_TO_A; ucell->tpiba = ModuleBase::TWO_PI/ucell->lat0; ucell->tpiba2 = ucell->tpiba * ucell->tpiba; ucell->latvec.e11 = this->latvec[0]; diff --git a/source/source_io/test_serial/prepare_unitcell.h b/source/source_io/test_serial/prepare_unitcell.h index 0ae5df18d3..7e73e71892 100644 --- a/source/source_io/test_serial/prepare_unitcell.h +++ b/source/source_io/test_serial/prepare_unitcell.h @@ -98,7 +98,7 @@ class UcellTestPrepare } //lattice info ucell->lat0 = this->lat0; - ucell->lat0_angstrom = ucell->lat0 * 0.529177; + ucell->lat0_angstrom = ucell->lat0 * ModuleBase::BOHR_TO_A; ucell->tpiba = ModuleBase::TWO_PI/ucell->lat0; ucell->tpiba2 = ucell->tpiba * ucell->tpiba; ucell->latvec.e11 = this->latvec[0]; diff --git a/source/source_io/write_dipole.cpp b/source/source_io/write_dipole.cpp index b10a754f44..d2f00631c0 100644 --- a/source/source_io/write_dipole.cpp +++ b/source/source_io/write_dipole.cpp @@ -37,9 +37,9 @@ void ModuleIO::write_dipole(const UnitCell& ucell, #ifndef __MPI double dipole_elec_x = 0.0, dipole_elec_y = 0.0, dipole_elec_z = 0.0; - double lat_factor_x = ucell.lat0 * 0.529177 / rhopw->nx; - double lat_factor_y = ucell.lat0 * 0.529177 / rhopw->ny; - double lat_factor_z = ucell.lat0 * 0.529177 / rhopw->nz; + double lat_factor_x = ucell.lat0 * ModuleBase::BOHR_TO_A / rhopw->nx; + double lat_factor_y = ucell.lat0 * ModuleBase::BOHR_TO_A / rhopw->ny; + double lat_factor_z = ucell.lat0 * ModuleBase::BOHR_TO_A / rhopw->nz; for (int k = 0; k < rhopw->nz; k++) { diff --git a/source/source_lcao/FORCE_STRESS.cpp b/source/source_lcao/FORCE_STRESS.cpp index ef455260f8..043c45a2fa 100644 --- a/source/source_lcao/FORCE_STRESS.cpp +++ b/source/source_lcao/FORCE_STRESS.cpp @@ -536,7 +536,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, GlobalV::ofs_running << " " << std::setw(8) << iat; for (int i = 0; i < 3; i++) { - if (std::abs(fcs(iat, i) * ModuleBase::Ry_to_eV / 0.529177) + if (std::abs(fcs(iat, i) * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A) < Force_Stress_LCAO::force_invalid_threshold_ev) { fcs(iat, i) = 0.0; diff --git a/source/source_lcao/module_deltaspin/test/lambda_loop_helper_test.cpp b/source/source_lcao/module_deltaspin/test/lambda_loop_helper_test.cpp index 2acf2d34eb..445832e924 100644 --- a/source/source_lcao/module_deltaspin/test/lambda_loop_helper_test.cpp +++ b/source/source_lcao/module_deltaspin/test/lambda_loop_helper_test.cpp @@ -83,7 +83,7 @@ TEST_F(spinconstrain::SpinConstrainTest, CheckRestriction) std::vector> search = { {0.0, 0.0, 40} }; - double alpha_trial = 0.1 / 13.605698; + double alpha_trial = 0.1 / ModuleBase::Ry_to_eV; testing::internal::CaptureStdout(); sc.check_restriction(search, alpha_trial); std::string output = testing::internal::GetCapturedStdout(); diff --git a/source/source_lcao/module_deltaspin/test/prepare_unitcell.h b/source/source_lcao/module_deltaspin/test/prepare_unitcell.h index 0b280cd2ac..7abe206dbf 100644 --- a/source/source_lcao/module_deltaspin/test/prepare_unitcell.h +++ b/source/source_lcao/module_deltaspin/test/prepare_unitcell.h @@ -95,7 +95,7 @@ class UcellTestPrepare } // lattice info ucell.lat0 = this->lat0; - ucell.lat0_angstrom = ucell.lat0 * 0.529177; + ucell.lat0_angstrom = ucell.lat0 * ModuleBase::BOHR_TO_A; ucell.tpiba = ModuleBase::TWO_PI / ucell.lat0; ucell.tpiba2 = ucell.tpiba * ucell.tpiba; ucell.latvec.e11 = this->latvec[0]; diff --git a/source/source_lcao/module_hcontainer/test/prepare_unitcell.h b/source/source_lcao/module_hcontainer/test/prepare_unitcell.h index 3f0f58cabd..3df427d83b 100644 --- a/source/source_lcao/module_hcontainer/test/prepare_unitcell.h +++ b/source/source_lcao/module_hcontainer/test/prepare_unitcell.h @@ -94,7 +94,7 @@ class UcellTestPrepare } // lattice info ucell.lat0 = this->lat0; - ucell.lat0_angstrom = ucell.lat0 * 0.529177; + ucell.lat0_angstrom = ucell.lat0 * ModuleBase::BOHR_TO_A; ucell.tpiba = ModuleBase::TWO_PI / ucell.lat0; ucell.tpiba2 = ucell.tpiba * ucell.tpiba; ucell.latvec.e11 = this->latvec[0]; diff --git a/source/source_md/test/setcell.h b/source/source_md/test/setcell.h index 4c633e53b9..3841b8501b 100644 --- a/source/source_md/test/setcell.h +++ b/source/source_md/test/setcell.h @@ -36,7 +36,7 @@ class Setcell ucell.atom_label[0] = "Ar"; ucell.lat0 = 1; - ucell.lat0_angstrom = ucell.lat0 * 0.529177; + ucell.lat0_angstrom = ucell.lat0 * ModuleBase::BOHR_TO_A; ucell.tpiba = ModuleBase::TWO_PI / ucell.lat0; ucell.tpiba2 = ucell.tpiba * ucell.tpiba; diff --git a/source/source_relax/bfgs.cpp b/source/source_relax/bfgs.cpp index 00f5de2cc4..d105dc2791 100644 --- a/source/source_relax/bfgs.cpp +++ b/source/source_relax/bfgs.cpp @@ -1,7 +1,6 @@ #include "bfgs.h" #include "source_pw/module_pwdft/global.h" #include "source_base/module_external/lapack_connector.h" -#include "source_base/matrix3.h" #include "source_io/module_parameter/parameter.h" #include "ions_move_basic.h" #include "source_cell/update_cell.h" @@ -23,13 +22,13 @@ void BFGS::allocate(const int _size) H[i][i] = alpha; } - pos = std::vector> (size, std::vector(3, 0.0)); + pos = std::vector> (size, ModuleBase::Vector3(0.0, 0.0, 0.0)); pos0 = std::vector(3*size, 0.0); - pos_taud = std::vector> (size, std::vector(3, 0.0)); + pos_taud = std::vector> (size, ModuleBase::Vector3(0.0, 0.0, 0.0)); pos_taud0 = std::vector(3*size, 0.0); - dpos = std::vector>(size, std::vector(3, 0.0)); + dpos = std::vector>(size, ModuleBase::Vector3(0.0, 0.0, 0.0)); force0 = std::vector(3*size, 0.0); - force = std::vector>(size, std::vector(3, 0.0)); + force = std::vector>(size, ModuleBase::Vector3(0.0, 0.0, 0.0)); steplength = std::vector(size, 0.0); } @@ -39,13 +38,16 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell) GetPos(ucell,pos); GetPostaud(ucell,pos_taud); ucell.ionic_position_updated = true; - assert(_force.nr == force.size() && _force.nc == force[0].size()); + assert(_force.nr == force.size() && _force.nc == 3); for(int i = 0; i < _force.nr; i++) { - for(int j=0;j<_force.nc;j++) - { - force[i][j]=_force(i,j)*ModuleBase::Ry_to_eV/ModuleBase::BOHR_TO_A; - } + force[i].x=_force(i,0)*ModuleBase::Ry_to_eV/ModuleBase::BOHR_TO_A; + force[i].y=_force(i,1)*ModuleBase::Ry_to_eV/ModuleBase::BOHR_TO_A; + force[i].z=_force(i,2)*ModuleBase::Ry_to_eV/ModuleBase::BOHR_TO_A; + // for(int j=0;j<_force.nc;j++) + // { + // force[i][j]=_force(i,j)*ModuleBase::Ry_to_eV/ModuleBase::BOHR_TO_A; + // } } int k=0; for(int i=0;iDetermineStep(steplength,dpos,maxstep); this->UpdatePos(ucell); this->CalculateLargestGrad(_force,ucell); - this->IsRestrain(dpos); + this->IsRestrain(); // print out geometry information during bfgs_trad relax unitcell::print_tau(ucell.atoms,ucell.Coordinate,ucell.ntype,ucell.lat0,GlobalV::ofs_running); } -void BFGS::GetPos(UnitCell& ucell,std::vector>& pos) +void BFGS::GetPos(UnitCell& ucell,std::vector>& pos) { assert(pos.size() == ucell.nat); int k=0; @@ -85,16 +87,16 @@ void BFGS::GetPos(UnitCell& ucell,std::vector>& pos) { for(int j=0;j>& pos_taud) + std::vector>& pos_taud) { assert(pos_taud.size() == ucell.nat); int k=0; @@ -102,21 +104,21 @@ void BFGS::GetPostaud(UnitCell& ucell, { for(int j=0;j>& force, - std::vector>& pos, +void BFGS::PrepareStep(std::vector>& force, + std::vector>& pos, std::vector>& H, std::vector& pos0, std::vector& force0, std::vector& steplength, - std::vector>& dpos, + std::vector>& dpos, int& size, UnitCell& ucell) { @@ -199,7 +201,7 @@ void BFGS::Update(std::vector& pos, //shortest_move=shortest_move*ModuleBase::BOHR_TO_A*ucell.lat0; dpos[i]=shortest_move; } - std::vector> c=ReshapeVToM(dpos); + std::vector> c=ReshapeVToM(dpos); for(int iat=0; iat& pos, //convert unit ModuleBase::Vector3 move_ion_cart; - move_ion_cart.x = c[iat][0] *ModuleBase::BOHR_TO_A * ucell.lat0; - move_ion_cart.y = c[iat][1] * ModuleBase::BOHR_TO_A * ucell.lat0; - move_ion_cart.z = c[iat][2] * ModuleBase::BOHR_TO_A * ucell.lat0; + move_ion_cart.x = c[iat].x * ModuleBase::BOHR_TO_A * ucell.lat0; + move_ion_cart.y = c[iat].y * ModuleBase::BOHR_TO_A * ucell.lat0; + move_ion_cart.z = c[iat].z * ModuleBase::BOHR_TO_A * ucell.lat0; //convert pos ModuleBase::Vector3 move_ion_dr = move_ion_cart* ucell.latvec; @@ -247,7 +249,7 @@ void BFGS::Update(std::vector& pos, } void BFGS::DetermineStep(std::vector& steplength, - std::vector>& dpos, + std::vector>& dpos, double& maxstep) { std::vector::iterator maxsteplength = max_element(steplength.begin(), steplength.end()); @@ -259,10 +261,9 @@ void BFGS::DetermineStep(std::vector& steplength, double scale = maxstep / a; for(int i = 0; i < size; i++) { - for(int j=0;j<3;j++) - { - dpos[i][j]*=scale; - } + dpos[i].x*=scale; + dpos[i].y*=scale; + dpos[i].z*=scale; } } } @@ -272,11 +273,9 @@ void BFGS::UpdatePos(UnitCell& ucell) double a[3*size]; for(int i=0;iMAddM(pos, dpos);*/ } -void BFGS::IsRestrain(std::vector>& dpos) -{ - Ions_Move_Basic::converged = largest_grad - * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A grad= std::vector(3*size, 0.0); @@ -362,3 +355,9 @@ void BFGS::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell << std::endl; } } + +void BFGS::IsRestrain() +{ + Ions_Move_Basic::converged = largest_grad + * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A steplength;//the length of atoms displacement std::vector> H;//Hessian matrix std::vector force0;//force in previous step - std::vector> force; + std::vector> force; std::vector pos0;//atom pos in previous step(cartesian coordinates) - std::vector> pos; + std::vector> pos; std::vector pos_taud0;//atom pos in previous step(relative coordinates) - std::vector> pos_taud; - std::vector> dpos; + std::vector> pos_taud; + std::vector> dpos; - void PrepareStep(std::vector>& force,std::vector>& pos,std::vector>& H,std::vector& pos0,std::vector& force0,std::vector& steplength,std::vector>& dpos,int& size,UnitCell& ucell);//calculate the atomic displacement in one iteration step - void IsRestrain(std::vector>& dpos);//check if converged + void PrepareStep(std::vector>& force,std::vector>& pos,std::vector>& H,std::vector& pos0,std::vector& force0,std::vector& steplength,std::vector>& dpos,int& size,UnitCell& ucell);//calculate the atomic displacement in one iteration step + void IsRestrain();//check if converged void CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell); - void GetPos(UnitCell& ucell,std::vector>& pos); - void GetPostaud(UnitCell& ucell,std::vector>& pos_taud); + void GetPos(UnitCell& ucell,std::vector>& pos); + void GetPostaud(UnitCell& ucell,std::vector>& pos_taud); void Update(std::vector& pos, std::vector& force,std::vector>& H,UnitCell& ucell);//update hessian matrix - void DetermineStep(std::vector& steplength,std::vector>& dpos,double& maxstep);//normalize large atomic displacements based on maxstep + void DetermineStep(std::vector& steplength,std::vector>& dpos,double& maxstep);//normalize large atomic displacements based on maxstep void UpdatePos(UnitCell& ucell);//update ucell with the new coordinates }; diff --git a/source/source_relax/ions_move_basic.cpp b/source/source_relax/ions_move_basic.cpp index 1bf1e025a6..e4a2caef9d 100644 --- a/source/source_relax/ions_move_basic.cpp +++ b/source/source_relax/ions_move_basic.cpp @@ -147,10 +147,10 @@ void Ions_Move_Basic::check_converged(const UnitCell &ucell, const double *grad) if (PARAM.inp.out_level == "ie") { std::cout << " ETOT DIFF (eV) : " << Ions_Move_Basic::ediff * ModuleBase::Ry_to_eV << std::endl; - std::cout << " LARGEST GRAD (eV/Angstrom) : " << Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.529177 + std::cout << " LARGEST GRAD (eV/Angstrom) : " << Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << std::endl; - GlobalV::ofs_running << "\n Largest force is " << largest_grad * ModuleBase::Ry_to_eV / 0.529177 + GlobalV::ofs_running << "\n Largest force is " << largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << " eV/Angstrom while threshold is " << PARAM.inp.force_thr_ev << " eV/Angstrom" << std::endl; } @@ -178,7 +178,7 @@ void Ions_Move_Basic::check_converged(const UnitCell &ucell, const double *grad) else { GlobalV::ofs_running << "\n Ion relaxation is not converged yet (threshold is " - << PARAM.inp.force_thr * ModuleBase::Ry_to_eV / 0.529177 << ")" << std::endl; + << PARAM.inp.force_thr * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << ")" << std::endl; // std::cout << "\n etot_diff=" << etot_diff << " etot_thr=" << etot_thr //<< " largest_grad=" << largest_grad << " force_thr=" << PARAM.inp.force_thr << std::endl; Ions_Move_Basic::converged = false; diff --git a/source/source_relax/ions_move_cg.cpp b/source/source_relax/ions_move_cg.cpp index 9b8aa4b5c0..eec0009764 100644 --- a/source/source_relax/ions_move_cg.cpp +++ b/source/source_relax/ions_move_cg.cpp @@ -165,7 +165,7 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const if (Ions_Move_Basic::relax_method[0] == "cg_bfgs") { - if (Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.529177 + if (Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A < RELAX_CG_THR) // cg to bfgs by pengfei 13-8-8 { Ions_Move_Basic::relax_method[0] = "bfgs"; diff --git a/source/source_relax/lbfgs.cpp b/source/source_relax/lbfgs.cpp index 9ad2497d30..31eaa2116a 100644 --- a/source/source_relax/lbfgs.cpp +++ b/source/source_relax/lbfgs.cpp @@ -1,6 +1,5 @@ #include "lbfgs.h" #include "source_pw/module_pwdft/global.h" -#include "source_base/matrix3.h" #include "source_io/module_parameter/parameter.h" #include "ions_move_basic.h" #include "source_cell/update_cell.h" @@ -15,13 +14,13 @@ void LBFGS::allocate(const int _size) // initialize H0、H、pos0、force0、for iteration=0; H = std::vector>(3*size, std::vector(3*size, 0.0)); H0=1/alpha; - pos = std::vector> (size, std::vector(3, 0.0)); + pos = std::vector> (size, ModuleBase::Vector3(0.0, 0.0, 0.0)); pos0 = std::vector(3*size, 0.0); - pos_taud = std::vector> (size, std::vector(3, 0.0)); + pos_taud = std::vector> (size, ModuleBase::Vector3(0.0, 0.0, 0.0)); pos_taud0 = std::vector(3*size, 0.0); - dpos = std::vector>(size, std::vector(3, 0.0)); + dpos = std::vector>(size, ModuleBase::Vector3(0.0, 0.0, 0.0)); force0 = std::vector(3*size, 0.0); - force = std::vector>(size, std::vector(3, 0.0)); + force = std::vector>(size, ModuleBase::Vector3(0.0, 0.0, 0.0)); steplength = std::vector(size, 0.0); //l_search.init_line_search(); } @@ -64,12 +63,12 @@ void LBFGS::relax_step(const ModuleBase::matrix _force,UnitCell& ucell,const dou this->determine_step(steplength,dpos,maxstep); this->update_pos(ucell); this->calculate_largest_grad(_force,ucell); - this->is_restrain(dpos); + this->is_restrain(); // mohan add 2025-06-22 unitcell::print_tau(ucell.atoms,ucell.Coordinate,ucell.ntype,ucell.lat0,GlobalV::ofs_running); } -void LBFGS::get_pos(UnitCell& ucell,std::vector>& pos) +void LBFGS::get_pos(UnitCell& ucell,std::vector>& pos) { int k=0; for(int i=0;i>& pos) } } -void LBFGS::get_pos_taud(UnitCell& ucell,std::vector>& pos_taud) +void LBFGS::get_pos_taud(UnitCell& ucell,std::vector>& pos_taud) { int k=0; for(int i=0;i>& pos_t } } -void LBFGS::prepare_step(std::vector>& force, - std::vector>& pos, +void LBFGS::prepare_step(std::vector>& force, + std::vector>& pos, std::vector>& H, std::vector& pos0, std::vector& force0, - std::vector>& dpos, + std::vector>& dpos, UnitCell& ucell, const double &etot) { @@ -130,7 +129,7 @@ void LBFGS::prepare_step(std::vector>& force, std::vector temp0=DotInVAndFloat(z,-1); dpos=ReshapeVToM(temp0); std::vector temp1=DotInVAndFloat(changedforce,-1); - std::vector> g=ReshapeVToM(temp1); + //std::vector> g=ReshapeVToM(temp1); energy=etot; //alpha_k=l_search.line_search(ucell,pos,g,energy,maxstep,size,dpos,pos,solver); //std::vector temp2=DotInVAndFloat(temp0,alpha_k); @@ -150,7 +149,7 @@ void LBFGS::prepare_step(std::vector>& force, pos_taud0=ReshapeMToV(pos_taud); force0 = changedforce; } -void LBFGS::update(std::vector>& pos_taud, +void LBFGS::update(std::vector>& pos_taud, std::vector& pos_taud0, std::vector& force, std::vector& force0, @@ -178,7 +177,7 @@ void LBFGS::update(std::vector>& pos_taud, } dpos[i]=shortest_move; } - std::vector> c=ReshapeVToM(dpos); + std::vector> c=ReshapeVToM(dpos); for(int iat=0; iat>& pos_taud, //convert unit ModuleBase::Vector3 move_ion_cart; - move_ion_cart.x = c[iat][0] *ModuleBase::BOHR_TO_A * ucell.lat0; - move_ion_cart.y = c[iat][1] * ModuleBase::BOHR_TO_A * ucell.lat0; - move_ion_cart.z = c[iat][2] * ModuleBase::BOHR_TO_A * ucell.lat0; + move_ion_cart.x = c[iat].x * ModuleBase::BOHR_TO_A * ucell.lat0; + move_ion_cart.y = c[iat].y * ModuleBase::BOHR_TO_A * ucell.lat0; + move_ion_cart.z = c[iat].z * ModuleBase::BOHR_TO_A * ucell.lat0; //convert pos ModuleBase::Vector3 move_ion_dr = move_ion_cart* ucell.latvec; @@ -222,7 +221,7 @@ void LBFGS::update(std::vector>& pos_taud, rho.erase(rho.begin()); } } -void LBFGS::determine_step(std::vector& steplength,std::vector>& dpos,double& maxstep) +void LBFGS::determine_step(std::vector& steplength,std::vector>& dpos,double& maxstep) { std::vector::iterator maxsteplength = max_element(steplength.begin(), steplength.end()); double a = *maxsteplength; @@ -252,9 +251,9 @@ void LBFGS::update_pos(UnitCell& ucell) unitcell::update_pos_tau(ucell.lat,a,ucell.ntype,ucell.nat,ucell.atoms); } -void LBFGS::is_restrain(std::vector>& dpos) +void LBFGS::is_restrain() { - Ions_Move_Basic::converged = Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.529177> H; ///< Inverse Hessian approximation - std::vector force0; ///< Previous step forces - std::vector> force; ///< Force history - std::vector pos0; ///< Previous positions - std::vector> pos; ///< Position history - std::vector pos_taud0; ///< Previous fractional positions - std::vector> pos_taud; ///< Fractional position history - std::vector> dpos; ///< Position displacements + std::vector steplength;//the length of atoms displacement + std::vector> H;//Hessian matrix + std::vector force0;//force in previous step + std::vector> force; + std::vector pos0;//atom pos in previous step(cartesian coordinates) + std::vector> pos; + std::vector pos_taud0;//atom pos in previous step(relative coordinates) + std::vector> pos_taud; + std::vector> dpos; std::vector> s; ///< Position difference vectors std::vector> y; ///< Force difference vectors std::vector rho; ///< Scalar products for L-BFGS update - std::vector steplength; ///< Step lengths for each atom /** * @brief Prepare optimization step parameters */ - void prepare_step(std::vector>& force, - std::vector>& pos, + void prepare_step(std::vector>& force, + std::vector>& pos, std::vector>& H, std::vector& pos0, std::vector& force0, - std::vector>& dpos, + std::vector>& dpos, UnitCell& ucell, const double &etot); @@ -78,7 +78,7 @@ class LBFGS * @brief Judge if the cell is restrain * @param dpos Position displacements to constrain */ - void is_restrain(std::vector>& dpos); + void is_restrain(); /** * @brief Calculate maximum gradient component @@ -94,7 +94,7 @@ class LBFGS * @param pos Output position vector */ void get_pos(UnitCell& ucell, - std::vector>& pos); + std::vector>& pos); /** * @brief Get fractional positions from unit cell @@ -102,7 +102,7 @@ class LBFGS * @param pos_taud Output fractional positions */ void get_pos_taud(UnitCell& ucell, - std::vector>& pos_taud); + std::vector>& pos_taud); /** * @brief Update L-BFGS history buffers @@ -117,7 +117,7 @@ class LBFGS * @param y Force differences buffer * @param rho Scalar products buffer */ - void update(std::vector>& pos_taud, + void update(std::vector>& pos_taud, std::vector& pos_taud0, std::vector& force, std::vector& force0, @@ -135,7 +135,7 @@ class LBFGS * @param maxstep Maximum allowed step length */ void determine_step(std::vector& steplength, - std::vector>& dpos, + std::vector>& dpos, double& maxstep); /** diff --git a/source/source_relax/matrix_methods.cpp b/source/source_relax/matrix_methods.cpp index 27f106ae56..96079d55e1 100644 --- a/source/source_relax/matrix_methods.cpp +++ b/source/source_relax/matrix_methods.cpp @@ -2,15 +2,18 @@ -std::vector ReshapeMToV(std::vector>& matrix) + +std::vector ReshapeMToV(std::vector>& matrix) { assert(!matrix.empty()); - assert(matrix[0].size() == 3); int size = matrix.size(); std::vector result; result.reserve(3*size); - for (const auto& row : matrix) { - result.insert(result.end(), row.begin(), row.end()); + for (const auto& v : matrix) + { + result.push_back(v.x); + result.push_back(v.y); + result.push_back(v.z); } return result; } @@ -42,16 +45,15 @@ std::vector VSubV(std::vector& a, std::vector& b) return result; } -std::vector> ReshapeVToM(std::vector& matrix) +std::vector> ReshapeVToM(std::vector& matrix) { assert(matrix.size() % 3 == 0); - std::vector> result = std::vector>(matrix.size() / 3, std::vector(3)); + std::vector> result = std::vector>(matrix.size() / 3, ModuleBase::Vector3(0.0, 0.0, 0.0)); for(int i = 0; i < result.size(); i++) { - for(int j = 0; j < 3; j++) - { - result[i][j] = matrix[i*3 + j]; - } + result[i].x = matrix[i*3 ]; + result[i].y = matrix[i*3 + 1]; + result[i].z = matrix[i*3 + 2]; } return result; } diff --git a/source/source_relax/matrix_methods.h b/source/source_relax/matrix_methods.h index c87979d52e..fb6bba42b9 100644 --- a/source/source_relax/matrix_methods.h +++ b/source/source_relax/matrix_methods.h @@ -3,14 +3,14 @@ #include #include +#include "source_base/vector3.h" - -std::vector ReshapeMToV(std::vector>& matrix); +std::vector ReshapeMToV(std::vector>& matrix); std::vector> MAddM(std::vector>& a, std::vector>& b); std::vector VSubV(std::vector& a, std::vector& b); std::vector VAddV(std::vector& a, std::vector& b); -std::vector> ReshapeVToM(std::vector& matrix); +std::vector> ReshapeVToM(std::vector& matrix); std::vector DotInMAndV1(std::vector>& matrix, std::vector& vec); std::vector DotInMAndV2(std::vector>& matrix, std::vector& vec); double DotInVAndV(std::vector& vec1, std::vector& vec2); From b0415800b250e7bacc4450db696b007fb6c808de Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Sat, 18 Oct 2025 12:47:17 +0800 Subject: [PATCH 3/3] remove vector in bfgs and remove const number --- source/source_relax/test/bfgs_test.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/source/source_relax/test/bfgs_test.cpp b/source/source_relax/test/bfgs_test.cpp index a9cc3efcf8..fb82b90ca0 100644 --- a/source/source_relax/test/bfgs_test.cpp +++ b/source/source_relax/test/bfgs_test.cpp @@ -68,9 +68,9 @@ TEST_F(BFGSTest, DetermineStepScaling) bfgs.allocate(size); std::vector steplength = {1.0, 0.1}; - std::vector> dpos = { - {1.0, 1.0, 1.0}, - {0.1, 0.1, 0.1} + std::vector> dpos = { + ModuleBase::Vector3(1.0, 1.0, 1.0), + ModuleBase::Vector3(0.1, 0.1, 0.1) }; double maxstep = 0.5; bfgs.DetermineStep(steplength, dpos, maxstep);