From 117cdd7775696c53035b95d4133b2cf05de76fdd Mon Sep 17 00:00:00 2001 From: SY Wang Date: Mon, 19 Oct 2026 01:13:08 +0800 Subject: [PATCH] Precommit: Prettify source_relax --- source/source_relax/bfgs_basic.cpp | 14 +- source/source_relax/bfgs_basic.h | 17 +- source/source_relax/cg_base.cpp | 35 +- source/source_relax/cg_base.h | 39 +- source/source_relax/ions_move_basic.cpp | 57 +- source/source_relax/ions_move_basic.h | 39 +- source/source_relax/ions_move_bfgs.cpp | 46 +- source/source_relax/ions_move_bfgs.h | 32 +- source/source_relax/ions_move_bfgs2.cpp | 284 ++++--- source/source_relax/ions_move_bfgs2.h | 72 +- source/source_relax/ions_move_cg.cpp | 33 +- source/source_relax/ions_move_cg.h | 17 +- source/source_relax/ions_move_lbfgs.cpp | 276 ++++--- source/source_relax/ions_move_lbfgs.h | 75 +- source/source_relax/ions_move_methods.cpp | 27 +- source/source_relax/ions_move_methods.h | 25 +- source/source_relax/ions_move_sd.cpp | 24 +- source/source_relax/ions_move_sd.h | 14 +- source/source_relax/lat_change_method.cpp | 10 +- source/source_relax/lat_change_method.h | 15 +- source/source_relax/lattice_change_basic.cpp | 31 +- source/source_relax/lattice_change_basic.h | 23 +- source/source_relax/lattice_change_cg.cpp | 16 +- source/source_relax/lattice_change_cg.h | 12 +- source/source_relax/line_search.cpp | 10 +- source/source_relax/line_search.h | 19 +- source/source_relax/matrix_methods.cpp | 83 +- source/source_relax/matrix_methods.h | 8 +- source/source_relax/relax_data.cpp | 3 +- source/source_relax/relax_data.h | 25 +- source/source_relax/relax_driver.cpp | 126 +-- source/source_relax/relax_driver.h | 56 +- source/source_relax/relax_nsync.cpp | 45 +- source/source_relax/relax_nsync.h | 27 +- source/source_relax/relax_stru_io.cpp | 10 +- source/source_relax/relax_stru_io.h | 112 +-- source/source_relax/relax_sync.cpp | 33 +- source/source_relax/relax_sync.h | 28 +- source/source_relax/socket_driver.cpp | 14 +- source/source_relax/socket_frame.cpp | 98 +-- source/source_relax/socket_frame.h | 10 +- source/source_relax/socket_frame_utils.cpp | 19 +- source/source_relax/socket_handlers.cpp | 89 +- source/source_relax/socket_handlers.h | 10 +- source/source_relax/socket_ipi.cpp | 3 +- source/source_relax/socket_utils.cpp | 24 +- source/source_relax/socket_utils.h | 17 +- source/source_relax/unittests/for_test.h | 20 +- .../unittests/mock_remake_cell.cpp | 92 +-- .../source_relax/unittests/mock_remake_cell.h | 21 +- source/source_relax/unittests/relax_test.h | 42 +- .../unittests/test_bfgs_basic.cpp | 13 +- .../source_relax/unittests/test_cg_base.cpp | 11 +- .../unittests/test_ions_move_basic.cpp | 54 +- .../unittests/test_ions_move_bfgs.cpp | 63 +- .../unittests/test_ions_move_bfgs2.cpp | 87 +- .../unittests/test_ions_move_cg.cpp | 16 +- .../unittests/test_ions_move_lbfgs.cpp | 12 +- .../unittests/test_ions_move_methods.cpp | 108 ++- .../unittests/test_ions_move_sd.cpp | 11 +- .../unittests/test_lat_change_method.cpp | 13 +- .../unittests/test_lattice_change_basic.cpp | 40 +- .../unittests/test_lattice_change_cg.cpp | 12 +- .../unittests/test_line_search.cpp | 5 +- .../unittests/test_matrix_methods.cpp | 1 - .../unittests/test_relax_history.cpp | 6 +- .../unittests/test_relax_nsync.cpp | 20 +- .../unittests/test_relax_stru_io.cpp | 10 +- .../unittests/test_relax_sync.cpp | 766 ++++++++++-------- .../unittests/test_socket_driver.cpp | 88 +- .../unittests/test_socket_frame.cpp | 135 ++- .../unittests/test_socket_ipi.cpp | 3 +- 72 files changed, 2041 insertions(+), 1710 deletions(-) diff --git a/source/source_relax/bfgs_basic.cpp b/source/source_relax/bfgs_basic.cpp index 7a939974611..db722705a13 100644 --- a/source/source_relax/bfgs_basic.cpp +++ b/source/source_relax/bfgs_basic.cpp @@ -1,9 +1,11 @@ #include "bfgs_basic.h" -#include + #include "ions_move_basic.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" +#include + using namespace Ions_Move_Basic; double BFGS_Basic::relax_bfgs_w1 = -1.0; // default is 0.01 @@ -31,7 +33,7 @@ void BFGS_Basic::allocate_basic(void) return; } -void BFGS_Basic::update_inverse_hessian(const double &lat0, std::ofstream& ofs) +void BFGS_Basic::update_inverse_hessian(const double& lat0, std::ofstream& ofs) { // ModuleBase::TITLE("Ions_Move_BFGS","update_inverse_hessian"); assert(dim > 0); @@ -157,7 +159,11 @@ void BFGS_Basic::save_bfgs(void) // a new bfgs step is done // we have already done well in the previous direction // we should get a new direction in this case -void BFGS_Basic::new_step(const double &lat0, int& update_iter, std::ofstream& ofs, std::vector& etot_info, const int test_relax_method) +void BFGS_Basic::new_step(const double& lat0, + int& update_iter, + std::ofstream& ofs, + std::vector& etot_info, + const int test_relax_method) { ModuleBase::TITLE("BFGS_Basic", "new_step"); @@ -341,7 +347,7 @@ void BFGS_Basic::compute_trust_radius(std::ofstream& ofs, std::vector& e return; } -double BFGS_Basic::check_move(const double &lat0, const double &pos, const double &pos_p) +double BFGS_Basic::check_move(const double& lat0, const double& pos, const double& pos_p) { // this must be careful. // unit is ucell.lat0. diff --git a/source/source_relax/bfgs_basic.h b/source/source_relax/bfgs_basic.h index 8b871b5e4ab..41892826bfb 100644 --- a/source/source_relax/bfgs_basic.h +++ b/source/source_relax/bfgs_basic.h @@ -1,9 +1,10 @@ #ifndef BFGS_BASIC #define BFGS_BASIC +#include "source_base/matrix.h" + #include #include -#include "source_base/matrix.h" #include // references @@ -129,7 +130,11 @@ class BFGS_Basic protected: void allocate_basic(void); - void new_step(const double& lat0, int& update_iter, std::ofstream& ofs, std::vector& etot_info, const int test_relax_method); + void new_step(const double& lat0, + int& update_iter, + std::ofstream& ofs, + std::vector& etot_info, + const int test_relax_method); void reset_hessian(void); void save_bfgs(void); @@ -146,15 +151,15 @@ class BFGS_Basic static double relax_bfgs_w2; // fixed: parameters for Wolfe conditions. protected: - bool save_flag=false; - bool tr_min_hit=false; //.TRUE. if the trust_radius has already been set - // to the minimum value at the previous step + bool save_flag = false; + bool tr_min_hit = false; //.TRUE. if the trust_radius has already been set + // to the minimum value at the previous step // mohan add 2010-07-27 double check_move(const double& lat0, const double& pos, const double& pos_p); private: - bool wolfe_flag=false; + bool wolfe_flag = false; ModuleBase::matrix inv_hess; int bfgs_ndim; diff --git a/source/source_relax/cg_base.cpp b/source/source_relax/cg_base.cpp index 9fee5024d95..c34119ebf8b 100644 --- a/source/source_relax/cg_base.cpp +++ b/source/source_relax/cg_base.cpp @@ -1,10 +1,16 @@ #include "cg_base.h" -#include + #include "source_base/global_function.h" -void CG_Base::setup_cg_grad(const int dim, double *grad, const double *grad0, - double *cg_grad, const double *cg_grad0, - const int &ncggrad, int &flag) +#include + +void CG_Base::setup_cg_grad(const int dim, + double* grad, + const double* grad0, + double* cg_grad, + const double* cg_grad0, + const int& ncggrad, + int& flag) { ModuleBase::TITLE("CG_Base", "setup_cg_grad"); double gamma = 0.0; @@ -32,7 +38,7 @@ void CG_Base::setup_cg_grad(const int dim, double *grad, const double *grad0, cgp_g += cg_grad0[i] * grad[i]; } - const double gamma1 = gg / gp_gp; // FR + const double gamma1 = gg / gp_gp; // FR const double gamma2 = (gg - g_gp) / gp_gp; // PRP if (gamma1 < 0.5) @@ -51,7 +57,7 @@ void CG_Base::setup_cg_grad(const int dim, double *grad, const double *grad0, } } -void CG_Base::setup_move(const int dim, double *move, double *cg_gradn, const double &trust_radius) +void CG_Base::setup_move(const int dim, double* move, double* cg_gradn, const double& trust_radius) { ModuleBase::TITLE("CG_Base", "setup_move"); for (int i = 0; i < dim; ++i) @@ -60,9 +66,7 @@ void CG_Base::setup_move(const int dim, double *move, double *cg_gradn, const do } } -void CG_Base::Brent(double &fa, double &fb, double &fc, - double &xa, double &xb, double &xc, - double &best_x, double &xpt) +void CG_Base::Brent(double& fa, double& fb, double& fc, double& xa, double& xb, double& xc, double& best_x, double& xpt) { ModuleBase::TITLE("CG_Base", "Brent"); double dmove = 0.0; @@ -129,7 +133,7 @@ void CG_Base::Brent(double &fa, double &fb, double &fc, xc = dmove; } -void CG_Base::f_cal(const int dim, const double *g0, const double *g1, double &f_value) +void CG_Base::f_cal(const int dim, const double* g0, const double* g1, double& f_value) { ModuleBase::TITLE("CG_Base", "f_cal"); double hv0 = 0.0; @@ -142,9 +146,12 @@ void CG_Base::f_cal(const int dim, const double *g0, const double *g1, double &f f_value = hel / sqrt(hv0); } -void CG_Base::third_order(const double &e0, const double &e1, - const double &fa, const double &fb, - const double x, double &best_x) +void CG_Base::third_order(const double& e0, + const double& e1, + const double& fa, + const double& fb, + const double x, + double& best_x) { ModuleBase::TITLE("CG_Base", "third_order"); double k3 = 0.0; @@ -193,7 +200,7 @@ void CG_Base::third_order(const double &e0, const double &e1, best_x = dmove; } -void CG_Base::normalize(const int dim, double *cg_gradn, const double *cg_grad) +void CG_Base::normalize(const int dim, double* cg_gradn, const double* cg_grad) { ModuleBase::TITLE("CG_Base", "normalize"); double norm = 0.0; diff --git a/source/source_relax/cg_base.h b/source/source_relax/cg_base.h index e2a3d929121..825f92b4bc6 100644 --- a/source/source_relax/cg_base.h +++ b/source/source_relax/cg_base.h @@ -4,13 +4,14 @@ /** * @class CG_Base * @brief Base class for conjugate gradient optimization algorithms. - * + * * This class provides common CG utility functions that can be shared between * different optimization implementations (ions movement and lattice change). * The functions are designed to be independent of specific data sources. */ -class CG_Base { -public: +class CG_Base +{ + public: /** * @brief Setup conjugate gradient direction. * @param dim Dimension of the optimization problem @@ -21,9 +22,13 @@ class CG_Base { * @param ncggrad CG iteration counter * @param flag Restart flag */ - void setup_cg_grad(const int dim, double *grad, const double *grad0, - double *cg_grad, const double *cg_grad0, - const int &ncggrad, int &flag); + void setup_cg_grad(const int dim, + double* grad, + const double* grad0, + double* cg_grad, + const double* cg_grad0, + const int& ncggrad, + int& flag); /** * @brief Brent's method for one-dimensional minimization. @@ -36,9 +41,7 @@ class CG_Base { * @param best_x Output best step length * @param xpt Previous best point */ - void Brent(double &fa, double &fb, double &fc, - double &xa, double &xb, double &xc, - double &best_x, double &xpt); + void Brent(double& fa, double& fb, double& fc, double& xa, double& xb, double& xc, double& best_x, double& xpt); /** * @brief Third-order polynomial interpolation for line search. @@ -49,9 +52,12 @@ class CG_Base { * @param x Current step length * @param best_x Output optimal step length */ - void third_order(const double &e0, const double &e1, - const double &fa, const double &fb, - const double x, double &best_x); + void third_order(const double& e0, + const double& e1, + const double& fa, + const double& fb, + const double x, + double& best_x); /** * @brief Calculate projection of gradient onto search direction. @@ -60,7 +66,7 @@ class CG_Base { * @param g1 Second vector (usually gradient) * @param f_value Output projection value */ - void f_cal(const int dim, const double *g0, const double *g1, double &f_value); + void f_cal(const int dim, const double* g0, const double* g1, double& f_value); /** * @brief Calculate movement vector from CG direction. @@ -69,17 +75,16 @@ class CG_Base { * @param cg_gradn Normalized CG direction * @param trust_radius Step size limit */ - void setup_move(const int dim, double *move, double *cg_gradn, const double &trust_radius); - -protected: + void setup_move(const int dim, double* move, double* cg_gradn, const double& trust_radius); + protected: /** * @brief Normalize vector. * @param dim Dimension of vector * @param cg_gradn Output normalized vector * @param cg_grad Input vector */ - void normalize(const int dim, double *cg_gradn, const double *cg_grad); + void normalize(const int dim, double* cg_gradn, const double* cg_grad); }; #endif \ No newline at end of file diff --git a/source/source_relax/ions_move_basic.cpp b/source/source_relax/ions_move_basic.cpp index 1e8764ae958..795770f85f3 100644 --- a/source/source_relax/ions_move_basic.cpp +++ b/source/source_relax/ions_move_basic.cpp @@ -1,10 +1,11 @@ #include "ions_move_basic.h" -#include #include "source_base/global_function.h" #include "source_base/global_variable.h" -#include "source_cell/update_cell.h" #include "source_cell/print_cell.h" +#include "source_cell/update_cell.h" + +#include // Ions-specific parameters (shared variables are in Relax_Data) double Ions_Move_Basic::trust_radius = 0.0; @@ -14,7 +15,11 @@ double Ions_Move_Basic::relax_bfgs_rmin = -1.0; // default is 1e-5 double Ions_Move_Basic::relax_bfgs_init = -1.0; // default is 0.5 double Ions_Move_Basic::best_xxx = 1.0; -void Ions_Move_Basic::setup_gradient(const UnitCell &ucell, const ModuleBase::matrix &force, double *pos, double *grad, std::ofstream& ofs) +void Ions_Move_Basic::setup_gradient(const UnitCell& ucell, + const ModuleBase::matrix& force, + double* pos, + double* grad, + std::ofstream& ofs) { ModuleBase::TITLE("Ions_Move_Basic", "setup_gradient"); @@ -33,7 +38,7 @@ void Ions_Move_Basic::setup_gradient(const UnitCell &ucell, const ModuleBase::ma int iat = 0; for (int it = 0; it < ucell.ntype; it++) { - Atom *atom = &ucell.atoms[it]; + Atom* atom = &ucell.atoms[it]; for (int ia = 0; ia < ucell.atoms[it].na; ia++) { for (int ik = 0; ik < 3; ++ik) @@ -51,7 +56,11 @@ void Ions_Move_Basic::setup_gradient(const UnitCell &ucell, const ModuleBase::ma return; } -void Ions_Move_Basic::move_atoms(UnitCell &ucell, double *move, double *pos, std::ofstream& ofs, const int test_relax_method) +void Ions_Move_Basic::move_atoms(UnitCell& ucell, + double* move, + double* pos, + std::ofstream& ofs, + const int test_relax_method) { ModuleBase::TITLE("Ions_Move_Basic", "move_atoms"); @@ -65,17 +74,16 @@ void Ions_Move_Basic::move_atoms(UnitCell &ucell, double *move, double *pos, std { int iat = 0; ofs << "\n movement of ions (unit is Bohr) : " << std::endl; - ofs << " " << std::setw(12) << "Atom" << std::setw(15) << "x" << std::setw(15) << "y" - << std::setw(15) << "z" << std::endl; + ofs << " " << std::setw(12) << "Atom" << std::setw(15) << "x" << std::setw(15) << "y" << std::setw(15) << "z" + << std::endl; for (int it = 0; it < ucell.ntype; it++) { for (int ia = 0; ia < ucell.atoms[it].na; ia++) { std::stringstream ss; ss << "move_" << ucell.atoms[it].label << ia + 1; - ofs << " " << std::setw(12) << ss.str().c_str() << std::setw(15) << move[3 * iat + 0] - << std::setw(15) << move[3 * iat + 1] << std::setw(15) << move[3 * iat + 2] - << std::endl; + ofs << " " << std::setw(12) << ss.str().c_str() << std::setw(15) << move[3 * iat + 0] << std::setw(15) + << move[3 * iat + 1] << std::setw(15) << move[3 * iat + 2] << std::endl; iat++; } } @@ -85,7 +93,7 @@ void Ions_Move_Basic::move_atoms(UnitCell &ucell, double *move, double *pos, std const double move_threshold = 1.0e-10; const int total_freedom = ucell.nat * 3; - if (ModuleSymmetry::Symmetry::symm_flag && ucell.symm.all_mbl && ucell.symm.nrotk > 0) + if (ModuleSymmetry::Symmetry::symm_flag && ucell.symm.all_mbl && ucell.symm.nrotk > 0) { ucell.symm.symmetrize_vec3_nat(move); } @@ -97,18 +105,18 @@ void Ions_Move_Basic::move_atoms(UnitCell &ucell, double *move, double *pos, std pos[i] += move[i]; } } - unitcell::update_pos_tau(ucell.lat,pos,ucell.ntype,ucell.nat,ucell.atoms); + unitcell::update_pos_tau(ucell.lat, pos, ucell.ntype, ucell.nat, ucell.atoms); //-------------------------------------------- // Print out the structure file. //-------------------------------------------- - unitcell::print_tau(ucell.atoms,ucell.Coordinate,ucell.ntype,ucell.lat0,ofs); + unitcell::print_tau(ucell.atoms, ucell.Coordinate, ucell.ntype, ucell.lat0, ofs); return; } -bool Ions_Move_Basic::check_converged(const UnitCell &ucell, - const double *grad, +bool Ions_Move_Basic::check_converged(const UnitCell& ucell, + const double* grad, int& update_iter, std::ofstream& ofs, std::vector& etot_info, @@ -147,9 +155,8 @@ bool Ions_Move_Basic::check_converged(const UnitCell &ucell, { const double ediff = etot_info[0] - etot_info[1]; std::cout << " ETOT DIFF (eV) : " << ediff * ModuleBase::Ry_to_eV << std::endl; - std::cout << " LARGEST GRAD (eV/Angstrom) : " - << Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A - << std::endl; + std::cout << " LARGEST GRAD (eV/Angstrom) : " + << Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << std::endl; } const double etot_diff = std::abs(etot_info[0] - etot_info[1]); @@ -162,7 +169,7 @@ bool Ions_Move_Basic::check_converged(const UnitCell &ucell, ofs << " it may converged, otherwise no movement of atom is allowed." << std::endl; return true; } - else if (etot_diff < etot_thr && Ions_Move_Basic::largest_grad < force_thr ) + else if (etot_diff < etot_thr && Ions_Move_Basic::largest_grad < force_thr) { ++update_iter; return true; @@ -173,7 +180,11 @@ bool Ions_Move_Basic::check_converged(const UnitCell &ucell, } } -void Ions_Move_Basic::terminate(const bool converged, const int update_iter, const UnitCell &ucell, const int istep, std::ofstream& ofs) +void Ions_Move_Basic::terminate(const bool converged, + const int update_iter, + const UnitCell& ucell, + const int istep, + std::ofstream& ofs) { ModuleBase::TITLE("Ions_Move_Basic", "terminate"); if (converged) @@ -204,11 +215,11 @@ void Ions_Move_Basic::terminate(const bool converged, const int update_iter, con //----------------------------------------------------------- // Print the structure. //----------------------------------------------------------- - unitcell::print_tau(ucell.atoms,ucell.Coordinate,ucell.ntype,ucell.lat0,ofs); + unitcell::print_tau(ucell.atoms, ucell.Coordinate, ucell.ntype, ucell.lat0, ofs); return; } -void Ions_Move_Basic::setup_etot(const double &energy_in, const int istep, std::vector& etot_info) +void Ions_Move_Basic::setup_etot(const double& energy_in, const int istep, std::vector& etot_info) { // etot_info[0] = etot (current total energy) // etot_info[1] = etot_p (previous total energy) @@ -226,7 +237,7 @@ void Ions_Move_Basic::setup_etot(const double &energy_in, const int istep, std:: } } -double Ions_Move_Basic::dot_func(const double *a, const double *b, const int &dim_in) +double Ions_Move_Basic::dot_func(const double* a, const double* b, const int& dim_in) { double result = 0.0; for (int i = 0; i < dim_in; i++) diff --git a/source/source_relax/ions_move_basic.h b/source/source_relax/ions_move_basic.h index d8cd9a4b839..20c72976698 100644 --- a/source/source_relax/ions_move_basic.h +++ b/source/source_relax/ions_move_basic.h @@ -1,16 +1,17 @@ #ifndef IONS_MOVE_BASIC_H #define IONS_MOVE_BASIC_H -#include -#include #include "relax_data.h" #include "source_base/matrix.h" #include "source_cell/unitcell.h" +#include +#include + /** * @namespace Ions_Move_Basic * @brief Basic utilities and shared state for ionic relaxation algorithms. - * + * * This namespace provides common functions and parameters used by all * ion movement methods (BFGS, CG, SD, etc.). It shares core state variables * through references to Relax_Data, ensuring consistent data across different @@ -18,16 +19,16 @@ */ namespace Ions_Move_Basic { -static int& dim = Relax_Data::dim; ///< Dimension of free variables (3 * number of atoms) +static int& dim = Relax_Data::dim; ///< Dimension of free variables (3 * number of atoms) static double& largest_grad = Relax_Data::largest_grad; ///< Largest gradient component // Ions-specific parameters (not shared with lattice change) -extern double trust_radius; ///< Current trust radius -extern double trust_radius_old; ///< Previous trust radius -extern double relax_bfgs_rmax; ///< Maximum trust radius (default: 0.8 Bohr) -extern double relax_bfgs_rmin; ///< Minimum trust radius (default: 1e-5 Bohr) -extern double relax_bfgs_init; ///< Initial trust radius (default: 0.5 Bohr) -extern double best_xxx; ///< Last step length from CG, used as BFGS initial guess +extern double trust_radius; ///< Current trust radius +extern double trust_radius_old; ///< Previous trust radius +extern double relax_bfgs_rmax; ///< Maximum trust radius (default: 0.8 Bohr) +extern double relax_bfgs_rmin; ///< Minimum trust radius (default: 1e-5 Bohr) +extern double relax_bfgs_init; ///< Initial trust radius (default: 0.5 Bohr) +extern double best_xxx; ///< Last step length from CG, used as BFGS initial guess /** * @brief Setup gradient from atomic forces. @@ -37,7 +38,11 @@ extern double best_xxx; ///< Last step length from CG, used as BFGS * @param grad Output gradient array (dimension: dim) * @param ofs Output stream for logging */ -void setup_gradient(const UnitCell &ucell, const ModuleBase::matrix &force, double *pos, double *grad, std::ofstream& ofs); +void setup_gradient(const UnitCell& ucell, + const ModuleBase::matrix& force, + double* pos, + double* grad, + std::ofstream& ofs); /** * @brief Move atoms according to displacement vector. @@ -47,7 +52,7 @@ void setup_gradient(const UnitCell &ucell, const ModuleBase::matrix &force, doub * @param ofs Output stream for logging * @param test_relax_method Verbosity level for relaxation debug output */ -void move_atoms(UnitCell &ucell, double *move, double *pos, std::ofstream& ofs, const int test_relax_method); +void move_atoms(UnitCell& ucell, double* move, double* pos, std::ofstream& ofs, const int test_relax_method); /** * @brief Check convergence based on gradient threshold. @@ -62,8 +67,8 @@ void move_atoms(UnitCell &ucell, double *move, double *pos, std::ofstream& ofs, * @param test_relax_method Verbosity level for relaxation debug output * @return true if converged, false otherwise */ -bool check_converged(const UnitCell &ucell, - const double *grad, +bool check_converged(const UnitCell& ucell, + const double* grad, int& update_iter, std::ofstream& ofs, std::vector& etot_info, @@ -80,7 +85,7 @@ bool check_converged(const UnitCell &ucell, * @param istep Current ionic step index * @param ofs Output stream for logging */ -void terminate(const bool converged, const int update_iter, const UnitCell &ucell, const int istep, std::ofstream& ofs); +void terminate(const bool converged, const int update_iter, const UnitCell& ucell, const int istep, std::ofstream& ofs); /** * @brief Update energy values and compute energy difference. @@ -89,7 +94,7 @@ void terminate(const bool converged, const int update_iter, const UnitCell &ucel * @param istep Current ionic step index * @param etot_info Energy information array [etot, etot_p, ediff] */ -void setup_etot(const double &energy_in, const int istep, std::vector& etot_info); +void setup_etot(const double& energy_in, const int istep, std::vector& etot_info); /** * @brief Compute dot product of two vectors. @@ -98,7 +103,7 @@ void setup_etot(const double &energy_in, const int istep, std::vector& e * @param dim_in Dimension of vectors * @return Dot product value */ -double dot_func(const double *a, const double *b, const int &dim_in); +double dot_func(const double* a, const double* b, const int& dim_in); /** * @brief Third-order polynomial interpolation for line search. diff --git a/source/source_relax/ions_move_bfgs.cpp b/source/source_relax/ions_move_bfgs.cpp index 7789c49bb93..87172425885 100644 --- a/source/source_relax/ions_move_bfgs.cpp +++ b/source/source_relax/ions_move_bfgs.cpp @@ -1,10 +1,11 @@ #include "ions_move_bfgs.h" -#include #include "ions_move_basic.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" +#include + //============= MAP OF BFGS =========================== // (1) start() -> BFGS_Basic::check_converged() // -> restart_bfgs() -> bfgs_routine() -> save_bfgs() @@ -18,12 +19,12 @@ Ions_Move_BFGS::Ions_Move_BFGS() init_done = false; } -Ions_Move_BFGS::~Ions_Move_BFGS(){}; +Ions_Move_BFGS::~Ions_Move_BFGS() {}; void Ions_Move_BFGS::allocate() { ModuleBase::TITLE("Ions_Move_BFGS", "init"); - if (init_done) + if (init_done) { return; } @@ -57,7 +58,14 @@ void Ions_Move_BFGS::reset() Ions_Move_Basic::trust_radius_old = 0.0; } -bool Ions_Move_BFGS::start(UnitCell& ucell, const ModuleBase::matrix& force, const double& energy_in, const int istep, int& update_iter, std::ofstream& ofs, std::vector& etot_info, const Relax_Criteria& criteria) +bool Ions_Move_BFGS::start(UnitCell& ucell, + const ModuleBase::matrix& force, + const double& energy_in, + const int istep, + int& update_iter, + std::ofstream& ofs, + std::vector& etot_info, + const Relax_Criteria& criteria) { ModuleBase::TITLE("Ions_Move_BFGS", "start"); @@ -72,7 +80,15 @@ bool Ions_Move_BFGS::start(UnitCell& ucell, const ModuleBase::matrix& force, con Ions_Move_Basic::setup_gradient(ucell, force, pos_tmp.data(), this->grad.data(), ofs); } Ions_Move_Basic::setup_etot(energy_in, istep, etot_info); - bool converged = Ions_Move_Basic::check_converged(ucell, this->grad.data(), update_iter, ofs, etot_info, criteria.force_thr, criteria.force_thr_ev, criteria.out_level, criteria.test_relax_method); + bool converged = Ions_Move_Basic::check_converged(ucell, + this->grad.data(), + update_iter, + ofs, + etot_info, + criteria.force_thr, + criteria.force_thr_ev, + criteria.out_level, + criteria.test_relax_method); if (converged) { @@ -82,7 +98,13 @@ bool Ions_Move_BFGS::start(UnitCell& ucell, const ModuleBase::matrix& force, con else { this->restart_bfgs(ucell.lat0, update_iter, ofs, criteria.test_relax_method); - this->bfgs_routine(ucell.lat0, istep, update_iter, ofs, etot_info, criteria.out_level, criteria.test_relax_method); + this->bfgs_routine(ucell.lat0, + istep, + update_iter, + ofs, + etot_info, + criteria.out_level, + criteria.test_relax_method); this->save_bfgs(); Ions_Move_Basic::move_atoms(ucell, move.data(), pos.data(), ofs, criteria.test_relax_method); @@ -130,8 +152,8 @@ void Ions_Move_BFGS::restart_bfgs(const double& lat0, int& update_iter, std::ofs { // bfgs initialization std::fill(pos_p.begin(), pos_p.end(), 0.0); - std::fill(grad_p.begin(), grad_p.end(), 0.0); - std::fill(move_p.begin(), move_p.end(), 0.0); + std::fill(grad_p.begin(), grad_p.end(), 0.0); + std::fill(move_p.begin(), move_p.end(), 0.0); update_iter = 0; @@ -166,7 +188,13 @@ void Ions_Move_BFGS::restart_bfgs(const double& lat0, int& update_iter, std::ofs return; } -void Ions_Move_BFGS::bfgs_routine(const double& lat0, const int istep, int& update_iter, std::ofstream& ofs, std::vector& etot_info, const std::string& out_level, const int test_relax_method) +void Ions_Move_BFGS::bfgs_routine(const double& lat0, + const int istep, + int& update_iter, + std::ofstream& ofs, + std::vector& etot_info, + const std::string& out_level, + const int test_relax_method) { ModuleBase::TITLE("Ions_Move_BFGS", "bfgs_routine"); using namespace Ions_Move_Basic; diff --git a/source/source_relax/ions_move_bfgs.h b/source/source_relax/ions_move_bfgs.h index 9f0a85c8b7d..fd3eeb33385 100644 --- a/source/source_relax/ions_move_bfgs.h +++ b/source/source_relax/ions_move_bfgs.h @@ -1,13 +1,14 @@ #ifndef IONS_MOVE_BFGS_H #define IONS_MOVE_BFGS_H -#include -#include -#include #include "bfgs_basic.h" #include "relax_criteria.h" #include "source_base/matrix.h" #include "source_cell/unitcell.h" + +#include +#include +#include class Ions_Move_BFGS : public BFGS_Basic { public: @@ -16,7 +17,14 @@ class Ions_Move_BFGS : public BFGS_Basic void allocate(void); void reset(void); - bool start(UnitCell& ucell, const ModuleBase::matrix& force, const double& energy_in, const int istep, int& update_iter, std::ofstream& ofs, std::vector& etot_info, const Relax_Criteria& criteria); + bool start(UnitCell& ucell, + const ModuleBase::matrix& force, + const double& energy_in, + const int istep, + int& update_iter, + std::ofstream& ofs, + std::vector& etot_info, + const Relax_Criteria& criteria); //==================================================================== // Test seam; see the equivalent block in BFGS_Basic. Production code @@ -44,19 +52,23 @@ class Ions_Move_BFGS : public BFGS_Basic { bfgs_routine(lat0, istep, update_iter, ofs, etot_info, out_level, test_relax_method); } - void restart_bfgs_for_testing(const double& lat0, - int& update_iter, - std::ofstream& ofs, - const int test_relax_method) + void restart_bfgs_for_testing(const double& lat0, int& update_iter, std::ofstream& ofs, const int test_relax_method) { restart_bfgs(lat0, update_iter, ofs, test_relax_method); } private: bool init_done; - void bfgs_routine(const double& lat0, const int istep, int& update_iter, std::ofstream& ofs, std::vector& etot_info, const std::string& out_level, const int test_relax_method); + void bfgs_routine(const double& lat0, + const int istep, + int& update_iter, + std::ofstream& ofs, + std::vector& etot_info, + const std::string& out_level, + const int test_relax_method); void restart_bfgs(const double& lat0, int& update_iter, std::ofstream& ofs, const int test_relax_method); - bool first_step=true; // If it is the first step of the relaxation. The pos is only generated from ucell in the first step, and in the following steps, the pos is generated from the previous step. + bool first_step = true; // If it is the first step of the relaxation. The pos is only generated from ucell in the + // first step, and in the following steps, the pos is generated from the previous step. }; #endif diff --git a/source/source_relax/ions_move_bfgs2.cpp b/source/source_relax/ions_move_bfgs2.cpp index 5510c586aef..25d731bcea7 100644 --- a/source/source_relax/ions_move_bfgs2.cpp +++ b/source/source_relax/ions_move_bfgs2.cpp @@ -1,34 +1,35 @@ #include "ions_move_bfgs2.h" -#include "source_base/module_external/lapack_connector.h" -#include "source_io/module_parameter/parameter.h" + #include "ions_move_basic.h" #include "matrix_methods.h" +#include "source_base/module_external/lapack_connector.h" +#include "source_cell/print_cell.h" // lanshuyue add 2025-06-19 #include "source_cell/update_cell.h" -#include "source_cell/print_cell.h" // lanshuyue add 2025-06-19 +#include "source_io/module_parameter/parameter.h" //! initialize H0、H、pos0、force0、force -void Ions_Move_BFGS2::allocate(const int _size) +void Ions_Move_BFGS2::allocate(const int _size) { assert(_size > 0); - alpha=70;//default value in ase is 70 - maxstep=PARAM.inp.relax_bfgs_rmax; - size=_size; - sign=true; - H = std::vector>(3*size, std::vector(3*size, 0.0)); - - for (int i = 0; i < 3*size; ++i) + alpha = 70; // default value in ase is 70 + maxstep = PARAM.inp.relax_bfgs_rmax; + size = _size; + sign = true; + H = std::vector>(3 * size, std::vector(3 * size, 0.0)); + + for (int i = 0; i < 3 * size; ++i) { - H[i][i] = alpha; + H[i][i] = alpha; } - - pos = std::vector> (size, ModuleBase::Vector3(0.0, 0.0, 0.0)); - pos0 = std::vector(3*size, 0.0); - pos_taud = std::vector> (size, ModuleBase::Vector3(0.0, 0.0, 0.0)); - pos_taud0 = std::vector(3*size, 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, ModuleBase::Vector3(0.0, 0.0, 0.0)); + pos_taud0 = std::vector(3 * size, 0.0); dpos = std::vector>(size, ModuleBase::Vector3(0.0, 0.0, 0.0)); - force0 = std::vector(3*size, 0.0); + force0 = std::vector(3 * size, 0.0); force = std::vector>(size, ModuleBase::Vector3(0.0, 0.0, 0.0)); - steplength = std::vector(size, 0.0); + steplength = std::vector(size, 0.0); is_initialized = true; } @@ -37,170 +38,168 @@ void Ions_Move_BFGS2::reset() is_initialized = false; } - -bool Ions_Move_BFGS2::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell, std::ofstream& ofs_running) +bool Ions_Move_BFGS2::relax_step(const ModuleBase::matrix& _force, UnitCell& ucell, std::ofstream& ofs_running) { - if(!is_initialized) + if (!is_initialized) { allocate(ucell.nat); } - GetPos(ucell,pos); - GetPostaud(ucell,pos_taud); + GetPos(ucell, pos); + GetPostaud(ucell, pos_taud); ucell.ionic_position_updated = true; assert(_force.nr == force.size() && _force.nc == 3); - for(int i = 0; i < _force.nr; i++) + for (int i = 0; i < _force.nr; i++) { - 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; + 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;iPrepareStep(force,pos,H,pos0,force0,steplength,dpos,size,ucell); - this->DetermineStep(steplength,dpos,maxstep); + + this->PrepareStep(force, pos, H, pos0, force0, steplength, dpos, size, ucell); + this->DetermineStep(steplength, dpos, maxstep); this->UpdatePos(ucell); - this->CalculateLargestGrad(_force,ucell); + this->CalculateLargestGrad(_force, ucell); bool converged = this->IsRestrain(); - // print out geometry information during bfgs_trad relax - unitcell::print_tau(ucell.atoms,ucell.Coordinate,ucell.ntype,ucell.lat0,ofs_running); + // print out geometry information during bfgs_trad relax + unitcell::print_tau(ucell.atoms, ucell.Coordinate, ucell.ntype, ucell.lat0, ofs_running); return converged; } -void Ions_Move_BFGS2::GetPos(UnitCell& ucell,std::vector>& pos) +void Ions_Move_BFGS2::GetPos(UnitCell& ucell, std::vector>& pos) { assert(pos.size() == ucell.nat); - int iat=0; + int iat = 0; const double scale = ModuleBase::BOHR_TO_A * ucell.lat0; - for(int i=0;i>& pos_taud) +void Ions_Move_BFGS2::GetPostaud(UnitCell& ucell, std::vector>& pos_taud) { assert(pos_taud.size() == ucell.nat); - int iat=0; - for(int i=0;i>& force, - std::vector>& pos, - std::vector>& H, - std::vector& pos0, - std::vector& force0, - std::vector& steplength, - std::vector>& dpos, - int& size, - UnitCell& ucell) + std::vector>& pos, + std::vector>& H, + std::vector& pos0, + std::vector& force0, + std::vector& steplength, + std::vector>& dpos, + int& size, + UnitCell& ucell) { std::vector changedforce = ReshapeMToV(force); std::vector changedpos = ReshapeMToV(pos); - this->Update(changedpos, changedforce,H,ucell); - + this->Update(changedpos, changedforce, H, ucell); + //! call dysev - std::vector omega(3*size); - std::vector work(3*size*3*size); - int lwork=3*size*3*size; - int info=0; + std::vector omega(3 * size); + std::vector work(3 * size * 3 * size); + int lwork = 3 * size * 3 * size; + int info = 0; std::vector H_flat; - - for(const auto& row : H) + + for (const auto& row: H) { H_flat.insert(H_flat.end(), row.begin(), row.end()); } - - int value=3*size; - int* ptr=&value; - dsyev_("V","U",ptr,H_flat.data(),ptr,omega.data(),work.data(),&lwork,&info); - std::vector> V(3*size, std::vector(3*size, 0.0)); - for(int i = 0; i < 3*size; i++) + + int value = 3 * size; + int* ptr = &value; + dsyev_("V", "U", ptr, H_flat.data(), ptr, omega.data(), work.data(), &lwork, &info); + std::vector> V(3 * size, std::vector(3 * size, 0.0)); + for (int i = 0; i < 3 * size; i++) { - for(int j = 0; j < 3*size; j++) + for (int j = 0; j < 3 * size; j++) { - V[j][i] = H_flat[3*size*i + j]; + V[j][i] = H_flat[3 * size * i + j]; } } - std::vector a=DotInMAndV2(V, changedforce); - double threshold=1e-8; - for(int i = 0; i < a.size(); i++) + std::vector a = DotInMAndV2(V, changedforce); + double threshold = 1e-8; + for (int i = 0; i < a.size(); i++) { assert(std::abs(omega[i]) > threshold); - a[i]/=std::abs(omega[i]); + a[i] /= std::abs(omega[i]); } std::vector tmpdpos = DotInMAndV1(V, a); dpos = ReshapeVToM(tmpdpos); - for(int i = 0; i < size; i++) + for (int i = 0; i < size; i++) { double k = 0; - for(int j = 0; j < 3; j++) + for (int j = 0; j < 3; j++) { k += dpos[i][j] * dpos[i][j]; } steplength[i] = sqrt(k); } pos0 = ReshapeMToV(pos); - pos_taud0=ReshapeMToV(pos_taud); + pos_taud0 = ReshapeMToV(pos_taud); force0 = ReshapeMToV(force); } -void Ions_Move_BFGS2::Update(std::vector& pos, - std::vector& force, - std::vector>& H, - UnitCell& ucell) +void Ions_Move_BFGS2::Update(std::vector& pos, + std::vector& force, + std::vector>& H, + UnitCell& ucell) { - if(sign) + if (sign) { - sign=false; + sign = false; return; } - //std::vector dpos=this->VSubV(pos,pos0); - std::vector term=ReshapeMToV(pos_taud); + // std::vector dpos=this->VSubV(pos,pos0); + std::vector term = ReshapeMToV(pos_taud); std::vector dpos = VSubV(term, pos_taud0); - for(int i=0;i<3*size;i++) + for (int i = 0; i < 3 * size; i++) { double shortest_move = dpos[i]; - //dpos[i]/=ModuleBase::BOHR_TO_A; - //dpos[i]/=ucell.lat0; + // dpos[i]/=ModuleBase::BOHR_TO_A; + // dpos[i]/=ucell.lat0; for (int cell = -1; cell <= 1; ++cell) { const double now_move = dpos[i] + cell; @@ -209,91 +208,91 @@ void Ions_Move_BFGS2::Update(std::vector& pos, shortest_move = now_move; } } - //shortest_move=shortest_move*ModuleBase::BOHR_TO_A*ucell.lat0; - dpos[i]=shortest_move; + // shortest_move=shortest_move*ModuleBase::BOHR_TO_A*ucell.lat0; + dpos[i] = shortest_move; } - std::vector> c=ReshapeVToM(dpos); - for(int iat=0; iat> c = ReshapeVToM(dpos); + for (int iat = 0; iat < size; iat++) { - //Cartesian coordinate - //convert from Angstrom to unit of latvec (Bohr) + // Cartesian coordinate + // convert from Angstrom to unit of latvec (Bohr) - //convert unit + // convert unit ModuleBase::Vector3 move_ion_cart; 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; + // convert pos + ModuleBase::Vector3 move_ion_dr = move_ion_cart * ucell.latvec; int it = ucell.iat2it[iat]; int ia = ucell.iat2ia[iat]; Atom* atom = &ucell.atoms[it]; - if(atom->mbl[ia].x == 1) + if (atom->mbl[ia].x == 1) { dpos[iat * 3] = move_ion_dr.x; } - if(atom->mbl[ia].y == 1) + if (atom->mbl[ia].y == 1) { - dpos[iat * 3 + 1] = move_ion_dr.y ; + dpos[iat * 3 + 1] = move_ion_dr.y; } - if(atom->mbl[ia].z == 1) + if (atom->mbl[ia].z == 1) { - dpos[iat * 3 + 2] = move_ion_dr.z ; + dpos[iat * 3 + 2] = move_ion_dr.z; } } - double threshold=1e-7; + double threshold = 1e-7; double max_abs_dpos = 0.0; - for (const double value : dpos) + for (const double value: dpos) { max_abs_dpos = std::max(max_abs_dpos, std::abs(value)); } if (max_abs_dpos < threshold) { return; - } + } std::vector dforce = VSubV(force, force0); double a = DotInVAndV(dpos, dforce); std::vector dg = DotInMAndV1(H, dpos); double b = DotInVAndV(dpos, dg); - std::vector> term1=OuterVAndV(dforce, dforce); - std::vector> term2=OuterVAndV(dg, dg); - std::vector> term3=MPlus(term1, a); - std::vector> term4=MPlus(term2, b); + std::vector> term1 = OuterVAndV(dforce, dforce); + std::vector> term2 = OuterVAndV(dg, dg); + std::vector> term3 = MPlus(term1, a); + std::vector> term4 = MPlus(term2, b); H = MSubM(H, term3); H = MSubM(H, term4); } void Ions_Move_BFGS2::DetermineStep(std::vector& steplength, - std::vector>& dpos, - double& maxstep) + std::vector>& dpos, + double& maxstep) { std::vector::iterator maxsteplength = max_element(steplength.begin(), steplength.end()); double a = *maxsteplength; - double threshold=1e-10; + double threshold = 1e-10; assert(a > threshold); - if(a >= maxstep) + if (a >= maxstep) { double scale = maxstep / a; - for(int i = 0; i < size; i++) + for (int i = 0; i < size; i++) { - dpos[i].x*=scale; - dpos[i].y*=scale; - dpos[i].z*=scale; + dpos[i].x *= scale; + dpos[i].y *= scale; + dpos[i].z *= scale; } } } void Ions_Move_BFGS2::UpdatePos(UnitCell& ucell) { - double a[3*size]; - for(int i=0;iMAddM(pos, dpos);*/ } -void Ions_Move_BFGS2::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell) +void Ions_Move_BFGS2::CalculateLargestGrad(const ModuleBase::matrix& _force, UnitCell& ucell) { - std::vector grad= std::vector(3*size, 0.0); + std::vector grad = std::vector(3 * size, 0.0); int iat = 0; for (int it = 0; it < ucell.ntype; it++) { - Atom *atom = &ucell.atoms[it]; + Atom* atom = &ucell.atoms[it]; for (int ia = 0; ia < ucell.atoms[it].na; ia++) { for (int ik = 0; ik < 3; ++ik) @@ -355,7 +354,7 @@ void Ions_Move_BFGS2::CalculateLargestGrad(const ModuleBase::matrix& _force,Unit } } Ions_Move_Basic::largest_grad = 0.0; - for (int i = 0; i < 3*size; i++) + for (int i = 0; i < 3 * size; i++) { if (Ions_Move_Basic::largest_grad < std::abs(grad[i])) { @@ -366,9 +365,8 @@ void Ions_Move_BFGS2::CalculateLargestGrad(const ModuleBase::matrix& _force,Unit Ions_Move_Basic::largest_grad /= ucell.lat0; if (PARAM.inp.out_level == "ie") { - std::cout << " LARGEST GRAD (eV/Angstrom) : " << Ions_Move_Basic::largest_grad - * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A - << std::endl; + std::cout << " LARGEST GRAD (eV/Angstrom) : " + << Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << std::endl; } } diff --git a/source/source_relax/ions_move_bfgs2.h b/source/source_relax/ions_move_bfgs2.h index e84a1f6ad07..69d1c308834 100644 --- a/source/source_relax/ions_move_bfgs2.h +++ b/source/source_relax/ions_move_bfgs2.h @@ -1,22 +1,24 @@ #ifndef IONS_MOVE_BFGS2_H #define IONS_MOVE_BFGS2_H -#include -#include -#include -#include -#include #include "source_base/matrix.h" #include "source_base/matrix3.h" #include "source_cell/unitcell.h" +#include +#include +#include +#include +#include class Ions_Move_BFGS2 { -public: - void allocate(const int _size);//initialize parameters + public: + void allocate(const int _size); // initialize parameters void reset(); - bool relax_step(const ModuleBase::matrix& _force,UnitCell& ucell, std::ofstream& ofs_running);//a full iteration step + bool relax_step(const ModuleBase::matrix& _force, + UnitCell& ucell, + std::ofstream& ofs_running); // a full iteration step /// @brief whether allocate() has run and the optimiser holds live state. /// reset() clears it; ions_move_methods_test checks that. @@ -25,36 +27,48 @@ class Ions_Move_BFGS2 return is_initialized; } -private: + private: // The unit test drives the private step machinery directly (seeding the // Hessian and the previous-step state, then calling one stage at a time). friend class BFGSTest; - bool sign;//check if this is the first iteration - double alpha;//initialize H,diagonal element is alpha - double maxstep;//every movement smaller than maxstep - int size;//number of atoms - bool is_initialized=false; + bool sign; // check if this is the first iteration + double alpha; // initialize H,diagonal element is alpha + double maxstep; // every movement smaller than maxstep + int size; // number of atoms + bool is_initialized = false; - std::vector steplength;//the length of atoms displacement - std::vector> H;//Hessian matrix - std::vector force0;//force in previous step + 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 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_taud0; // atom pos in previous step(relative coordinates) 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 - bool 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 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 UpdatePos(UnitCell& ucell);//update ucell with the new coordinates - + + 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 + bool 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 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 UpdatePos(UnitCell& ucell); // update ucell with the new coordinates }; #endif // IONS_MOVE_BFGS2_H diff --git a/source/source_relax/ions_move_cg.cpp b/source/source_relax/ions_move_cg.cpp index 043944dfd98..0aa3d6071b5 100644 --- a/source/source_relax/ions_move_cg.cpp +++ b/source/source_relax/ions_move_cg.cpp @@ -1,7 +1,9 @@ #include "ions_move_cg.h" + #include "ions_move_basic.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" + #include double Ions_Move_CG::RELAX_CG_THR = -1.0; // default is 0.5 @@ -37,7 +39,15 @@ void Ions_Move_CG::allocate(const int dim) this->fmax = 0.0; } -bool Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const double &etot_in, const int istep, int& update_iter, std::ofstream& ofs, std::vector& etot_info, std::string& relax_method, const Relax_Criteria& criteria) +bool Ions_Move_CG::start(UnitCell& ucell, + const ModuleBase::matrix& force, + const double& etot_in, + const int istep, + int& update_iter, + std::ofstream& ofs, + std::vector& etot_info, + std::string& relax_method, + const Relax_Criteria& criteria) { ModuleBase::TITLE("Ions_Move_CG", "start"); assert(Ions_Move_Basic::dim > 0); @@ -77,7 +87,15 @@ bool Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const bool converged = false; if (flag == 0) { - converged = Ions_Move_Basic::check_converged(ucell, grad.data(), update_iter, ofs, etot_info, criteria.force_thr, criteria.force_thr_ev, criteria.out_level, criteria.test_relax_method); + converged = Ions_Move_Basic::check_converged(ucell, + grad.data(), + update_iter, + ofs, + etot_info, + criteria.force_thr, + criteria.force_thr_ev, + criteria.out_level, + criteria.test_relax_method); } if (converged) { @@ -88,7 +106,13 @@ bool Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const if (this->sd) { e0 = etot_in; - CG_Base::setup_cg_grad(dim, grad.data(), grad0.data(), cg_grad.data(), cg_grad0.data(), this->ncggrad, flag); + CG_Base::setup_cg_grad(dim, + grad.data(), + grad0.data(), + cg_grad.data(), + cg_grad0.data(), + this->ncggrad, + flag); this->ncggrad++; CG_Base::normalize(dim, cg_gradn.data(), cg_grad.data()); @@ -108,8 +132,7 @@ bool Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const if (relax_method == "cg_bfgs") { - if (Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A - < RELAX_CG_THR) + if (Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A < RELAX_CG_THR) { relax_method = "bfgs_trad"; } diff --git a/source/source_relax/ions_move_cg.h b/source/source_relax/ions_move_cg.h index efcce60d81a..1d91070630c 100644 --- a/source/source_relax/ions_move_cg.h +++ b/source/source_relax/ions_move_cg.h @@ -1,12 +1,13 @@ #ifndef IONS_MOVE_CG_H #define IONS_MOVE_CG_H +#include "cg_base.h" #include "relax_criteria.h" -#include -#include #include "source_base/matrix.h" #include "source_cell/unitcell.h" -#include "cg_base.h" + +#include +#include #include class Ions_Move_CG : public CG_Base @@ -16,7 +17,15 @@ class Ions_Move_CG : public CG_Base ~Ions_Move_CG() = default; void allocate(const int dim); - bool start(UnitCell &ucell, const ModuleBase::matrix &force, const double &etot, const int istep, int& update_iter, std::ofstream& ofs, std::vector& etot_info, std::string& relax_method, const Relax_Criteria& criteria); + bool start(UnitCell& ucell, + const ModuleBase::matrix& force, + const double& etot, + const int istep, + int& update_iter, + std::ofstream& ofs, + std::vector& etot_info, + std::string& relax_method, + const Relax_Criteria& criteria); static double RELAX_CG_THR; diff --git a/source/source_relax/ions_move_lbfgs.cpp b/source/source_relax/ions_move_lbfgs.cpp index fed26a3c82f..5fbb01c96c7 100644 --- a/source/source_relax/ions_move_lbfgs.cpp +++ b/source/source_relax/ions_move_lbfgs.cpp @@ -1,170 +1,177 @@ #include "ions_move_lbfgs.h" -#include "matrix_methods.h" + #include "ions_move_basic.h" +#include "matrix_methods.h" +#include "source_cell/print_cell.h" // mohan add 2025-06-19 #include "source_cell/update_cell.h" -#include "source_cell/print_cell.h" // mohan add 2025-06-19 -void Ions_Move_LBFGS::allocate(const int _size, const double relax_bfgs_rmax, const std::string& out_level) // initialize H0、H、pos0、force0、force +void Ions_Move_LBFGS::allocate(const int _size, + const double relax_bfgs_rmax, + const std::string& out_level) // initialize H0、H、pos0、force0、force { - alpha=70;//default value in ase is 70 - maxstep=relax_bfgs_rmax; - this->out_level=out_level; - size=_size; - memory=100; - iteration=0; - H = std::vector>(3*size, std::vector(3*size, 0.0)); - H0=1/alpha; - pos = std::vector> (size, ModuleBase::Vector3(0.0, 0.0, 0.0)); - pos0 = std::vector(3*size, 0.0); - pos_taud = std::vector> (size, ModuleBase::Vector3(0.0, 0.0, 0.0)); - pos_taud0 = std::vector(3*size, 0.0); + alpha = 70; // default value in ase is 70 + maxstep = relax_bfgs_rmax; + this->out_level = out_level; + size = _size; + memory = 100; + iteration = 0; + H = std::vector>(3 * size, std::vector(3 * size, 0.0)); + H0 = 1 / alpha; + pos = std::vector>(size, ModuleBase::Vector3(0.0, 0.0, 0.0)); + pos0 = std::vector(3 * size, 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, ModuleBase::Vector3(0.0, 0.0, 0.0)); - force0 = std::vector(3*size, 0.0); + force0 = std::vector(3 * size, 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(); + // l_search.init_line_search(); } -bool Ions_Move_LBFGS::relax_step(const ModuleBase::matrix _force,UnitCell& ucell,const double &etot, std::ofstream& ofs_running, const Relax_Criteria& criteria) +bool Ions_Move_LBFGS::relax_step(const ModuleBase::matrix _force, + UnitCell& ucell, + const double& etot, + std::ofstream& ofs_running, + const Relax_Criteria& criteria) { - get_pos(ucell,pos); - get_pos_taud(ucell,pos_taud); - //solver=p_esolver; + get_pos(ucell, pos); + get_pos_taud(ucell, pos_taud); + // solver=p_esolver; ucell.ionic_position_updated = true; - for(int i = 0; i < _force.nr; i++) + for (int i = 0; i < _force.nr; i++) { - for(int j=0;j<_force.nc;j++) + for (int j = 0; j < _force.nc; j++) { - force[i][j]=_force(i,j)*ModuleBase::Ry_to_eV/ModuleBase::BOHR_TO_A; + force[i][j] = _force(i, j) * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A; } } - int k=0; - for(int i=0;iprepare_step(force,pos,H,pos0,force0,dpos,ucell,etot); - this->determine_step(steplength,dpos,maxstep); + this->prepare_step(force, pos, H, pos0, force0, dpos, ucell, etot); + this->determine_step(steplength, dpos, maxstep); this->update_pos(ucell); - this->calculate_largest_grad(_force,ucell); + this->calculate_largest_grad(_force, ucell); bool converged = this->is_restrain(criteria); // mohan add 2025-06-22 - unitcell::print_tau(ucell.atoms,ucell.Coordinate,ucell.ntype,ucell.lat0,ofs_running); + unitcell::print_tau(ucell.atoms, ucell.Coordinate, ucell.ntype, ucell.lat0, ofs_running); return converged; } -void Ions_Move_LBFGS::get_pos(UnitCell& ucell,std::vector>& pos) +void Ions_Move_LBFGS::get_pos(UnitCell& ucell, std::vector>& pos) { - int k=0; - for(int i=0;i>& pos_taud) +void Ions_Move_LBFGS::get_pos_taud(UnitCell& ucell, std::vector>& pos_taud) { - int k=0; - for(int i=0;i>& force, - std::vector>& pos, - std::vector>& H, - std::vector& pos0, - std::vector& force0, - std::vector>& dpos, - UnitCell& ucell, - const double &etot) + std::vector>& pos, + std::vector>& H, + std::vector& pos0, + std::vector& force0, + std::vector>& dpos, + UnitCell& ucell, + const double& etot) { std::vector changedforce = ReshapeMToV(force); std::vector changedpos = ReshapeMToV(pos); - this->update(pos_taud,pos_taud0,changedforce,force0,ucell,iteration,memory,s,y,rho); - std::vector q=DotInVAndFloat(changedforce,-1); - int loopmax=std::min(memory,iteration); + this->update(pos_taud, pos_taud0, changedforce, force0, ucell, iteration, memory, s, y, rho); + std::vector q = DotInVAndFloat(changedforce, -1); + int loopmax = std::min(memory, iteration); std::vector a(loopmax); - for(int i=loopmax-1;i>=0;i--) + for (int i = loopmax - 1; i >= 0; i--) { - a[i]=rho[i]*DotInVAndV(s[i],q); - std::vector temp=DotInVAndFloat(y[i],a[i]); - q=VSubV(q,temp); + a[i] = rho[i] * DotInVAndV(s[i], q); + std::vector temp = DotInVAndFloat(y[i], a[i]); + q = VSubV(q, temp); } - std::vector z=DotInVAndFloat(q,H0); - for(int i=0;i z = DotInVAndFloat(q, H0); + for (int i = 0; i < loopmax; i++) { - double b=rho[i]*DotInVAndV(y[i],z); - std::vector temp=DotInVAndFloat(s[i],a[i]-b); - z=VAddV(z,temp); + double b = rho[i] * DotInVAndV(y[i], z); + std::vector temp = DotInVAndFloat(s[i], a[i] - b); + z = VAddV(z, temp); } - std::vector temp0=DotInVAndFloat(z,-1); - dpos=ReshapeVToM(temp0); - std::vector temp1=DotInVAndFloat(changedforce,-1); - //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); - std::vector temp2=DotInVAndFloat(temp0,1); - dpos=ReshapeVToM(temp2); - for(int i = 0; i < size; i++) + std::vector temp0 = DotInVAndFloat(z, -1); + dpos = ReshapeVToM(temp0); + std::vector temp1 = DotInVAndFloat(changedforce, -1); + // 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); + std::vector temp2 = DotInVAndFloat(temp0, 1); + dpos = ReshapeVToM(temp2); + for (int i = 0; i < size; i++) { double k = 0; - for(int j = 0; j < 3; j++) + for (int j = 0; j < 3; j++) { k += dpos[i][j] * dpos[i][j]; } steplength[i] = sqrt(k); } - iteration+=1; + iteration += 1; pos0 = ReshapeMToV(pos); - pos_taud0=ReshapeMToV(pos_taud); + pos_taud0 = ReshapeMToV(pos_taud); force0 = changedforce; } -void Ions_Move_LBFGS::update(std::vector>& pos_taud, - std::vector& pos_taud0, - std::vector& force, - std::vector& force0, - UnitCell& ucell, - int iteration, - int memory, - std::vector>& s, - std::vector>& y, - std::vector& rho) +void Ions_Move_LBFGS::update(std::vector>& pos_taud, + std::vector& pos_taud0, + std::vector& force, + std::vector& force0, + UnitCell& ucell, + int iteration, + int memory, + std::vector>& s, + std::vector>& y, + std::vector& rho) { - if(iteration>0) + if (iteration > 0) { - std::vector term=ReshapeMToV(pos_taud); - std::vector dpos =VSubV(term, pos_taud0); - for(int i=0;i<3*size;i++) + std::vector term = ReshapeMToV(pos_taud); + std::vector dpos = VSubV(term, pos_taud0); + for (int i = 0; i < 3 * size; i++) { double shortest_move = dpos[i]; for (int cell = -1; cell <= 1; ++cell) @@ -175,64 +182,66 @@ void Ions_Move_LBFGS::update(std::vector>& pos_taud, shortest_move = now_move; } } - dpos[i]=shortest_move; + dpos[i] = shortest_move; } - std::vector> c=ReshapeVToM(dpos); - for(int iat=0; iat> c = ReshapeVToM(dpos); + for (int iat = 0; iat < size; iat++) { - //Cartesian coordinate - //convert from Angstrom to unit of latvec (Bohr) + // Cartesian coordinate + // convert from Angstrom to unit of latvec (Bohr) - //convert unit + // convert unit ModuleBase::Vector3 move_ion_cart; 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; + // convert pos + ModuleBase::Vector3 move_ion_dr = move_ion_cart * ucell.latvec; int it = ucell.iat2it[iat]; int ia = ucell.iat2ia[iat]; Atom* atom = &ucell.atoms[it]; - if(atom->mbl[ia].x == 1) + if (atom->mbl[ia].x == 1) { dpos[iat * 3] = move_ion_dr.x; } - if(atom->mbl[ia].y == 1) + if (atom->mbl[ia].y == 1) { - dpos[iat * 3 + 1] = move_ion_dr.y ; + dpos[iat * 3 + 1] = move_ion_dr.y; } - if(atom->mbl[ia].z == 1) + if (atom->mbl[ia].z == 1) { - dpos[iat * 3 + 2] = move_ion_dr.z ; + dpos[iat * 3 + 2] = move_ion_dr.z; } } - std::vector dforce =VSubV(force0, force); - double rho0=1.0/DotInVAndV(dpos,dforce); + std::vector dforce = VSubV(force0, force); + double rho0 = 1.0 / DotInVAndV(dpos, dforce); s.push_back(dpos); y.push_back(dforce); rho.push_back(rho0); } - if(iteration>memory) + if (iteration > memory) { s.erase(s.begin()); y.erase(y.begin()); rho.erase(rho.begin()); } } -void Ions_Move_LBFGS::determine_step(std::vector& steplength,std::vector>& dpos,double& maxstep) +void Ions_Move_LBFGS::determine_step(std::vector& steplength, + std::vector>& dpos, + double& maxstep) { std::vector::iterator maxsteplength = max_element(steplength.begin(), steplength.end()); double a = *maxsteplength; - if(a >= maxstep) + if (a >= maxstep) { double scale = maxstep / a; - for(int i = 0; i < size; i++) + for (int i = 0; i < size; i++) { - for(int j=0;j<3;j++) + for (int j = 0; j < 3; j++) { - dpos[i][j]*=scale; + dpos[i][j] *= scale; } } } @@ -240,15 +249,15 @@ void Ions_Move_LBFGS::determine_step(std::vector& steplength,std::vector void Ions_Move_LBFGS::update_pos(UnitCell& ucell) { std::vector a(3 * size, 0.0); - for(int i=0;i grad= std::vector(3*size, 0.0); + std::vector grad = std::vector(3 * size, 0.0); int iat = 0; for (int it = 0; it < ucell.ntype; it++) { - Atom *atom = &ucell.atoms[it]; + Atom* atom = &ucell.atoms[it]; for (int ia = 0; ia < ucell.atoms[it].na; ia++) { for (int ik = 0; ik < 3; ++ik) @@ -276,7 +285,7 @@ void Ions_Move_LBFGS::calculate_largest_grad(const ModuleBase::matrix& _force,Un } } Ions_Move_Basic::largest_grad = 0.0; - for (int i = 0; i < 3*size; i++) + for (int i = 0; i < 3 * size; i++) { if (Ions_Move_Basic::largest_grad < std::abs(grad[i])) { @@ -286,10 +295,7 @@ void Ions_Move_LBFGS::calculate_largest_grad(const ModuleBase::matrix& _force,Un Ions_Move_Basic::largest_grad /= ucell.lat0; if (out_level == "ie") { - std::cout << " LARGEST GRAD (eV/Angstrom) : " - << Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A - << std::endl; + std::cout << " LARGEST GRAD (eV/Angstrom) : " + << Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << std::endl; } - } - diff --git a/source/source_relax/ions_move_lbfgs.h b/source/source_relax/ions_move_lbfgs.h index 02d2447e5ce..edf7a585f6a 100644 --- a/source/source_relax/ions_move_lbfgs.h +++ b/source/source_relax/ions_move_lbfgs.h @@ -1,19 +1,19 @@ #ifndef IONS_MOVE_LBFGS_H #define IONS_MOVE_LBFGS_H -#include -#include #include #include #include #include -//#include "line_search.h" +#include +#include +// #include "line_search.h" +#include "relax_criteria.h" #include "source_base/matrix.h" #include "source_base/matrix3.h" #include "source_cell/unitcell.h" #include "source_esolver/esolver.h" #include "source_esolver/esolver_ks.h" -#include "relax_criteria.h" /** * @class Ions_Move_LBFGS @@ -21,7 +21,7 @@ */ class Ions_Move_LBFGS { -public: + public: /** * @brief Initialize L-BFGS parameters * @param _size Number of atoms in system @@ -39,35 +39,35 @@ class Ions_Move_LBFGS */ bool relax_step(const ModuleBase::matrix _force, UnitCell& ucell, - const double &etot, + const double& etot, std::ofstream& ofs_running, const Relax_Criteria& criteria); -private: - //LineSearch l_search; - double alpha; ///< Initial Hessian diagonal element - double maxstep; ///< Maximum allowed step length - int size; ///< Number of atoms in system - int memory; ///< Number of previous steps to store - double H0; ///< Initial inverse Hessian approximation - int iteration; ///< Current iteration count - double energy; ///< Current system energy - double alpha_k; ///< Step size parameter - std::string out_level; ///< Output verbosity level + private: + // LineSearch l_search; + double alpha; ///< Initial Hessian diagonal element + double maxstep; ///< Maximum allowed step length + int size; ///< Number of atoms in system + int memory; ///< Number of previous steps to store + double H0; ///< Initial inverse Hessian approximation + int iteration; ///< Current iteration count + double energy; ///< Current system energy + double alpha_k; ///< Step size parameter + std::string out_level; ///< Output verbosity level - ModuleESolver::ESolver* solver = nullptr; ///< Structure solver - std::vector steplength;//the length of atoms displacement - std::vector> H;//Hessian matrix - std::vector force0;//force in previous step + ModuleESolver::ESolver* solver = nullptr; ///< Structure solver + 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 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_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> s; ///< Position difference vectors + std::vector> y; ///< Force difference vectors + std::vector rho; ///< Scalar products for L-BFGS update /** * @brief Prepare optimization step parameters @@ -79,7 +79,7 @@ class Ions_Move_LBFGS std::vector& force0, std::vector>& dpos, UnitCell& ucell, - const double &etot); + const double& etot); /** * @brief Judge if the relaxation has converged @@ -93,24 +93,21 @@ class Ions_Move_LBFGS * @param _force Current force matrix * @param ucell Unit cell being optimized */ - void calculate_largest_grad(const ModuleBase::matrix& _force, - UnitCell& ucell); + void calculate_largest_grad(const ModuleBase::matrix& _force, UnitCell& ucell); /** * @brief Extract atomic positions from unit cell * @param ucell Unit cell to read * @param pos Output position vector */ - void get_pos(UnitCell& ucell, - std::vector>& pos); + void get_pos(UnitCell& ucell, std::vector>& pos); /** * @brief Get fractional positions from unit cell * @param ucell Unit cell to read * @param pos_taud Output fractional positions */ - void get_pos_taud(UnitCell& ucell, - std::vector>& pos_taud); + void get_pos_taud(UnitCell& ucell, std::vector>& pos_taud); /** * @brief Update L-BFGS history buffers @@ -125,10 +122,10 @@ class Ions_Move_LBFGS * @param y Force differences buffer * @param rho Scalar products buffer */ - void update(std::vector>& pos_taud, - std::vector& pos_taud0, + void update(std::vector>& pos_taud, + std::vector& pos_taud0, std::vector& force, - std::vector& force0, + std::vector& force0, UnitCell& ucell, int iteration, int memory, @@ -143,14 +140,14 @@ class Ions_Move_LBFGS * @param maxstep Maximum allowed step length */ void determine_step(std::vector& steplength, - std::vector>& dpos, - double& maxstep); + std::vector>& dpos, + double& maxstep); /** * @brief Update atomic positions in unit cell * @param ucell Unit cell to update */ - void update_pos(UnitCell& ucell); + void update_pos(UnitCell& ucell); }; #endif \ No newline at end of file diff --git a/source/source_relax/ions_move_methods.cpp b/source/source_relax/ions_move_methods.cpp index 0cc9de4c1b1..6929921380a 100644 --- a/source/source_relax/ions_move_methods.cpp +++ b/source/source_relax/ions_move_methods.cpp @@ -1,12 +1,11 @@ #include "ions_move_methods.h" -#include - #include "ions_move_basic.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" #include "source_io/module_parameter/parameter.h" +#include Ions_Move_Methods::Ions_Move_Methods() { @@ -15,7 +14,7 @@ Ions_Move_Methods::~Ions_Move_Methods() { } -void Ions_Move_Methods::allocate(const int &natom, const std::string& relax_method) +void Ions_Move_Methods::allocate(const int& natom, const std::string& relax_method) { if (natom <= 0) { @@ -40,13 +39,13 @@ void Ions_Move_Methods::allocate(const int &natom, const std::string& relax_meth this->cg.allocate(Ions_Move_Basic::dim); this->bfgs.allocate(); } - else if(relax_method == "bfgs_trad") + else if (relax_method == "bfgs_trad") { - this->bfgs_trad.allocate(natom); + this->bfgs_trad.allocate(natom); } - else if(relax_method == "lbfgs") + else if (relax_method == "lbfgs") { - this->lbfgs.allocate(natom, PARAM.inp.relax_bfgs_rmax, PARAM.inp.out_level); + this->lbfgs.allocate(natom, PARAM.inp.relax_bfgs_rmax, PARAM.inp.out_level); } else { @@ -56,14 +55,14 @@ void Ions_Move_Methods::allocate(const int &natom, const std::string& relax_meth } // void Ions_Move_Methods::cal_movement(const int &istep, const ModuleBase::matrix &f, const double &etot) -void Ions_Move_Methods::cal_movement(const int &istep, - const int &force_step, - const ModuleBase::matrix &f, - const double &etot, - UnitCell &ucell, +void Ions_Move_Methods::cal_movement(const int& istep, + const int& force_step, + const ModuleBase::matrix& f, + const double& etot, + UnitCell& ucell, std::ofstream& ofs, std::string& relax_method, - const Relax_Criteria& criteria) + const Relax_Criteria& criteria) { ModuleBase::TITLE("Ions_Move_Methods", "init"); if (relax_method == "bfgs") @@ -84,7 +83,7 @@ void Ions_Move_Methods::cal_movement(const int &istep, } else if (relax_method == "bfgs_trad") { - converged_ = bfgs_trad.relax_step(f, ucell, ofs); + converged_ = bfgs_trad.relax_step(f, ucell, ofs); } else if (relax_method == "lbfgs") { diff --git a/source/source_relax/ions_move_methods.h b/source/source_relax/ions_move_methods.h index 00d0106585d..d52d4116ed9 100644 --- a/source/source_relax/ions_move_methods.h +++ b/source/source_relax/ions_move_methods.h @@ -1,16 +1,17 @@ #ifndef IONS_MOVE_METHODS_H #define IONS_MOVE_METHODS_H -#include -#include -#include #include "ions_move_basic.h" -#include "relax_criteria.h" #include "ions_move_bfgs.h" -#include "ions_move_cg.h" -#include "ions_move_sd.h" #include "ions_move_bfgs2.h" +#include "ions_move_cg.h" #include "ions_move_lbfgs.h" +#include "ions_move_sd.h" +#include "relax_criteria.h" + +#include +#include +#include class Ions_Move_Methods { @@ -18,12 +19,12 @@ class Ions_Move_Methods Ions_Move_Methods(); ~Ions_Move_Methods(); - void allocate(const int &natom, const std::string& relax_method); - void cal_movement(const int &istep, - const int &force_step, - const ModuleBase::matrix &f, - const double &etot, - UnitCell &ucell, + void allocate(const int& natom, const std::string& relax_method); + void cal_movement(const int& istep, + const int& force_step, + const ModuleBase::matrix& f, + const double& etot, + UnitCell& ucell, std::ofstream& ofs, std::string& relax_method, const Relax_Criteria& criteria); diff --git a/source/source_relax/ions_move_sd.cpp b/source/source_relax/ions_move_sd.cpp index 60466185d33..32b61bb25e5 100644 --- a/source/source_relax/ions_move_sd.cpp +++ b/source/source_relax/ions_move_sd.cpp @@ -1,10 +1,11 @@ #include "ions_move_sd.h" -#include #include "ions_move_basic.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" +#include + using namespace Ions_Move_Basic; Ions_Move_SD::Ions_Move_SD() : energy_saved(1.0e10) @@ -19,7 +20,14 @@ void Ions_Move_SD::allocate() pos_saved.resize(dim, 0.0); } -bool Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const double& etot_in, const int istep, int& update_iter, std::ofstream& ofs, std::vector& etot_info, const Relax_Criteria& criteria) +bool Ions_Move_SD::start(UnitCell& ucell, + const ModuleBase::matrix& force, + const double& etot_in, + const int istep, + int& update_iter, + std::ofstream& ofs, + std::vector& etot_info, + const Relax_Criteria& criteria) { ModuleBase::TITLE("Ions_Move_SD", "start"); @@ -39,7 +47,7 @@ bool Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const printf("in cheak_converged"); printf("pos[0]: %f\n", pos[0]); energy_saved = etot_in; - for (int i = 0; i < dim; i++) + for (int i = 0; i < dim; i++) { pos_saved[i] = pos[i]; } @@ -55,7 +63,15 @@ bool Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const } } - bool converged = Ions_Move_Basic::check_converged(ucell, grad.data(), update_iter, ofs, etot_info, criteria.force_thr, criteria.force_thr_ev, criteria.out_level, criteria.test_relax_method); + bool converged = Ions_Move_Basic::check_converged(ucell, + grad.data(), + update_iter, + ofs, + etot_info, + criteria.force_thr, + criteria.force_thr_ev, + criteria.out_level, + criteria.test_relax_method); if (converged) { Ions_Move_Basic::terminate(converged, update_iter, ucell, istep, ofs); diff --git a/source/source_relax/ions_move_sd.h b/source/source_relax/ions_move_sd.h index 8f62b436eb5..44b7bfa6887 100644 --- a/source/source_relax/ions_move_sd.h +++ b/source/source_relax/ions_move_sd.h @@ -2,10 +2,11 @@ #define IONS_MOVE_SD_H #include "relax_criteria.h" -#include -#include #include "source_base/matrix.h" #include "source_cell/unitcell.h" + +#include +#include #include #include @@ -31,7 +32,14 @@ class Ions_Move_SD ~Ions_Move_SD() = default; void allocate(void); - bool start(UnitCell& ucell, const ModuleBase::matrix& force, const double& etot, const int istep, int& update_iter, std::ofstream& ofs, std::vector& etot_info, const Relax_Criteria& criteria); + bool start(UnitCell& ucell, + const ModuleBase::matrix& force, + const double& etot, + const int istep, + int& update_iter, + std::ofstream& ofs, + std::vector& etot_info, + const Relax_Criteria& criteria); /// @brief Energy of the last accepted step. double get_energy_saved() const diff --git a/source/source_relax/lat_change_method.cpp b/source/source_relax/lat_change_method.cpp index 22ba35e257c..29900684e3c 100644 --- a/source/source_relax/lat_change_method.cpp +++ b/source/source_relax/lat_change_method.cpp @@ -17,11 +17,11 @@ void Lattice_Change_Methods::allocate() return; } -void Lattice_Change_Methods::cal_lattice_change(const int &istep, - const int &stress_step, - const ModuleBase::matrix &stress, - const double &etot, - UnitCell &ucell, +void Lattice_Change_Methods::cal_lattice_change(const int& istep, + const int& stress_step, + const ModuleBase::matrix& stress, + const double& etot, + UnitCell& ucell, std::ofstream& ofs, const Relax_Criteria& criteria) { diff --git a/source/source_relax/lat_change_method.h b/source/source_relax/lat_change_method.h index 9235a63457d..304cfae7b85 100644 --- a/source/source_relax/lat_change_method.h +++ b/source/source_relax/lat_change_method.h @@ -1,11 +1,12 @@ #ifndef LAT_CHANGE_METHOD_H #define LAT_CHANGE_METHOD_H +#include "lattice_change_basic.h" +#include "lattice_change_cg.h" #include "relax_criteria.h" + #include #include -#include "lattice_change_basic.h" -#include "lattice_change_cg.h" class Lattice_Change_Methods { @@ -16,11 +17,11 @@ class Lattice_Change_Methods void allocate(void); - void cal_lattice_change(const int &istep, - const int &stress_step, - const ModuleBase::matrix &stress, - const double &etot, - UnitCell &ucell, + void cal_lattice_change(const int& istep, + const int& stress_step, + const ModuleBase::matrix& stress, + const double& etot, + UnitCell& ucell, std::ofstream& ofs, const Relax_Criteria& criteria); diff --git a/source/source_relax/lattice_change_basic.cpp b/source/source_relax/lattice_change_basic.cpp index 9838980e44c..dca2699cdd5 100644 --- a/source/source_relax/lattice_change_basic.cpp +++ b/source/source_relax/lattice_change_basic.cpp @@ -11,7 +11,7 @@ int Lattice_Change_Basic::stress_step = 0; double Lattice_Change_Basic::lattice_change_ini = 0.01; // default is 0.01 std::string Lattice_Change_Basic::fixed_axes = "None"; -void Lattice_Change_Basic::setup_gradient(const UnitCell &ucell, double *lat, double *grad, ModuleBase::matrix &stress) +void Lattice_Change_Basic::setup_gradient(const UnitCell& ucell, double* lat, double* grad, ModuleBase::matrix& stress) { ModuleBase::TITLE("Lattice_Change_Basic", "setup_gradient"); @@ -97,7 +97,7 @@ void Lattice_Change_Basic::setup_gradient(const UnitCell &ucell, double *lat, do return; } -void Lattice_Change_Basic::change_lattice(UnitCell &ucell, double *move, double *lat, const bool fixed_ibrav) +void Lattice_Change_Basic::change_lattice(UnitCell& ucell, double* move, double* lat, const bool fixed_ibrav) { ModuleBase::TITLE("Lattice_Change_Basic", "change_lattice"); @@ -107,22 +107,23 @@ void Lattice_Change_Basic::change_lattice(UnitCell &ucell, double *move, double if (ModuleSymmetry::Symmetry::symm_flag && ucell.symm.nrotk > 0) { ModuleBase::matrix move_mat_t(3, 3); - for (int i = 0;i < 3;++i) + for (int i = 0; i < 3; ++i) { - for (int j = 0;j < 3;++j) + for (int j = 0; j < 3; ++j) { - move_mat_t(j, i) = move[i * 3 + j] / ucell.lat0; //transpose + move_mat_t(j, i) = move[i * 3 + j] / ucell.lat0; // transpose } } - ModuleBase::matrix symm_move_mat_t = (move_mat_t * ucell.G.to_matrix());//symmetrize (latvec^{-1} * move_mat)^T + ModuleBase::matrix symm_move_mat_t = (move_mat_t * ucell.G.to_matrix()); // symmetrize (latvec^{-1} * + // move_mat)^T ucell.symm.symmetrize_mat3(symm_move_mat_t, ucell.lat); - move_mat_t = symm_move_mat_t * ucell.latvec.Transpose().to_matrix();//G^{-1}=latvec^T + move_mat_t = symm_move_mat_t * ucell.latvec.Transpose().to_matrix(); // G^{-1}=latvec^T - for (int i = 0;i < 3;++i) + for (int i = 0; i < 3; ++i) { - for (int j = 0;j < 3;++j) + for (int j = 0; j < 3; ++j) { - move[i * 3 + j] = move_mat_t(j, i) * ucell.lat0;//transpose back + move[i * 3 + j] = move_mat_t(j, i) * ucell.lat0; // transpose back } } } @@ -230,7 +231,11 @@ void Lattice_Change_Basic::change_lattice(UnitCell &ucell, double *move, double return; } -bool Lattice_Change_Basic::check_converged(const UnitCell &ucell, ModuleBase::matrix &stress, double *grad, std::ofstream& ofs, const double& stress_thr) +bool Lattice_Change_Basic::check_converged(const UnitCell& ucell, + ModuleBase::matrix& stress, + double* grad, + std::ofstream& ofs, + const double& stress_thr) { ModuleBase::TITLE("Lattice_Change_Basic", "check_converged"); @@ -241,7 +246,7 @@ bool Lattice_Change_Basic::check_converged(const UnitCell &ucell, ModuleBase::ma { for (int i = 0; i < 3; i++) { - if (stress_ii_max < std::abs(stress(i, i))) + if (stress_ii_max < std::abs(stress(i, i))) { stress_ii_max = std::abs(stress(i, i)); } @@ -332,7 +337,7 @@ void Lattice_Change_Basic::terminate(const bool converged, std::ofstream& ofs) return; } -void Lattice_Change_Basic::setup_etot(const double &energy_in, std::vector& etot_info) +void Lattice_Change_Basic::setup_etot(const double& energy_in, std::vector& etot_info) { // etot_info[0] = etot (current total energy) // etot_info[1] = etot_p (previous total energy) diff --git a/source/source_relax/lattice_change_basic.h b/source/source_relax/lattice_change_basic.h index 651f328e88a..e00827fba79 100644 --- a/source/source_relax/lattice_change_basic.h +++ b/source/source_relax/lattice_change_basic.h @@ -1,12 +1,13 @@ #ifndef LATTICE_CHANGE_BASIC_H #define LATTICE_CHANGE_BASIC_H -#include -#include #include "relax_data.h" #include "source_base/matrix.h" #include "source_cell/unitcell.h" +#include +#include + /** * @namespace Lattice_Change_Basic * @brief Basic utilities and shared state for lattice relaxation algorithms. @@ -19,12 +20,12 @@ namespace Lattice_Change_Basic { // Shared state variables (referenced from Relax_Data for unified data sharing) -static int& dim = Relax_Data::dim_lattice; ///< Dimension of free variables (9 for full lattice) +static int& dim = Relax_Data::dim_lattice; ///< Dimension of free variables (9 for full lattice) static double& largest_grad = Relax_Data::largest_grad; ///< Largest gradient component // Lattice-specific parameters (not shared with ion movement) -extern int update_iter; ///< Number of successfully updated iterations -extern int stress_step; ///< Index of stress optimization step +extern int update_iter; ///< Number of successfully updated iterations +extern int stress_step; ///< Index of stress optimization step extern double lattice_change_ini; ///< Initial trust radius for lattice change (default: 0.01) extern std::string fixed_axes; ///< Fixed axes constraint ("None", "shape", "volume", or specific axes) @@ -35,7 +36,7 @@ extern std::string fixed_axes; ///< Fixed axes constraint ("None", "shape", " * @param grad Output gradient array (9 elements) * @param stress Stress tensor (3x3 matrix) */ -void setup_gradient(const UnitCell &ucell, double *lat, double *grad, ModuleBase::matrix &stress); +void setup_gradient(const UnitCell& ucell, double* lat, double* grad, ModuleBase::matrix& stress); /** * @brief Update lattice vectors according to displacement. @@ -43,7 +44,7 @@ void setup_gradient(const UnitCell &ucell, double *lat, double *grad, ModuleBase * @param move Displacement vector for lattice change (9 elements) * @param lat Current lattice vectors (9 elements) */ -void change_lattice(UnitCell &ucell, double *move, double *lat, const bool fixed_ibrav); +void change_lattice(UnitCell& ucell, double* move, double* lat, const bool fixed_ibrav); /** * @brief Check convergence based on stress threshold. @@ -53,7 +54,11 @@ void change_lattice(UnitCell &ucell, double *move, double *lat, const bool fixed * @param ofs Output stream for logging * @return true if converged, false otherwise */ -bool check_converged(const UnitCell &ucell, ModuleBase::matrix &stress, double *grad, std::ofstream& ofs, const double& stress_thr); +bool check_converged(const UnitCell& ucell, + ModuleBase::matrix& stress, + double* grad, + std::ofstream& ofs, + const double& stress_thr); /** * @brief Terminate lattice optimization and output results. @@ -68,6 +73,6 @@ void terminate(const bool converged, std::ofstream& ofs); * @param judgement Flag for SD method (true) or BFGS (false) * @param etot_info Vector containing [etot, etot_p] */ -void setup_etot(const double &energy_in, std::vector& etot_info); +void setup_etot(const double& energy_in, std::vector& etot_info); } // namespace Lattice_Change_Basic #endif diff --git a/source/source_relax/lattice_change_cg.cpp b/source/source_relax/lattice_change_cg.cpp index bc9a52e1407..9ffb91e9602 100644 --- a/source/source_relax/lattice_change_cg.cpp +++ b/source/source_relax/lattice_change_cg.cpp @@ -3,6 +3,7 @@ #include "lattice_change_basic.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" + #include using namespace Lattice_Change_Basic; @@ -38,7 +39,12 @@ void Lattice_Change_CG::allocate(void) this->fmax = 0.0; } -bool Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_in, const double &etot_in, std::ofstream& ofs, std::vector& etot_info, const Relax_Criteria& criteria) +bool Lattice_Change_CG::start(UnitCell& ucell, + const ModuleBase::matrix& stress_in, + const double& etot_in, + std::ofstream& ofs, + std::vector& etot_info, + const Relax_Criteria& criteria) { ModuleBase::TITLE("Lattice_Change_CG", "start"); @@ -94,7 +100,13 @@ bool Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_ if (this->sd) { e0 = etot_in; - CG_Base::setup_cg_grad(dim, grad.data(), grad0.data(), cg_grad.data(), cg_grad0.data(), this->ncggrad, flag); + CG_Base::setup_cg_grad(dim, + grad.data(), + grad0.data(), + cg_grad.data(), + cg_grad0.data(), + this->ncggrad, + flag); this->ncggrad++; CG_Base::normalize(dim, cg_gradn.data(), cg_grad.data()); diff --git a/source/source_relax/lattice_change_cg.h b/source/source_relax/lattice_change_cg.h index a786808fee0..baecf3f161c 100644 --- a/source/source_relax/lattice_change_cg.h +++ b/source/source_relax/lattice_change_cg.h @@ -1,11 +1,12 @@ #ifndef LATTICE_CHANGE_CG_H #define LATTICE_CHANGE_CG_H +#include "cg_base.h" #include "relax_criteria.h" -#include #include "source_base/matrix.h" #include "source_cell/unitcell.h" -#include "cg_base.h" + +#include #include class Lattice_Change_CG : public CG_Base @@ -16,7 +17,12 @@ class Lattice_Change_CG : public CG_Base ~Lattice_Change_CG() = default; void allocate(void); - bool start(UnitCell &ucell, const ModuleBase::matrix &stress_in, const double &etot, std::ofstream& ofs, std::vector& etot_info, const Relax_Criteria& criteria); + bool start(UnitCell& ucell, + const ModuleBase::matrix& stress_in, + const double& etot, + std::ofstream& ofs, + std::vector& etot_info, + const Relax_Criteria& criteria); private: // Unit tests need to inspect the saved previous-step state and to seed diff --git a/source/source_relax/line_search.cpp b/source/source_relax/line_search.cpp index b557e9d60ba..a66fdfb3bff 100644 --- a/source/source_relax/line_search.cpp +++ b/source/source_relax/line_search.cpp @@ -2,6 +2,7 @@ #include "source_base/global_function.h" #include "source_base/global_variable.h" + #include #include #include @@ -107,7 +108,7 @@ bool Line_Search::third_order(const double x, const double y, const double f, do { dmoveh = 4.0; } - + if (dmove > 2.0 * dmoveh || dmoveh > 2.0 * dmove || (fa * fb > 0 && dmove < 1.0) || (fa * fb < 0 && dmove > 1.0)) { @@ -195,7 +196,12 @@ void Line_Search::update_brent(const double x, const double y, const double f) } } -bool Line_Search::brent(const double x, const double y, const double f, double& xnew, const double conv_thr, std::ofstream& ofs_running) +bool Line_Search::brent(const double x, + const double y, + const double f, + double& xnew, + const double conv_thr, + std::ofstream& ofs_running) { ls_step++; diff --git a/source/source_relax/line_search.h b/source/source_relax/line_search.h index 35794d379af..0de66188785 100644 --- a/source/source_relax/line_search.h +++ b/source/source_relax/line_search.h @@ -9,12 +9,12 @@ class Line_Search Line_Search() {}; ~Line_Search() {}; - bool line_search(const bool restart, // whether to restart line search (when cg direction has changed) - const double x, // position for function at x - const double y, // value of function at x - const double f, // gradient of function at x - double& xnew, // postion where function has to be evaluated - const double conv_thr, // predicted change of function value of function at xnew + bool line_search(const bool restart, // whether to restart line search (when cg direction has changed) + const double x, // position for function at x + const double y, // value of function at x + const double f, // gradient of function at x + double& xnew, // postion where function has to be evaluated + const double conv_thr, // predicted change of function value of function at xnew std::ofstream& ofs_running); bool first_order(const double x, const double y, const double f, double& xnew); @@ -25,7 +25,12 @@ class Line_Search void update_brent(const double x, const double y, const double f); - bool brent(const double x, const double y, const double f, double& xnew, const double conv_thr, std::ofstream& ofs_running); + bool brent(const double x, + const double y, + const double f, + double& xnew, + const double conv_thr, + std::ofstream& ofs_running); private: int ls_step = 0; // step of line search diff --git a/source/source_relax/matrix_methods.cpp b/source/source_relax/matrix_methods.cpp index 96079d55e1d..a2b0df1e4d6 100644 --- a/source/source_relax/matrix_methods.cpp +++ b/source/source_relax/matrix_methods.cpp @@ -1,15 +1,12 @@ #include "matrix_methods.h" - - - -std::vector ReshapeMToV(std::vector>& matrix) +std::vector ReshapeMToV(std::vector>& matrix) { assert(!matrix.empty()); int size = matrix.size(); std::vector result; - result.reserve(3*size); - for (const auto& v : matrix) + result.reserve(3 * size); + for (const auto& v: matrix) { result.push_back(v.x); result.push_back(v.y); @@ -18,15 +15,15 @@ std::vector ReshapeMToV(std::vector>& matrix return result; } -std::vector> MAddM(std::vector>& a, - std::vector>& b) +std::vector> MAddM(std::vector>& a, std::vector>& b) { assert(!a.empty() && !b.empty()); assert(a.size() == b.size() && a[0].size() == b[0].size()); - std::vector> result = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); - for(int i = 0; i < a.size(); i++) + std::vector> result + = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); + for (int i = 0; i < a.size(); i++) { - for(int j = 0; j < a[0].size(); j++) + for (int j = 0; j < a[0].size(); j++) { result[i][j] = a[i][j] + b[i][j]; } @@ -34,52 +31,53 @@ std::vector> MAddM(std::vector>& a, return result; } -std::vector VSubV(std::vector& a, std::vector& b) +std::vector VSubV(std::vector& a, std::vector& b) { assert(a.size() == b.size()); std::vector result = std::vector(a.size(), 0.0); - for(int i = 0; i < a.size(); i++) + for (int i = 0; i < a.size(); i++) { result[i] = a[i] - b[i]; } 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, ModuleBase::Vector3(0.0, 0.0, 0.0)); - for(int i = 0; i < result.size(); i++) + std::vector> result + = std::vector>(matrix.size() / 3, ModuleBase::Vector3(0.0, 0.0, 0.0)); + for (int i = 0; i < result.size(); i++) { - result[i].x = matrix[i*3 ]; - result[i].y = matrix[i*3 + 1]; - result[i].z = matrix[i*3 + 2]; + result[i].x = matrix[i * 3]; + result[i].y = matrix[i * 3 + 1]; + result[i].z = matrix[i * 3 + 2]; } return result; } -std::vector DotInMAndV1(std::vector>& matrix, std::vector& vec) +std::vector DotInMAndV1(std::vector>& matrix, std::vector& vec) { assert(!matrix.empty()); assert(matrix[0].size() == vec.size()); std::vector result(matrix.size(), 0.0); - for(int i = 0; i < result.size(); i++) + for (int i = 0; i < result.size(); i++) { - for(int j = 0; j < vec.size(); j++) + for (int j = 0; j < vec.size(); j++) { result[i] += matrix[i][j] * vec[j]; } } return result; } -std::vector DotInMAndV2(std::vector>& matrix, std::vector& vec) +std::vector DotInMAndV2(std::vector>& matrix, std::vector& vec) { assert(!matrix.empty()); assert(matrix.size() == vec.size()); std::vector result(matrix.size(), 0.0); - for(int i = 0; i < result.size(); i++) + for (int i = 0; i < result.size(); i++) { - for(int j = 0; j < vec.size(); j++) + for (int j = 0; j < vec.size(); j++) { result[i] += matrix[j][i] * vec[j]; } @@ -87,24 +85,25 @@ std::vector DotInMAndV2(std::vector>& matrix, std::v return result; } -double DotInVAndV(std::vector& vec1, std::vector& vec2) +double DotInVAndV(std::vector& vec1, std::vector& vec2) { assert(vec1.size() == vec2.size()); double result = 0.0; - for(int i = 0; i < vec1.size(); i++) + for (int i = 0; i < vec1.size(); i++) { result += vec1[i] * vec2[i]; } return result; } -std::vector> OuterVAndV(std::vector& a, std::vector& b) +std::vector> OuterVAndV(std::vector& a, std::vector& b) { assert(a.size() == b.size()); - std::vector> result = std::vector>(a.size(), std::vector(b.size(), 0.0)); - for(int i = 0; i < a.size(); i++) + std::vector> result + = std::vector>(a.size(), std::vector(b.size(), 0.0)); + for (int i = 0; i < a.size(); i++) { - for(int j = 0; j < b.size(); j++) + for (int j = 0; j < b.size(); j++) { result[i][j] = a[i] * b[j]; } @@ -116,10 +115,11 @@ std::vector> MPlus(std::vector>& a, doub { assert(!a.empty()); assert(b != 0); - std::vector> result = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); - for(int i = 0; i < a.size(); i++) + std::vector> result + = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); + for (int i = 0; i < a.size(); i++) { - for(int j = 0; j < a[0].size(); j++) + for (int j = 0; j < a[0].size(); j++) { result[i][j] = a[i][j] / b; } @@ -131,10 +131,11 @@ std::vector> MSubM(std::vector>& a, std: { assert(!a.empty() && !b.empty()); assert(a.size() == b.size() && a[0].size() == b[0].size()); - std::vector> result = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); - for(int i = 0; i < a.size(); i++) + std::vector> result + = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); + for (int i = 0; i < a.size(); i++) { - for(int j = 0; j < a[0].size(); j++) + for (int j = 0; j < a[0].size(); j++) { result[i][j] = a[i][j] - b[i][j]; } @@ -142,22 +143,22 @@ std::vector> MSubM(std::vector>& a, std: return result; } -std::vector DotInVAndFloat(std::vector& vec, double b) +std::vector DotInVAndFloat(std::vector& vec, double b) { assert(b != 0); std::vector result(vec.size(), 0.0); - for(int i = 0; i < vec.size(); i++) + for (int i = 0; i < vec.size(); i++) { result[i] = vec[i] * b; } return result; } -std::vector VAddV(std::vector& a, std::vector& b) +std::vector VAddV(std::vector& a, std::vector& b) { assert(a.size() == b.size()); std::vector result = std::vector(a.size(), 0.0); - for(int i = 0; i < a.size(); i++) + for (int i = 0; i < a.size(); i++) { result[i] = a[i] + b[i]; } diff --git a/source/source_relax/matrix_methods.h b/source/source_relax/matrix_methods.h index fb6bba42b92..8057aa3265a 100644 --- a/source/source_relax/matrix_methods.h +++ b/source/source_relax/matrix_methods.h @@ -1,10 +1,10 @@ #ifndef MATRIX_METHODS #define MATRIX_METHODS -#include -#include #include "source_base/vector3.h" +#include +#include std::vector ReshapeMToV(std::vector>& matrix); std::vector> MAddM(std::vector>& a, std::vector>& b); @@ -17,8 +17,6 @@ double DotInVAndV(std::vector& vec1, std::vector& vec2); std::vector> OuterVAndV(std::vector& a, std::vector& b); std::vector> MPlus(std::vector>& a, double b); std::vector> MSubM(std::vector>& a, std::vector>& b); -std::vector DotInVAndFloat(std::vector& vec, double b); - - +std::vector DotInVAndFloat(std::vector& vec, double b); #endif \ No newline at end of file diff --git a/source/source_relax/relax_data.cpp b/source/source_relax/relax_data.cpp index bddffb81e02..da571e11171 100644 --- a/source/source_relax/relax_data.cpp +++ b/source/source_relax/relax_data.cpp @@ -5,7 +5,8 @@ int Relax_Data::dim_lattice = 0; double Relax_Data::largest_grad = 0.0; -void Relax_Data::allocate(const int dim_in) { +void Relax_Data::allocate(const int dim_in) +{ pos.resize(dim_in, 0.0); grad.resize(dim_in, 0.0); move.resize(dim_in, 0.0); diff --git a/source/source_relax/relax_data.h b/source/source_relax/relax_data.h index ec8ed2b1f52..429b1b5f404 100644 --- a/source/source_relax/relax_data.h +++ b/source/source_relax/relax_data.h @@ -5,14 +5,15 @@ /** * @brief Unified data structure for geometry optimization algorithms. - * + * * This class provides a common data container for all relaxation methods * (BFGS, CG, SD, etc.) to share and reuse key variables like positions, * gradients, and move vectors. It serves as the single source of truth * for optimization state across different algorithms. */ -class Relax_Data { -public: +class Relax_Data +{ + public: /** * @brief Default constructor. */ @@ -30,19 +31,19 @@ class Relax_Data { void allocate(const int dim_in); // Static members - shared global state across all relaxation instances - static int dim; ///< Dimension of free variables (3 * number of atoms) for ion movement - static int dim_lattice; ///< Dimension of free variables (9) for lattice change + static int dim; ///< Dimension of free variables (3 * number of atoms) for ion movement + static int dim_lattice; ///< Dimension of free variables (9) for lattice change - static double largest_grad; ///< Largest gradient component (force) in current step + static double largest_grad; ///< Largest gradient component (force) in current step // Instance members - per-iteration data - std::vector pos; ///< Current atomic positions in Bohr - std::vector grad; ///< Current gradient (negative force) in Ry/Bohr - std::vector move; ///< Displacement vector for next step + std::vector pos; ///< Current atomic positions in Bohr + std::vector grad; ///< Current gradient (negative force) in Ry/Bohr + std::vector move; ///< Displacement vector for next step - std::vector pos_p; ///< Previous atomic positions - std::vector grad_p; ///< Previous gradient - std::vector move_p; ///< Previous displacement + std::vector pos_p; ///< Previous atomic positions + std::vector grad_p; ///< Previous gradient + std::vector move_p; ///< Previous displacement }; #endif \ No newline at end of file diff --git a/source/source_relax/relax_driver.cpp b/source/source_relax/relax_driver.cpp index 1b672d3f736..7a0b124ee0d 100644 --- a/source/source_relax/relax_driver.cpp +++ b/source/source_relax/relax_driver.cpp @@ -1,20 +1,20 @@ #include "relax_driver.h" + #include "relax_history.h" #include "relax_stru_io.h" #include "socket_driver.h" #include "source_base/formatter.h" #include "source_base/global_file.h" +#include "source_base/module_out/read_exit_file.h" #include "source_io/module_json/output_info.h" #include "source_io/module_output/output_log.h" #include "source_io/module_output/print_info.h" -#include "source_base/module_out/read_exit_file.h" #include "source_io/module_parameter/parameter.h" -void Relax_Driver::relax_driver( - ModuleESolver::ESolver* p_esolver, - UnitCell& ucell, - const Input_para& inp, - std::ofstream& ofs_running) +void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, + UnitCell& ucell, + const Input_para& inp, + std::ofstream& ofs_running) { ModuleBase::TITLE("Relax_Driver", "relax_driver"); ModuleBase::timer::start("Relax_Driver", "relax_driver"); @@ -104,13 +104,11 @@ void Relax_Driver::init_relax(const int nat, const Input_para& inp) void Relax_Driver::iter_info(const std::vector& steps, const Input_para& inp) { if (inp.out_level == "ie" - && (inp.calculation == "relax" - || inp.calculation == "cell-relax" - || inp.calculation == "scf" - || inp.calculation == "nscf") - && (inp.esolver_type != "lr")) + && (inp.calculation == "relax" || inp.calculation == "cell-relax" || inp.calculation == "scf" + || inp.calculation == "nscf") + && (inp.esolver_type != "lr")) { - ModuleIO::print_screen(steps[2], steps[1], steps[0]+1); + ModuleIO::print_screen(steps[2], steps[1], steps[0] + 1); } #ifdef __JSON @@ -124,12 +122,12 @@ void Relax_Driver::iter_info(const std::vector& steps, const Input_para& in } void Relax_Driver::esolve(const int istep, - ModuleESolver::ESolver* p_esolver, - UnitCell& ucell, - const Input_para& inp, - ModuleBase::matrix& force, - ModuleBase::matrix& stress, - double& etot) + ModuleESolver::ESolver* p_esolver, + UnitCell& ucell, + const Input_para& inp, + ModuleBase::matrix& force, + ModuleBase::matrix& stress, + double& etot) { p_esolver->runner(ucell, istep); @@ -147,13 +145,13 @@ void Relax_Driver::esolve(const int istep, } bool Relax_Driver::relax_step(std::vector& steps, - ModuleESolver::ESolver* p_esolver, - UnitCell& ucell, - const Input_para& inp, - const ModuleBase::matrix& force, - const ModuleBase::matrix& stress, - const double etot, - std::ofstream& ofs_running) + ModuleESolver::ESolver* p_esolver, + UnitCell& ucell, + const Input_para& inp, + const ModuleBase::matrix& force, + const ModuleBase::matrix& stress, + const double etot, + std::ofstream& ofs_running) { // Guard: For non-relaxation calculations (scf, nscf, etc.), return true immediately // to ensure the main loop exits after one iteration. This provides robustness @@ -175,8 +173,7 @@ bool Relax_Driver::relax_step(std::vector& steps, } else { - converged = this->rl_old.relax_step(steps[0]+1, etot, ucell, force, - stress, steps[1], steps[2], ofs_running); + converged = this->rl_old.relax_step(steps[0] + 1, etot, ucell, force, stress, steps[1], steps[2], ofs_running); } ModuleIO::output_after_relax(converged, p_esolver->conv_esolver, ofs_running); @@ -184,11 +181,16 @@ bool Relax_Driver::relax_step(std::vector& steps, return converged; } -void Relax_Driver::stru_out(const int istep, UnitCell& ucell, const Input_para& inp, const double etot, const ModuleBase::matrix& stress, const ModuleBase::matrix& force) +void Relax_Driver::stru_out(const int istep, + UnitCell& ucell, + const Input_para& inp, + const double etot, + const ModuleBase::matrix& stress, + const ModuleBase::matrix& force) { // out_stru is effective for scf/nscf/relax/cell-relax (md writes STRU_MD_* via md_restartfreq) - if (inp.calculation != "relax" && inp.calculation != "cell-relax" - && inp.calculation != "scf" && inp.calculation != "nscf") + if (inp.calculation != "relax" && inp.calculation != "cell-relax" && inp.calculation != "scf" + && inp.calculation != "nscf") { return; } @@ -215,34 +217,45 @@ void Relax_Driver::stru_out(const int istep, UnitCell& ucell, const Input_para& if (is_relax) { const std::string now_file = out_dir_ + (inp.out_stru == 1 ? "STRU_NOW" : "STRU_NOW.cif"); - relax_stru_io::write_stru(ucell, inp, now_file, header, force, need_orb, - deepks_setorb_, my_rank_, inp.cal_force); + relax_stru_io::write_stru(ucell, + inp, + now_file, + header, + force, + need_orb, + deepks_setorb_, + my_rank_, + inp.cal_force); } // Numbered files per out_freq_ion: only meaningful for relaxation calculations if (is_relax && freq_ok) { - const std::string step_file = out_dir_ + "STRU" + std::to_string(istep + 1) - + (inp.out_stru == 1 ? "" : ".cif"); - relax_stru_io::write_stru(ucell, inp, step_file, header, force, need_orb, - deepks_setorb_, my_rank_, inp.cal_force); + const std::string step_file = out_dir_ + "STRU" + std::to_string(istep + 1) + (inp.out_stru == 1 ? "" : ".cif"); + relax_stru_io::write_stru(ucell, + inp, + step_file, + header, + force, + need_orb, + deepks_setorb_, + my_rank_, + inp.cal_force); } } -void Relax_Driver::json_out(ModuleESolver::ESolver* p_esolver, UnitCell& ucell, const Input_para& inp, const ModuleBase::matrix& force, const ModuleBase::matrix& stress) +void Relax_Driver::json_out(ModuleESolver::ESolver* p_esolver, + UnitCell& ucell, + const Input_para& inp, + const ModuleBase::matrix& force, + const ModuleBase::matrix& stress) { #ifdef __JSON Json::add_output_energy(p_esolver->cal_energy() * ModuleBase::Ry_to_eV); double unit_transform = ModuleBase::RYDBERG_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; double fac = ModuleBase::Ry_to_eV / 0.529177; - Json::add_output_cell_coo_stress_force(ucell, - force, - fac, - stress, - unit_transform, - inp.cal_force, - inp.cal_stress); + Json::add_output_cell_coo_stress_force(ucell, force, fac, stress, unit_transform, inp.cal_force, inp.cal_stress); #endif } @@ -270,8 +283,15 @@ void Relax_Driver::final_out(const int istep, // Only write forces when they were actually computed and belong // to the geometry being written. const bool write_force = inp.cal_force && geometry_evaluated; - relax_stru_io::write_stru(ucell, inp, final_file, header, force, need_orb, - deepks_setorb_, my_rank_, write_force); + relax_stru_io::write_stru(ucell, + inp, + final_file, + header, + force, + need_orb, + deepks_setorb_, + my_rank_, + write_force); } // relax_nmax == 0 is a valid dry-run mode: no relaxation step was ever @@ -309,17 +329,15 @@ void Relax_Driver::print_not_converged_summary(const Input_para& inp, std::ofstr // Unified not-converged summary for both relaxation paths so ASE and // users can read the method, the actual number of steps taken, and // the per-step force / stress history from the running log. - const std::vector& force_hist = inp.uses_simultaneous_relaxation() - ? rl.get_max_force_history() - : rl_old.get_max_force_history(); - const std::vector& stress_hist = inp.uses_simultaneous_relaxation() - ? rl.get_max_stress_history() - : rl_old.get_max_stress_history(); + const std::vector& force_hist + = inp.uses_simultaneous_relaxation() ? rl.get_max_force_history() : rl_old.get_max_force_history(); + const std::vector& stress_hist + = inp.uses_simultaneous_relaxation() ? rl.get_max_stress_history() : rl_old.get_max_stress_history(); const int ionic_steps = static_cast(force_hist.size()); const std::string method = inp.relax_method.empty() ? "unknown" : inp.relax_method; ofs_running << " Relaxation method: " << method << std::endl; - ofs_running << " Relaxation stopped after " << ionic_steps << " ionic step(s) (relax_nmax = " - << inp.relax_nmax << " reached)." << std::endl; + ofs_running << " Relaxation stopped after " << ionic_steps << " ionic step(s) (relax_nmax = " << inp.relax_nmax + << " reached)." << std::endl; if (!force_hist.empty()) { ofs_running << " Largest force per step (eV/Angstrom):" << format_relax_history(force_hist); diff --git a/source/source_relax/relax_driver.h b/source/source_relax/relax_driver.h index be4b99ca478..7050eba29b9 100644 --- a/source/source_relax/relax_driver.h +++ b/source/source_relax/relax_driver.h @@ -1,15 +1,16 @@ #ifndef RELAX_DRIVER_H #define RELAX_DRIVER_H +#include "relax_nsync.h" +#include "relax_sync.h" +#include "source_base/matrix.h" #include "source_cell/unitcell.h" #include "source_esolver/esolver.h" -#include "relax_sync.h" -#include "relax_nsync.h" #include "source_io/module_parameter/input_parameter.h" -#include "source_base/matrix.h" -#include + #include #include +#include /** * @brief Driver class for geometry relaxation calculations. @@ -21,8 +22,8 @@ class Relax_Driver { public: - Relax_Driver(){}; - ~Relax_Driver(){}; + Relax_Driver() {}; + ~Relax_Driver() {}; /** * @brief Main driver function for relaxation calculations. @@ -37,9 +38,9 @@ class Relax_Driver * @param ofs_running Output stream for running log. */ void relax_driver(ModuleESolver::ESolver* p_esolver, - UnitCell& ucell, - const Input_para& inp, - std::ofstream& ofs_running); + UnitCell& ucell, + const Input_para& inp, + std::ofstream& ofs_running); private: /// New relaxation optimizer (Relax class) @@ -81,13 +82,19 @@ class Relax_Driver * @param stress Output matrix for calculated stress. * @param etot Output total energy. */ - void esolve(const int istep, ModuleESolver::ESolver* p_esolver, UnitCell& ucell, - const Input_para& inp, ModuleBase::matrix& force, ModuleBase::matrix& stress, double& etot); + void esolve(const int istep, + ModuleESolver::ESolver* p_esolver, + UnitCell& ucell, + const Input_para& inp, + ModuleBase::matrix& force, + ModuleBase::matrix& stress, + double& etot); /** * @brief Perform one relaxation step. * - * @param steps Vector containing step counters (modified in place): steps[0]=istep, steps[1]=force_step, steps[2]=stress_step. + * @param steps Vector containing step counters (modified in place): steps[0]=istep, steps[1]=force_step, + * steps[2]=stress_step. * @param p_esolver Pointer to the energy solver. * @param ucell Reference to the unit cell. * @param inp Input parameters for the calculation. @@ -97,9 +104,14 @@ class Relax_Driver * @param ofs_running Output stream for running log. * @return True if relaxation converged, false otherwise. */ - bool relax_step(std::vector& steps, ModuleESolver::ESolver* p_esolver, UnitCell& ucell, - const Input_para& inp, const ModuleBase::matrix& force, const ModuleBase::matrix& stress, - const double etot, std::ofstream& ofs_running); + bool relax_step(std::vector& steps, + ModuleESolver::ESolver* p_esolver, + UnitCell& ucell, + const Input_para& inp, + const ModuleBase::matrix& force, + const ModuleBase::matrix& stress, + const double etot, + std::ofstream& ofs_running); /** * @brief Output structure files before relaxation move. @@ -110,7 +122,12 @@ class Relax_Driver * @param etot Total energy in Ry corresponding to this structure. * @param stress Stress matrix (3x3) in Ry/Bohr^3 corresponding to this structure. */ - void stru_out(const int istep, UnitCell& ucell, const Input_para& inp, const double etot, const ModuleBase::matrix& stress, const ModuleBase::matrix& force); + void stru_out(const int istep, + UnitCell& ucell, + const Input_para& inp, + const double etot, + const ModuleBase::matrix& stress, + const ModuleBase::matrix& force); /** * @brief Output JSON format results. @@ -121,8 +138,11 @@ class Relax_Driver * @param force Matrix of calculated forces. * @param stress Matrix of calculated stress. */ - void json_out(ModuleESolver::ESolver* p_esolver, UnitCell& ucell, const Input_para& inp, - const ModuleBase::matrix& force, const ModuleBase::matrix& stress); + void json_out(ModuleESolver::ESolver* p_esolver, + UnitCell& ucell, + const Input_para& inp, + const ModuleBase::matrix& force, + const ModuleBase::matrix& stress); /** * @brief Output final results after relaxation. diff --git a/source/source_relax/relax_nsync.cpp b/source/source_relax/relax_nsync.cpp index b577df2e070..0bd4732b295 100644 --- a/source/source_relax/relax_nsync.cpp +++ b/source/source_relax/relax_nsync.cpp @@ -1,4 +1,5 @@ #include "relax_nsync.h" + #include "relax_history.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" @@ -27,9 +28,9 @@ void IonCellOptimizer::print_converged_summary(const int istep, * * Allocates memory and initializes the appropriate relaxation methods: * - For "relax" calculation: only initializes Ions_Move_Methods - * - For "cell-relax" calculation: initializes both Ions_Move_Methods and + * - For "cell-relax" calculation: initializes both Ions_Move_Methods and * Lattice_Change_Methods - * + * * @param natom Number of atoms in the system */ void IonCellOptimizer::init_relax(const int& natom, const Input_para& inp) @@ -51,17 +52,17 @@ void IonCellOptimizer::init_relax(const int& natom, const Input_para& inp) /** * @brief Perform one step of relaxation (atomic and/or cell). - * + * * Main relaxation loop that coordinates atomic and cell relaxation: * 1. Check for maximum iteration limit * 2. Determine calculation mode (relax vs cell-relax) * 3. Perform atomic relaxation if needed and atoms can move * 4. If in cell-relax mode and atomic relaxation converged, perform cell relaxation - * + * * Convergence behavior: * - Returns false if relaxation is still in progress * - Returns true if relaxation has converged or maximum iterations reached - * + * * @param istep Current total iteration step * @param energy Total energy of the system * @param ucell Unit cell containing atomic positions and lattice vectors @@ -72,13 +73,13 @@ void IonCellOptimizer::init_relax(const int& natom, const Input_para& inp) * @return true if relaxation is converged, false otherwise */ bool IonCellOptimizer::relax_step(const int& istep, - const double& energy, - UnitCell& ucell, - ModuleBase::matrix force, - ModuleBase::matrix stress, - int& force_step, - int& stress_step, - std::ofstream& ofs_running) + const double& energy, + UnitCell& ucell, + ModuleBase::matrix force, + ModuleBase::matrix stress, + int& force_step, + int& stress_step, + std::ofstream& ofs_running) { ModuleBase::TITLE("IonCellOptimizer", "relax_step"); @@ -101,8 +102,8 @@ bool IonCellOptimizer::relax_step(const int& istep, // the actual number of ionic steps taken is the history length. const int ionic_steps = static_cast(max_force_history_.size()); ofs_running << " Relaxation method: " << inp_->relax_method << std::endl; - ofs_running << " Relaxation stopped after " << ionic_steps << " ionic step(s) (relax_nmax = " - << inp_->relax_nmax << " reached)." << std::endl; + ofs_running << " Relaxation stopped after " << ionic_steps + << " ionic step(s) (relax_nmax = " << inp_->relax_nmax << " reached)." << std::endl; if (!max_force_history_.empty()) { ofs_running << " Largest force per step (eV/Angstrom):" << format_relax_history(max_force_history_); @@ -138,7 +139,7 @@ bool IonCellOptimizer::relax_step(const int& istep, if (need_atom_relax) { assert(inp_->cal_force == 1); - + // Calculate and apply atomic movement std::string relax_method = inp_->relax_method; @@ -157,8 +158,8 @@ bool IonCellOptimizer::relax_step(const int& istep, const double max_force_ev = IMM.get_largest_grad() * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A; max_force_history_.push_back(max_force_ev); - ofs_running << " Largest force is " << max_force_ev - << " eV/Angstrom while threshold is " << inp_->force_thr_ev << " eV/Angstrom" << std::endl; + ofs_running << " Largest force is " << max_force_ev << " eV/Angstrom while threshold is " << inp_->force_thr_ev + << " eV/Angstrom" << std::endl; if (!atom_converged) { ofs_running << "\n Relaxation is not converged yet!" << std::endl; @@ -193,7 +194,7 @@ bool IonCellOptimizer::relax_step(const int& istep, if (need_cell_relax) { assert(inp_->cal_stress == 1); - + // Calculate and apply lattice change Relax_Criteria criteria; criteria.force_thr = inp_->force_thr; @@ -208,8 +209,8 @@ bool IonCellOptimizer::relax_step(const int& istep, const double max_stress_kbar = LCM.get_largest_grad(); max_stress_history_.push_back(max_stress_kbar); - ofs_running << " Largest stress is " << max_stress_kbar - << " kbar while threshold is " << inp_->stress_thr << " kbar" << std::endl; + ofs_running << " Largest stress is " << max_stress_kbar << " kbar while threshold is " << inp_->stress_thr + << " kbar" << std::endl; if (converged) { print_converged_summary(istep, true, ofs_running); @@ -222,12 +223,12 @@ bool IonCellOptimizer::relax_step(const int& istep, stress_step++; IMM.reset_after_cell_change(inp_->relax_method, ofs_running); ucell.cell_parameter_updated = true; - + // Update cell-related parameters after volume change unitcell::setup_cell_after_vc(ucell, ofs_running, inp_->nspin); ModuleBase::GlobalFunc::DONE(ofs_running, "SETUP UNITCELL"); } - + return converged; } else if (is_cell_relax && !unitcell::if_cell_can_change(ucell.lat_axis_free)) diff --git a/source/source_relax/relax_nsync.h b/source/source_relax/relax_nsync.h index fbe3269a03f..ad06b04a01f 100644 --- a/source/source_relax/relax_nsync.h +++ b/source/source_relax/relax_nsync.h @@ -5,18 +5,19 @@ #include "lat_change_method.h" #include "source_cell/unitcell.h" #include "source_io/module_parameter/input_parameter.h" + #include #include /** * @class IonCellOptimizer * @brief Optimizer for ionic position relaxation and cell parameter optimization. - * + * * This class handles both atomic relaxation and cell relaxation calculations. * It manages the geometry optimization process by coordinating between * Ions_Move_Methods (for atomic position updates) and Lattice_Change_Methods * (for cell parameter updates). - * + * * The optimization follows a sequential approach: 1. First perform atomic relaxation until convergence 2. Then perform cell relaxation (only in cell-relax mode) @@ -32,11 +33,11 @@ class IonCellOptimizer /** * @brief Perform one step of relaxation (atomic and/or cell). - * + * * This is the main interface for relaxation. Depending on the calculation type * ("relax" or "cell-relax"), it will perform atomic relaxation, cell relaxation, * or both in sequence. - * + * * @param istep Current total iteration step * @param energy Total energy of the system * @param ucell Unit cell containing atomic positions and lattice vectors @@ -57,8 +58,8 @@ class IonCellOptimizer std::ofstream& ofs_running); private: - Ions_Move_Methods IMM; ///< Ionic movement methods for atom relaxation - Lattice_Change_Methods LCM; ///< Lattice change methods for cell relaxation + Ions_Move_Methods IMM; ///< Ionic movement methods for atom relaxation + Lattice_Change_Methods LCM; ///< Lattice change methods for cell relaxation const Input_para* inp_ = nullptr; /// Largest force of each ionic step (eV/Angstrom), kept across cell changes. @@ -77,14 +78,18 @@ class IonCellOptimizer * (true for cell-relax with a movable lattice). * @param ofs_running Output stream for the running log. */ - void print_converged_summary(const int istep, - const bool with_stress, - std::ofstream& ofs_running) const; + void print_converged_summary(const int istep, const bool with_stress, std::ofstream& ofs_running) const; public: /// Read-only observers of the per-step convergence history, for the final summary. - const std::vector& get_max_force_history() const { return max_force_history_; } - const std::vector& get_max_stress_history() const { return max_stress_history_; } + const std::vector& get_max_force_history() const + { + return max_force_history_; + } + const std::vector& get_max_stress_history() const + { + return max_stress_history_; + } }; #endif \ No newline at end of file diff --git a/source/source_relax/relax_stru_io.cpp b/source/source_relax/relax_stru_io.cpp index 8d3ec0efe0b..0f191ac644b 100644 --- a/source/source_relax/relax_stru_io.cpp +++ b/source/source_relax/relax_stru_io.cpp @@ -49,10 +49,9 @@ std::string build_stru_header(const int istep, // stress in kbar: Ry/Bohr^3 -> kbar, always exactly 3 lines for // downstream parsers. N/A marks uncomputed stress; no inline comment is // appended so that every stress line has the same token count. - const double stress_transform = ModuleBase::RYDBERG_SI - / (ModuleBase::BOHR_RADIUS_SI * ModuleBase::BOHR_RADIUS_SI - * ModuleBase::BOHR_RADIUS_SI) - * 1.0e-8; + const double stress_transform + = ModuleBase::RYDBERG_SI + / (ModuleBase::BOHR_RADIUS_SI * ModuleBase::BOHR_RADIUS_SI * ModuleBase::BOHR_RADIUS_SI) * 1.0e-8; const bool stress_valid = inp.cal_stress && geometry_evaluated; if (stress_valid) { @@ -81,7 +80,8 @@ std::string build_stru_header(const int istep, } else if (!inp.cal_stress && !inp.cal_force) { - note = "# NOTE: stress not computed (cal_stress=0); forces not computed (cal_force=0); per-atom f fields omitted intentionally"; + note = "# NOTE: stress not computed (cal_stress=0); forces not computed (cal_force=0); per-atom f fields " + "omitted intentionally"; } else if (!inp.cal_stress) { diff --git a/source/source_relax/relax_stru_io.h b/source/source_relax/relax_stru_io.h index 205da87c76a..2830f4b09b7 100644 --- a/source/source_relax/relax_stru_io.h +++ b/source/source_relax/relax_stru_io.h @@ -16,64 +16,64 @@ */ namespace relax_stru_io { - /** - * @brief Build the header comment for a structure file. - * - * Contains the ABACUS version, a timestamp, the relaxation step number, - * the total energy in eV, the 3x3 stress tensor in kbar and a trailing - * NOTE line. The header is always exactly 7 lines so that downstream - * parsers can skip a fixed-size comment block. - * - * @param istep Current (zero-based) relaxation step; printed as istep + 1. - * @param etot Total energy in Ry. - * @param stress Stress tensor (3x3) in Ry/Bohr^3. - * @param inp Input parameters (cal_force, cal_stress select the NOTE text). - * @param is_final If true, the step label is marked "(FINAL)". - * @param geometry_evaluated false when the geometry was proposed by the - * optimizer but never evaluated; the energy and stress are then - * marked N/A and the NOTE explains that forces are omitted. - * @return The formatted header string. - */ - std::string build_stru_header(const int istep, - const double etot, - const ModuleBase::matrix& stress, - const Input_para& inp, - const bool is_final, - const bool geometry_evaluated); +/** + * @brief Build the header comment for a structure file. + * + * Contains the ABACUS version, a timestamp, the relaxation step number, + * the total energy in eV, the 3x3 stress tensor in kbar and a trailing + * NOTE line. The header is always exactly 7 lines so that downstream + * parsers can skip a fixed-size comment block. + * + * @param istep Current (zero-based) relaxation step; printed as istep + 1. + * @param etot Total energy in Ry. + * @param stress Stress tensor (3x3) in Ry/Bohr^3. + * @param inp Input parameters (cal_force, cal_stress select the NOTE text). + * @param is_final If true, the step label is marked "(FINAL)". + * @param geometry_evaluated false when the geometry was proposed by the + * optimizer but never evaluated; the energy and stress are then + * marked N/A and the NOTE explains that forces are omitted. + * @return The formatted header string. + */ +std::string build_stru_header(const int istep, + const double etot, + const ModuleBase::matrix& stress, + const Input_para& inp, + const bool is_final, + const bool geometry_evaluated); - /** - * @brief Whether orbital output is required in the structure file. - * - * True for lcao / lcao_in_pw bases, and for pw with nao-initialized - * wavefunctions. - */ - bool need_orbital(const Input_para& inp); +/** + * @brief Whether orbital output is required in the structure file. + * + * True for lcao / lcao_in_pw bases, and for pw with nao-initialized + * wavefunctions. + */ +bool need_orbital(const Input_para& inp); - /** - * @brief Write one structure file in STRU (out_stru == 1) or CIF - * (out_stru == 2) format. - * - * @param ucell Unit cell to write. - * @param inp Input parameters (nspin, out_mul, calculation). - * @param filename Output file path. - * @param header Header comment built by build_stru_header. - * @param force Atomic forces (used for STRU output when has_force is true). - * @param need_orb Whether to write orbitals (from need_orbital). - * @param deepks_setorb DeePKS setorb flag forwarded to the STRU writer. - * @param my_rank MPI rank of the calling process. - * @param has_force true to emit per-atom f fields; forwarded to - * print_stru_file. Callers must pass false when the forces do not - * belong to the geometry being written (e.g. early relax exit). - */ - void write_stru(UnitCell& ucell, - const Input_para& inp, - const std::string& filename, - const std::string& header, - const ModuleBase::matrix& force, - const bool need_orb, - const bool deepks_setorb, - const int my_rank, - const bool has_force); +/** + * @brief Write one structure file in STRU (out_stru == 1) or CIF + * (out_stru == 2) format. + * + * @param ucell Unit cell to write. + * @param inp Input parameters (nspin, out_mul, calculation). + * @param filename Output file path. + * @param header Header comment built by build_stru_header. + * @param force Atomic forces (used for STRU output when has_force is true). + * @param need_orb Whether to write orbitals (from need_orbital). + * @param deepks_setorb DeePKS setorb flag forwarded to the STRU writer. + * @param my_rank MPI rank of the calling process. + * @param has_force true to emit per-atom f fields; forwarded to + * print_stru_file. Callers must pass false when the forces do not + * belong to the geometry being written (e.g. early relax exit). + */ +void write_stru(UnitCell& ucell, + const Input_para& inp, + const std::string& filename, + const std::string& header, + const ModuleBase::matrix& force, + const bool need_orb, + const bool deepks_setorb, + const int my_rank, + const bool has_force); } // namespace relax_stru_io #endif diff --git a/source/source_relax/relax_sync.cpp b/source/source_relax/relax_sync.cpp index 1b8081ea847..34df2b21cce 100644 --- a/source/source_relax/relax_sync.cpp +++ b/source/source_relax/relax_sync.cpp @@ -1,12 +1,11 @@ #include "relax_sync.h" #include "relax_history.h" - #include "source_base/matrix3.h" #include "source_base/parallel_common.h" #include "source_base/tool_title.h" -#include "source_cell/update_cell.h" #include "source_cell/print_cell.h" +#include "source_cell/update_cell.h" #include "source_io/module_parameter/input_parameter.h" #include "source_relax/ions_move_basic.h" @@ -103,7 +102,10 @@ bool Relax::relax_step(UnitCell& ucell, return relax_done; } -bool Relax::setup_gradient(const UnitCell& ucell, const ModuleBase::matrix& force, const ModuleBase::matrix& stress, std::ofstream& ofs_running) +bool Relax::setup_gradient(const UnitCell& ucell, + const ModuleBase::matrix& force, + const ModuleBase::matrix& stress, + std::ofstream& ofs_running) { ModuleBase::TITLE("Relax", "setup_gradient"); @@ -136,8 +138,10 @@ bool Relax::setup_gradient(const UnitCell& ucell, const ModuleBase::matrix& forc return force_converged; } -double Relax::setup_ion_gradient(const UnitCell& ucell, const ModuleBase::matrix& force, - bool& force_converged, std::ofstream& ofs_running) +double Relax::setup_ion_gradient(const UnitCell& ucell, + const ModuleBase::matrix& force, + bool& force_converged, + std::ofstream& ofs_running) { grad_ion.zero_out(); ModuleBase::matrix force_eva = force * (ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A); // convert to eV/Angstrom @@ -199,15 +203,16 @@ double Relax::setup_ion_gradient(const UnitCell& ucell, const ModuleBase::matrix etot_p = etot; } - max_force_history_.push_back(max_grad); - ofs_running << "\n Largest force is " << max_grad << - " eV/Angstrom while threshold is " << inp_->force_thr_ev << " eV/Angstrom" << std::endl; + ofs_running << "\n Largest force is " << max_grad << " eV/Angstrom while threshold is " << inp_->force_thr_ev + << " eV/Angstrom" << std::endl; return max_grad; } -void Relax::setup_cell_gradient(const UnitCell& ucell, const ModuleBase::matrix& stress, - bool& force_converged, std::ofstream& ofs_running) +void Relax::setup_cell_gradient(const UnitCell& ucell, + const ModuleBase::matrix& stress, + bool& force_converged, + std::ofstream& ofs_running) { grad_cell.zero_out(); ModuleBase::matrix stress_ev = stress * (ucell.omega * ModuleBase::Ry_to_eV); @@ -286,8 +291,8 @@ void Relax::setup_cell_gradient(const UnitCell& ucell, const ModuleBase::matrix& } max_stress_history_.push_back(largest_grad); - ofs_running << " Largest stress is " << largest_grad << " kbar while threshold is " - << inp_->stress_thr << " kbar" << std::endl; + ofs_running << " Largest stress is " << largest_grad << " kbar while threshold is " << inp_->stress_thr << " kbar" + << std::endl; } void Relax::print_gradient_summary(bool force_converged, std::ofstream& ofs_running) @@ -688,10 +693,10 @@ void Relax::update_ion_positions(UnitCell& ucell, double fac, std::ofstream& ofs ucell.symm.symmetrize_vec3_nat(move_ion.data()); } - unitcell::update_pos_taud(ucell.lat,move_ion.data(),ucell.ntype,ucell.nat,ucell.atoms); + unitcell::update_pos_taud(ucell.lat, move_ion.data(), ucell.ntype, ucell.nat, ucell.atoms); // Print the structure file. - unitcell::print_tau(ucell.atoms,ucell.Coordinate,ucell.ntype,ucell.lat0,ofs_running); + unitcell::print_tau(ucell.atoms, ucell.Coordinate, ucell.ntype, ucell.lat0, ofs_running); } void Relax::update_reciprocal_cell(UnitCell& ucell, std::ofstream& ofs_running) diff --git a/source/source_relax/relax_sync.h b/source/source_relax/relax_sync.h index ebb2e6dc3f2..8cdf70a27e8 100644 --- a/source/source_relax/relax_sync.h +++ b/source/source_relax/relax_sync.h @@ -8,6 +8,7 @@ #include "source_base/matrix3.h" #include "source_cell/unitcell.h" #include "source_io/module_parameter/input_parameter.h" + #include #include @@ -33,21 +34,28 @@ class Relax // constraints are considered here // also check if relaxation has converged // based on threshold in force & stress - bool setup_gradient(const UnitCell& ucell, const ModuleBase::matrix& force, const ModuleBase::matrix& stress, std::ofstream& ofs_running); + bool setup_gradient(const UnitCell& ucell, + const ModuleBase::matrix& force, + const ModuleBase::matrix& stress, + std::ofstream& ofs_running); // Compute the ionic gradient from the force (in eV/Angstrom), honoring // per-atom move flags, and return the largest component. Records the value // in the force history and the running log, prints the optional out_level // diagnostics, and updates force_converged. - double setup_ion_gradient(const UnitCell& ucell, const ModuleBase::matrix& force, - bool& force_converged, std::ofstream& ofs_running); + double setup_ion_gradient(const UnitCell& ucell, + const ModuleBase::matrix& force, + bool& force_converged, + std::ofstream& ofs_running); // Compute the cell gradient from the stress, applying the fixed_axes // constraints (shape / volume / per-axis), and update force_converged. // Records the largest stress in the history and the running log. Only runs // when if_cell_moves is true. - void setup_cell_gradient(const UnitCell& ucell, const ModuleBase::matrix& stress, - bool& force_converged, std::ofstream& ofs_running); + void setup_cell_gradient(const UnitCell& ucell, + const ModuleBase::matrix& stress, + bool& force_converged, + std::ofstream& ofs_running); // Print the converged / not-converged summary to the running log. On // convergence, reports the method, the converged step (istep + 1, since @@ -152,8 +160,14 @@ class Relax public: /// Read-only observers of the per-step convergence history, for the final summary. - const std::vector& get_max_force_history() const { return max_force_history_; } - const std::vector& get_max_stress_history() const { return max_stress_history_; } + const std::vector& get_max_force_history() const + { + return max_force_history_; + } + const std::vector& get_max_stress_history() const + { + return max_stress_history_; + } }; #endif \ No newline at end of file diff --git a/source/source_relax/socket_driver.cpp b/source/source_relax/socket_driver.cpp index d63c83441d5..03d19cb2764 100644 --- a/source/source_relax/socket_driver.cpp +++ b/source/source_relax/socket_driver.cpp @@ -1,23 +1,23 @@ #include "socket_driver.h" -#include "source_relax/socket_handlers.h" -#include "source_relax/socket_utils.h" #include "source_base/timer.h" #include "source_cell/unitcell.h" #include "source_esolver/esolver.h" #include "source_io/module_parameter/input_parameter.h" +#include "source_relax/socket_handlers.h" +#include "source_relax/socket_utils.h" #include #include #include +using SocketHandlers::DriverContext; +using SocketUtils::address; using SocketUtils::ComputedFrame; using SocketUtils::DriverState; using SocketUtils::ipi_cell_bohr_from_unitcell; using SocketUtils::is_root; using SocketUtils::quit_if_root_failed; -using SocketUtils::address; -using SocketHandlers::DriverContext; namespace { @@ -79,8 +79,7 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, while (true) { - const std::string header - = SocketHandlers::read_header_bcast(socket, context.state); + const std::string header = SocketHandlers::read_header_bcast(socket, context.state); if (header.empty()) { log_peer_closed(ofs_running); @@ -109,8 +108,7 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, } else { - quit_if_root_failed(is_root() ? 1 : 0, - is_root() ? "unknown i-PI header: " + header : ""); + quit_if_root_failed(is_root() ? 1 : 0, is_root() ? "unknown i-PI header: " + header : ""); } } } diff --git a/source/source_relax/socket_frame.cpp b/source/source_relax/socket_frame.cpp index 1a009879ddb..9dc745bb1a3 100644 --- a/source/source_relax/socket_frame.cpp +++ b/source/source_relax/socket_frame.cpp @@ -24,9 +24,7 @@ CellValidation make_failed_cell_validation() return result; } -bool validate_cell_entries(const Matrix9& cell_wire, - const Matrix9& inverse_wire, - std::string& message) +bool validate_cell_entries(const Matrix9& cell_wire, const Matrix9& inverse_wire, std::string& message) { if (!is_finite_matrix(cell_wire) || !is_finite_matrix(inverse_wire)) { @@ -51,8 +49,7 @@ bool validate_cell_tolerances(const double max_condition_number, return true; } -bool compute_cell_scale(const Matrix9& cell_wire, double& scale, Matrix9& scaled_cell, - std::string& message) +bool compute_cell_scale(const Matrix9& cell_wire, double& scale, Matrix9& scaled_cell, std::string& message) { scale = 0.0; for (std::size_t index = 0; index < cell_wire.size(); ++index) @@ -83,10 +80,8 @@ bool compute_cell_determinant(const double scale, return false; } const long double scale_long = scale; - const long double determinant - = determinant_scaled * scale_long * scale_long * scale_long; - if (!std::isfinite(determinant) - || determinant > static_cast(std::numeric_limits::max())) + const long double determinant = determinant_scaled * scale_long * scale_long * scale_long; + if (!std::isfinite(determinant) || determinant > static_cast(std::numeric_limits::max())) { message = "cell determinant is not representable as a finite double"; return false; @@ -137,8 +132,7 @@ bool compute_condition_number(const double singular_values[kMatrixDimension], return false; } condition_number_2 = largest_singular / smallest_singular; - if (!std::isfinite(condition_number_2) - || condition_number_2 >= max_condition_number) + if (!std::isfinite(condition_number_2) || condition_number_2 >= max_condition_number) { message = "cell condition number is not below the configured maximum"; return false; @@ -160,13 +154,11 @@ void compute_cell_inverse(const double scale, for (int singular = 0; singular < kMatrixDimension; ++singular) { const long double sigma = singular_values[singular]; - inverse_value - += static_cast(right_vectors[row * kMatrixDimension + singular]) - * orthogonal_columns[column * kMatrixDimension + singular] - / (static_cast(scale) * sigma * sigma); + inverse_value += static_cast(right_vectors[row * kMatrixDimension + singular]) + * orthogonal_columns[column * kMatrixDimension + singular] + / (static_cast(scale) * sigma * sigma); } - computed_inverse[row * kMatrixDimension + column] - = static_cast(inverse_value); + computed_inverse[row * kMatrixDimension + column] = static_cast(inverse_value); } } } @@ -179,15 +171,12 @@ bool check_inverse_consistency(const Matrix9& cell_wire, double& inverse_residual, std::string& message) { - const double direct_inverse_residual - = received_inverse_residual(cell_wire, inverse_wire, false); - const double transposed_inverse_residual - = received_inverse_residual(cell_wire, inverse_wire, true); + const double direct_inverse_residual = received_inverse_residual(cell_wire, inverse_wire, false); + const double transposed_inverse_residual = received_inverse_residual(cell_wire, inverse_wire, true); inverse_residual = std::min(direct_inverse_residual, transposed_inverse_residual); const double residual_limit = inverse_absolute_tolerance - + inverse_relative_tolerance * condition_number_2 - * std::numeric_limits::epsilon(); + + inverse_relative_tolerance * condition_number_2 * std::numeric_limits::epsilon(); if (!std::isfinite(inverse_residual) || inverse_residual > residual_limit) { message = "received cell inverse is inconsistent with the cell"; @@ -199,9 +188,7 @@ bool check_inverse_consistency(const Matrix9& cell_wire, Matrix9 transpose_matrix9(const Matrix9& values) { - return {{values[0], values[3], values[6], - values[1], values[4], values[7], - values[2], values[5], values[8]}}; + return {{values[0], values[3], values[6], values[1], values[4], values[7], values[2], values[5], values[8]}}; } CellValidation validate_ipi_cell(const Matrix9& cell_wire, @@ -238,22 +225,27 @@ CellValidation validate_ipi_cell(const Matrix9& cell_wire, double singular_values[kMatrixDimension]; Matrix9 orthogonal_columns; Matrix9 right_vectors; - if (!compute_cell_svd(scaled_cell, singular_values, orthogonal_columns, right_vectors, - result.message)) + if (!compute_cell_svd(scaled_cell, singular_values, orthogonal_columns, right_vectors, result.message)) { return result; } - if (!compute_condition_number(singular_values, max_condition_number, - result.condition_number_2, result.message)) + if (!compute_condition_number(singular_values, max_condition_number, result.condition_number_2, result.message)) { return result; } - compute_cell_inverse(scale, singular_values, orthogonal_columns, right_vectors, + compute_cell_inverse(scale, + singular_values, + orthogonal_columns, + right_vectors, result.computed_inverse_wire_bohr_inv); - if (!check_inverse_consistency(cell_wire, inverse_wire, result.condition_number_2, - inverse_absolute_tolerance, inverse_relative_tolerance, - result.inverse_residual, result.message)) + if (!check_inverse_consistency(cell_wire, + inverse_wire, + result.condition_number_2, + inverse_absolute_tolerance, + inverse_relative_tolerance, + result.inverse_residual, + result.message)) { return result; } @@ -262,9 +254,7 @@ CellValidation validate_ipi_cell(const Matrix9& cell_wire, return result; } -bool validate_positions(const std::vector& positions_bohr, - std::size_t coordinate_count, - std::string& message) +bool validate_positions(const std::vector& positions_bohr, std::size_t coordinate_count, std::string& message) { if (positions_bohr.size() != coordinate_count) { @@ -337,10 +327,8 @@ bool validate_virial_inputs(const Matrix9& stress_ry_per_bohr3, message = "cell volume must be finite and positive"; return false; } - if (!std::isfinite(antisymmetric_absolute_tolerance) - || antisymmetric_absolute_tolerance < 0.0 - || !std::isfinite(antisymmetric_relative_tolerance) - || antisymmetric_relative_tolerance < 0.0) + if (!std::isfinite(antisymmetric_absolute_tolerance) || antisymmetric_absolute_tolerance < 0.0 + || !std::isfinite(antisymmetric_relative_tolerance) || antisymmetric_relative_tolerance < 0.0) { message = "stress symmetry tolerances must be finite and nonnegative"; return false; @@ -363,17 +351,13 @@ bool check_stress_symmetry(const Matrix9& stress_ry_per_bohr3, { for (int column = row + 1; column < kMatrixDimension; ++column) { - const double difference - = std::fabs(stress_ry_per_bohr3[row * kMatrixDimension + column] - - stress_ry_per_bohr3[column * kMatrixDimension + row]); - max_antisymmetric_component - = std::max(max_antisymmetric_component, difference); + const double difference = std::fabs(stress_ry_per_bohr3[row * kMatrixDimension + column] + - stress_ry_per_bohr3[column * kMatrixDimension + row]); + max_antisymmetric_component = std::max(max_antisymmetric_component, difference); } } - const double symmetry_limit - = antisymmetric_absolute_tolerance + antisymmetric_relative_tolerance * maximum_stress; - if (!std::isfinite(max_antisymmetric_component) - || max_antisymmetric_component > symmetry_limit) + const double symmetry_limit = antisymmetric_absolute_tolerance + antisymmetric_relative_tolerance * maximum_stress; + if (!std::isfinite(max_antisymmetric_component) || max_antisymmetric_component > symmetry_limit) { message = "stress tensor is not symmetric within tolerance"; return false; @@ -397,8 +381,7 @@ bool compute_symmetric_virial(const Matrix9& stress_ry_per_bohr3, + stress_ry_per_bohr3[column * kMatrixDimension + row]); const long double converted = 0.5L * volume_bohr3 * symmetric_stress; if (!std::isfinite(converted) - || std::fabs(converted) - > static_cast(std::numeric_limits::max())) + || std::fabs(converted) > static_cast(std::numeric_limits::max())) { message = "converted virial is not representable as finite doubles"; return false; @@ -412,9 +395,9 @@ bool compute_symmetric_virial(const Matrix9& stress_ry_per_bohr3, } // namespace VirialConversion make_ipi_virial(const Matrix9& stress_ry_per_bohr3, - double volume_bohr3, - double antisymmetric_absolute_tolerance, - double antisymmetric_relative_tolerance) + double volume_bohr3, + double antisymmetric_absolute_tolerance, + double antisymmetric_relative_tolerance) { VirialConversion result = make_failed_virial_conversion(); @@ -434,10 +417,7 @@ VirialConversion make_ipi_virial(const Matrix9& stress_ry_per_bohr3, { return result; } - if (!compute_symmetric_virial(stress_ry_per_bohr3, - volume_bohr3, - result.wire_virial_hartree, - result.message)) + if (!compute_symmetric_virial(stress_ry_per_bohr3, volume_bohr3, result.wire_virial_hartree, result.message)) { return result; } diff --git a/source/source_relax/socket_frame.h b/source/source_relax/socket_frame.h index 759a4ee13a2..bee302af79c 100644 --- a/source/source_relax/socket_frame.h +++ b/source/source_relax/socket_frame.h @@ -35,17 +35,15 @@ CellValidation validate_ipi_cell(const Matrix9& cell_wire, double max_condition_number, double inverse_absolute_tolerance, double inverse_relative_tolerance); -bool validate_positions(const std::vector& positions_bohr, - std::size_t coordinate_count, - std::string& message); +bool validate_positions(const std::vector& positions_bohr, std::size_t coordinate_count, std::string& message); bool checked_position_count(std::int32_t nat_socket, int nat_expected, std::size_t& coordinate_count, std::string& message); VirialConversion make_ipi_virial(const Matrix9& stress_ry_per_bohr3, - double volume_bohr3, - double antisymmetric_absolute_tolerance, - double antisymmetric_relative_tolerance); + double volume_bohr3, + double antisymmetric_absolute_tolerance, + double antisymmetric_relative_tolerance); } // namespace SocketFrame #endif diff --git a/source/source_relax/socket_frame_utils.cpp b/source/source_relax/socket_frame_utils.cpp index d4507572c2e..e048ab09ba8 100644 --- a/source/source_relax/socket_frame_utils.cpp +++ b/source/source_relax/socket_frame_utils.cpp @@ -75,9 +75,7 @@ void rotate_columns(SocketFrame::Matrix9& values, int first, int second, double bool one_sided_jacobi(SocketFrame::Matrix9& columns, SocketFrame::Matrix9& right_vectors) { - right_vectors = {{1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0}}; + right_vectors = {{1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}}; const double multiplier = 32.0 * std::numeric_limits::epsilon(); const int pairs[3][2] = {{0, 1}, {0, 2}, {1, 2}}; @@ -97,8 +95,7 @@ bool one_sided_jacobi(SocketFrame::Matrix9& columns, SocketFrame::Matrix9& right } const double tau = (second_norm - first_norm) / (2.0 * dot); - const double tangent - = std::copysign(1.0 / (std::fabs(tau) + std::hypot(1.0, tau)), tau); + const double tangent = std::copysign(1.0 / (std::fabs(tau) + std::hypot(1.0, tau)), tau); const double cosine = 1.0 / std::sqrt(1.0 + tangent * tangent); const double sine = tangent * cosine; rotate_columns(columns, first, second, cosine, sine); @@ -124,9 +121,7 @@ long double scaled_determinant(const SocketFrame::Matrix9& values) const long double a20 = values[6]; const long double a21 = values[7]; const long double a22 = values[8]; - return a00 * (a11 * a22 - a12 * a21) - - a01 * (a10 * a22 - a12 * a20) - + a02 * (a10 * a21 - a11 * a20); + return a00 * (a11 * a22 - a12 * a21) - a01 * (a10 * a22 - a12 * a20) + a02 * (a10 * a21 - a11 * a20); } double received_inverse_residual(const SocketFrame::Matrix9& cell, @@ -141,11 +136,9 @@ double received_inverse_residual(const SocketFrame::Matrix9& cell, long double product = 0.0L; for (int inner = 0; inner < kMatrixDimension; ++inner) { - const int inverse_index = transpose_inverse - ? column * kMatrixDimension + inner - : inner * kMatrixDimension + column; - product += static_cast(cell[row * kMatrixDimension + inner]) - * inverse[inverse_index]; + const int inverse_index + = transpose_inverse ? column * kMatrixDimension + inner : inner * kMatrixDimension + column; + product += static_cast(cell[row * kMatrixDimension + inner]) * inverse[inverse_index]; } const long double expected = row == column ? 1.0L : 0.0L; maximum = std::max(maximum, std::fabs(product - expected)); diff --git a/source/source_relax/socket_handlers.cpp b/source/source_relax/socket_handlers.cpp index 94fa504beb7..71ceb4c4eab 100644 --- a/source/source_relax/socket_handlers.cpp +++ b/source/source_relax/socket_handlers.cpp @@ -1,10 +1,10 @@ #include "socket_handlers.h" -#include "source_relax/socket_frame.h" #include "source_base/global_function.h" #include "source_cell/unitcell.h" #include "source_esolver/esolver.h" #include "source_io/module_parameter/input_parameter.h" +#include "source_relax/socket_frame.h" #ifdef __JSON #include "source_io/module_json/output_info.h" #endif @@ -19,8 +19,16 @@ namespace SocketHandlers { +using SocketUtils::all_ranks_converged; +using SocketUtils::bcast_double_vector; +using SocketUtils::bcast_header; +using SocketUtils::bcast_int32; using SocketUtils::ComputedFrame; using SocketUtils::DriverState; +using SocketUtils::fail_during_collective_stage; +using SocketUtils::flatten_forces_hartree_per_bohr; +using SocketUtils::ipi_cell_bohr_from_unitcell; +using SocketUtils::is_root; using SocketUtils::kInverseAbsoluteTolerance; using SocketUtils::kInverseRelativeTolerance; using SocketUtils::kMaxCellCondition; @@ -28,14 +36,6 @@ using SocketUtils::kMaxInitBytes; using SocketUtils::kRyToHartree; using SocketUtils::kStressAbsoluteTolerance; using SocketUtils::kStressRelativeTolerance; -using SocketUtils::all_ranks_converged; -using SocketUtils::bcast_double_vector; -using SocketUtils::bcast_header; -using SocketUtils::bcast_int32; -using SocketUtils::fail_during_collective_stage; -using SocketUtils::flatten_forces_hartree_per_bohr; -using SocketUtils::ipi_cell_bohr_from_unitcell; -using SocketUtils::is_root; using SocketUtils::matrix9_from_stress; using SocketUtils::max_abs_delta; using SocketUtils::max_wrapped_direct_delta_from_unitcell; @@ -105,9 +105,7 @@ void handle_status(IpiSocket& socket, const DriverState state) quit_if_root_failed(io_failed, io_message); } -void handle_init(IpiSocket& socket, - DriverContext& context, - std::ofstream& ofs_running) +void handle_init(IpiSocket& socket, DriverContext& context, std::ofstream& ofs_running) { std::int32_t rid = 0; std::int32_t nbytes = 0; @@ -223,9 +221,7 @@ PosdataPayload read_posdata(IpiSocket& socket, const UnitCell& ucell, const Driv if (io_failed == 0) { payload.positions = socket.read_doubles(coordinate_count); - if (!SocketFrame::validate_positions(payload.positions, - coordinate_count, - io_message)) + if (!SocketFrame::validate_positions(payload.positions, coordinate_count, io_message)) { io_failed = 1; } @@ -263,29 +259,23 @@ void bcast_posdata(PosdataPayload& payload) bcast_double_vector(payload.positions); } -void check_posdata_geometry(const DriverContext& context, - const PosdataPayload& payload, - DriverContext& mutable_context) +void check_posdata_geometry(const DriverContext& context, const PosdataPayload& payload, DriverContext& mutable_context) { - const double max_cell_delta_bohr = max_abs_delta( - std::vector(payload.cell.begin(), payload.cell.end()), - context.reference_cell); + const double max_cell_delta_bohr + = max_abs_delta(std::vector(payload.cell.begin(), payload.cell.end()), context.reference_cell); if (max_cell_delta_bohr > unchanged_cell_tolerance(payload.cell)) { - ModuleBase::WARNING_QUIT("ABACUS socket", - "variable-cell socket updates are not supported yet."); + ModuleBase::WARNING_QUIT("ABACUS socket", "variable-cell socket updates are not supported yet."); } if (!mutable_context.checked_initial_positions) { mutable_context.checked_initial_positions = true; - if (max_wrapped_direct_delta_from_unitcell(*context.ucell, payload.positions) > 1.0e-5 - && is_root()) + if (max_wrapped_direct_delta_from_unitcell(*context.ucell, payload.positions) > 1.0e-5 && is_root()) { - ModuleBase::WARNING( - "ABACUS socket", - "first POSDATA positions are not PBC-equivalent to STRU atom order; " - "i-PI POSDATA carries no species, so the client atoms should use the same atom order " - "as STRU."); + ModuleBase::WARNING("ABACUS socket", + "first POSDATA positions are not PBC-equivalent to STRU atom order; " + "i-PI POSDATA carries no species, so the client atoms should use the same atom order " + "as STRU."); } } } @@ -321,9 +311,7 @@ void run_esolver_for_positions(UnitCell& ucell, } } -void compute_energy_hartree(ModuleESolver::ESolver* esolver, - ComputedFrame& computed, - std::ofstream& ofs_running) +void compute_energy_hartree(ModuleESolver::ESolver* esolver, ComputedFrame& computed, std::ofstream& ofs_running) { double energy_ry = 0.0; try @@ -339,8 +327,7 @@ void compute_energy_hartree(ModuleESolver::ESolver* esolver, fail_during_collective_stage("cal_energy", "unknown socket energy failure"); } const int local_failed = std::isfinite(energy_ry) ? 0 : 1; - throw_if_any_rank_failed(local_failed, - local_failed == 0 ? "" : "socket energy is not finite"); + throw_if_any_rank_failed(local_failed, local_failed == 0 ? "" : "socket energy is not finite"); if (!std::isfinite(energy_ry)) { ModuleBase::WARNING_QUIT("ABACUS socket", "socket energy is not finite."); @@ -348,16 +335,12 @@ void compute_energy_hartree(ModuleESolver::ESolver* esolver, computed.energy_hartree = energy_ry * kRyToHartree; if (is_root()) { - ofs_running << " ABACUS socket return energy " - << energy_ry << " Ry, " - << energy_ry * ModuleBase::Ry_to_eV << " eV, " - << computed.energy_hartree << " Ha" << std::endl; + ofs_running << " ABACUS socket return energy " << energy_ry << " Ry, " << energy_ry * ModuleBase::Ry_to_eV + << " eV, " << computed.energy_hartree << " Ha" << std::endl; } } -void compute_forces(UnitCell& ucell, - ModuleESolver::ESolver* esolver, - ComputedFrame& computed) +void compute_forces(UnitCell& ucell, ModuleESolver::ESolver* esolver, ComputedFrame& computed) { ModuleBase::matrix force; try @@ -392,9 +375,7 @@ void compute_forces(UnitCell& ucell, computed.forces_present = true; } -void compute_stress(UnitCell& ucell, - ModuleESolver::ESolver* esolver, - ComputedFrame& computed) +void compute_stress(UnitCell& ucell, ModuleESolver::ESolver* esolver, ComputedFrame& computed) { ModuleBase::matrix stress; try @@ -413,11 +394,10 @@ void compute_stress(UnitCell& ucell, std::string local_message; try { - const SocketFrame::VirialConversion virial - = SocketFrame::make_ipi_virial(matrix9_from_stress(stress), - ucell.omega, - kStressAbsoluteTolerance, - kStressRelativeTolerance); + const SocketFrame::VirialConversion virial = SocketFrame::make_ipi_virial(matrix9_from_stress(stress), + ucell.omega, + kStressAbsoluteTolerance, + kStressRelativeTolerance); if (!virial.ok) { throw std::runtime_error(virial.message); @@ -439,9 +419,7 @@ void compute_stress(UnitCell& ucell, } } // namespace -void handle_posdata(IpiSocket& socket, - DriverContext& context, - std::ofstream& ofs_running) +void handle_posdata(IpiSocket& socket, DriverContext& context, std::ofstream& ofs_running) { PosdataPayload payload = read_posdata(socket, *context.ucell, context.state); bcast_posdata(payload); @@ -456,9 +434,8 @@ void handle_posdata(IpiSocket& socket, computed.scf_converged = all_ranks_converged(context.esolver->conv_esolver); if (!computed.scf_converged && is_root()) { - ModuleBase::WARNING( - "ABACUS socket", - "SCF did not converge; returning the available frame and marking it in i-PI extras."); + ModuleBase::WARNING("ABACUS socket", + "SCF did not converge; returning the available frame and marking it in i-PI extras."); } compute_energy_hartree(context.esolver, computed, ofs_running); if (context.inp->cal_force) diff --git a/source/source_relax/socket_handlers.h b/source/source_relax/socket_handlers.h index 4681491733e..f73224503ca 100644 --- a/source/source_relax/socket_handlers.h +++ b/source/source_relax/socket_handlers.h @@ -1,8 +1,8 @@ #ifndef ABACUS_SOURCE_RELAX_SOCKET_HANDLERS_H #define ABACUS_SOURCE_RELAX_SOCKET_HANDLERS_H -#include "source_relax/socket_utils.h" #include "source_relax/socket_ipi.h" +#include "source_relax/socket_utils.h" #include #include @@ -38,13 +38,9 @@ std::string read_header_bcast(IpiSocket& socket, const SocketUtils::DriverState void handle_status(IpiSocket& socket, const SocketUtils::DriverState state); -void handle_init(IpiSocket& socket, - DriverContext& context, - std::ofstream& ofs_running); +void handle_init(IpiSocket& socket, DriverContext& context, std::ofstream& ofs_running); -void handle_posdata(IpiSocket& socket, - DriverContext& context, - std::ofstream& ofs_running); +void handle_posdata(IpiSocket& socket, DriverContext& context, std::ofstream& ofs_running); void handle_getforce(IpiSocket& socket, DriverContext& context); diff --git a/source/source_relax/socket_ipi.cpp b/source/source_relax/socket_ipi.cpp index e84d8637ba2..521a1b96789 100644 --- a/source/source_relax/socket_ipi.cpp +++ b/source/source_relax/socket_ipi.cpp @@ -13,8 +13,7 @@ static_assert(sizeof(std::int32_t) == 4, "i-PI requires a 4-byte integer"); static_assert(sizeof(double) == 8, "i-PI requires an 8-byte float"); -static_assert(std::numeric_limits::is_iec559, - "i-PI requires IEEE-754 double precision"); +static_assert(std::numeric_limits::is_iec559, "i-PI requires IEEE-754 double precision"); namespace { diff --git a/source/source_relax/socket_utils.cpp b/source/source_relax/socket_utils.cpp index ccdfa38b369..8f41d005460 100644 --- a/source/source_relax/socket_utils.cpp +++ b/source/source_relax/socket_utils.cpp @@ -50,17 +50,12 @@ void throw_if_any_rank_failed(int local_failed, std::string local_message) /// no collective: the exception propagates to Socket_Driver::socket_driver, /// which reports a clean WARNING_QUIT into the log. Do not replace the /// serial throw with abort(); serial builds are user-facing. -[[noreturn]] void fail_during_collective_stage(const char* stage, - const std::string& message) +[[noreturn]] void fail_during_collective_stage(const char* stage, const std::string& message) { #ifdef __MPI int rank = -1; MPI_Comm_rank(MPI_COMM_WORLD, &rank); - std::fprintf(stderr, - "ABACUS_SOCKET_MPI_FATAL stage=%s rank=%d message=%s\n", - stage, - rank, - message.c_str()); + std::fprintf(stderr, "ABACUS_SOCKET_MPI_FATAL stage=%s rank=%d message=%s\n", stage, rank, message.c_str()); std::fflush(stderr); MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE); std::abort(); @@ -82,8 +77,7 @@ std::string properties_extra(const ComputedFrame& frame) { extra << ",\"stress\""; } - extra << "],\"scf_converged\":" - << (frame.scf_converged ? "true" : "false") << "}"; + extra << "],\"scf_converged\":" << (frame.scf_converged ? "true" : "false") << "}"; return extra.str(); } @@ -167,9 +161,15 @@ std::vector ipi_cell_bohr_from_unitcell(const UnitCell& ucell) // ASE/i-PI sends POSDATA cell as cell.T in C order. ABACUS stores // lattice vectors as rows in latvec, so use the transposed order here. return { - ucell.latvec.e11 * lat0, ucell.latvec.e21 * lat0, ucell.latvec.e31 * lat0, - ucell.latvec.e12 * lat0, ucell.latvec.e22 * lat0, ucell.latvec.e32 * lat0, - ucell.latvec.e13 * lat0, ucell.latvec.e23 * lat0, ucell.latvec.e33 * lat0, + ucell.latvec.e11 * lat0, + ucell.latvec.e21 * lat0, + ucell.latvec.e31 * lat0, + ucell.latvec.e12 * lat0, + ucell.latvec.e22 * lat0, + ucell.latvec.e32 * lat0, + ucell.latvec.e13 * lat0, + ucell.latvec.e23 * lat0, + ucell.latvec.e33 * lat0, }; } diff --git a/source/source_relax/socket_utils.h b/source/source_relax/socket_utils.h index ccac9761967..4314d165302 100644 --- a/source/source_relax/socket_utils.h +++ b/source/source_relax/socket_utils.h @@ -1,8 +1,8 @@ #ifndef ABACUS_SOURCE_RELAX_SOCKET_UTILS_H #define ABACUS_SOURCE_RELAX_SOCKET_UTILS_H -#include "source_relax/socket_frame.h" #include "source_base/matrix.h" +#include "source_relax/socket_frame.h" #include #include @@ -16,8 +16,7 @@ namespace SocketUtils constexpr double kRyToHartree = 0.5; constexpr int kIpiRankRoot = 0; constexpr double kMaxCellCondition = 1.0e12; -constexpr double kInverseAbsoluteTolerance - = 64.0 * std::numeric_limits::epsilon(); +constexpr double kInverseAbsoluteTolerance = 64.0 * std::numeric_limits::epsilon(); constexpr double kInverseRelativeTolerance = 64.0; constexpr double kStressAbsoluteTolerance = 1.0e-10; constexpr double kStressRelativeTolerance = 1.0e-8; @@ -43,8 +42,7 @@ struct ComputedFrame bool all_ranks_converged(bool local_converged); void throw_if_any_rank_failed(int local_failed, std::string local_message); -[[noreturn]] void fail_during_collective_stage(const char* stage, - const std::string& message); +[[noreturn]] void fail_during_collective_stage(const char* stage, const std::string& message); std::string properties_extra(const ComputedFrame& frame); bool is_root(); void bcast_double_vector(std::vector& values); @@ -56,14 +54,11 @@ void quit_if_root_failed(int root_failed, std::string root_message); std::string bcast_header(std::string header); std::string address(); std::vector ipi_cell_bohr_from_unitcell(const UnitCell& ucell); -double max_wrapped_direct_delta_from_unitcell(const UnitCell& ucell, - const std::vector& positions_bohr); +double max_wrapped_direct_delta_from_unitcell(const UnitCell& ucell, const std::vector& positions_bohr); double max_abs_delta(const std::vector& a, const std::vector& b); double unchanged_cell_tolerance(const SocketFrame::Matrix9& cell); -void set_positions_from_ipi_bohr(UnitCell& ucell, - const std::vector& positions_bohr); -std::vector flatten_forces_hartree_per_bohr(const ModuleBase::matrix& force, - int nat); +void set_positions_from_ipi_bohr(UnitCell& ucell, const std::vector& positions_bohr); +std::vector flatten_forces_hartree_per_bohr(const ModuleBase::matrix& force, int nat); SocketFrame::Matrix9 matrix9_from_stress(const ModuleBase::matrix& stress); std::vector vector_from_matrix9(const SocketFrame::Matrix9& values); } // namespace SocketUtils diff --git a/source/source_relax/unittests/for_test.h b/source/source_relax/unittests/for_test.h index 698f57f3d4f..521dfcf83e9 100644 --- a/source/source_relax/unittests/for_test.h +++ b/source/source_relax/unittests/for_test.h @@ -99,12 +99,20 @@ pseudo::pseudo() pseudo::~pseudo() { } -SepPot::SepPot(){} -SepPot::~SepPot(){} -Sep_Cell::Sep_Cell() noexcept {} -Sep_Cell::~Sep_Cell() noexcept {} +SepPot::SepPot() +{ +} +SepPot::~SepPot() +{ +} +Sep_Cell::Sep_Cell() noexcept +{ +} +Sep_Cell::~Sep_Cell() noexcept +{ +} int ModuleSymmetry::Symmetry::symm_flag = 0; -void ModuleSymmetry::Symmetry::symmetrize_mat3(ModuleBase::matrix& sigma, const Lattice& lat)const {}; -void ModuleSymmetry::Symmetry::symmetrize_vec3_nat(double* v)const {}; +void ModuleSymmetry::Symmetry::symmetrize_mat3(ModuleBase::matrix& sigma, const Lattice& lat) const {}; +void ModuleSymmetry::Symmetry::symmetrize_vec3_nat(double* v) const {}; #endif diff --git a/source/source_relax/unittests/mock_remake_cell.cpp b/source/source_relax/unittests/mock_remake_cell.cpp index ee51ba8653b..06fe5aef274 100644 --- a/source/source_relax/unittests/mock_remake_cell.cpp +++ b/source/source_relax/unittests/mock_remake_cell.cpp @@ -1,58 +1,58 @@ #include "mock_remake_cell.h" + #include namespace unitcell { - // Mock state variables - bool remake_cell_called = false; - std::string remake_cell_latName = ""; - ModuleBase::Matrix3 remake_cell_latvec; - - void remake_cell(Lattice& lat) - { - remake_cell_called = true; - remake_cell_latName = lat.latName; - remake_cell_latvec = lat.latvec; +// Mock state variables +bool remake_cell_called = false; +std::string remake_cell_latName = ""; +ModuleBase::Matrix3 remake_cell_latvec; - // Mock implementation: enforce simple cubic structure for "sc" - if (lat.latName == "sc") - { - double celldm = std::sqrt(lat.latvec.e11 * lat.latvec.e11 + - lat.latvec.e12 * lat.latvec.e12 + - lat.latvec.e13 * lat.latvec.e13); - lat.latvec.Zero(); - lat.latvec.e11 = celldm; - lat.latvec.e22 = celldm; - lat.latvec.e33 = celldm; - } - // Mock implementation: enforce FCC structure for "fcc" - else if (lat.latName == "fcc") - { - double celldm = std::sqrt(lat.latvec.e11 * lat.latvec.e11 + - lat.latvec.e12 * lat.latvec.e12 + - lat.latvec.e13 * lat.latvec.e13) / std::sqrt(2.0); - lat.latvec.e11 = -celldm; - lat.latvec.e12 = 0.0; - lat.latvec.e13 = celldm; - lat.latvec.e21 = 0.0; - lat.latvec.e22 = celldm; - lat.latvec.e23 = celldm; - lat.latvec.e31 = -celldm; - lat.latvec.e32 = celldm; - lat.latvec.e33 = 0.0; - } - // For other lattice types, do nothing (just track the call) - } +void remake_cell(Lattice& lat) +{ + remake_cell_called = true; + remake_cell_latName = lat.latName; + remake_cell_latvec = lat.latvec; - void reset_remake_cell_mock() + // Mock implementation: enforce simple cubic structure for "sc" + if (lat.latName == "sc") { - remake_cell_called = false; - remake_cell_latName = ""; - remake_cell_latvec.Zero(); + double celldm = std::sqrt(lat.latvec.e11 * lat.latvec.e11 + lat.latvec.e12 * lat.latvec.e12 + + lat.latvec.e13 * lat.latvec.e13); + lat.latvec.Zero(); + lat.latvec.e11 = celldm; + lat.latvec.e22 = celldm; + lat.latvec.e33 = celldm; } - - bool was_remake_cell_called() + // Mock implementation: enforce FCC structure for "fcc" + else if (lat.latName == "fcc") { - return remake_cell_called; + double celldm = std::sqrt(lat.latvec.e11 * lat.latvec.e11 + lat.latvec.e12 * lat.latvec.e12 + + lat.latvec.e13 * lat.latvec.e13) + / std::sqrt(2.0); + lat.latvec.e11 = -celldm; + lat.latvec.e12 = 0.0; + lat.latvec.e13 = celldm; + lat.latvec.e21 = 0.0; + lat.latvec.e22 = celldm; + lat.latvec.e23 = celldm; + lat.latvec.e31 = -celldm; + lat.latvec.e32 = celldm; + lat.latvec.e33 = 0.0; } + // For other lattice types, do nothing (just track the call) +} + +void reset_remake_cell_mock() +{ + remake_cell_called = false; + remake_cell_latName = ""; + remake_cell_latvec.Zero(); +} + +bool was_remake_cell_called() +{ + return remake_cell_called; } +} // namespace unitcell diff --git a/source/source_relax/unittests/mock_remake_cell.h b/source/source_relax/unittests/mock_remake_cell.h index c2d4d63bd81..c09774effd1 100644 --- a/source/source_relax/unittests/mock_remake_cell.h +++ b/source/source_relax/unittests/mock_remake_cell.h @@ -2,21 +2,22 @@ #define MOCK_REMAKE_CELL_H #include "source_cell/unitcell_data.h" + #include namespace unitcell { - // Mock state tracking - extern bool remake_cell_called; - extern std::string remake_cell_latName; - extern ModuleBase::Matrix3 remake_cell_latvec; +// Mock state tracking +extern bool remake_cell_called; +extern std::string remake_cell_latName; +extern ModuleBase::Matrix3 remake_cell_latvec; - // Mock implementation of remake_cell - void remake_cell(Lattice& lat); +// Mock implementation of remake_cell +void remake_cell(Lattice& lat); - // Helper functions for testing - void reset_remake_cell_mock(); - bool was_remake_cell_called(); -} +// Helper functions for testing +void reset_remake_cell_mock(); +bool was_remake_cell_called(); +} // namespace unitcell #endif // MOCK_REMAKE_CELL_H diff --git a/source/source_relax/unittests/relax_test.h b/source/source_relax/unittests/relax_test.h index 8e83d7cc252..c35ce0dffe2 100644 --- a/source/source_relax/unittests/relax_test.h +++ b/source/source_relax/unittests/relax_test.h @@ -1,22 +1,30 @@ -#include "source_cell/unitcell.h" #include "source_base/matrix.h" +#include "source_cell/unitcell.h" -UnitCell::UnitCell(){}; -UnitCell::~UnitCell(){}; +UnitCell::UnitCell() {}; +UnitCell::~UnitCell() {}; -Magnetism::Magnetism(){}; -Magnetism::~Magnetism(){}; +Magnetism::Magnetism() {}; +Magnetism::~Magnetism() {}; -Atom::Atom(){}; -Atom::~Atom(){}; -Atom_pseudo::Atom_pseudo(){}; -Atom_pseudo::~Atom_pseudo(){}; -pseudo::pseudo(){}; -pseudo::~pseudo(){}; -SepPot::SepPot(){} -SepPot::~SepPot(){} -Sep_Cell::Sep_Cell() noexcept {} -Sep_Cell::~Sep_Cell() noexcept {} +Atom::Atom() {}; +Atom::~Atom() {}; +Atom_pseudo::Atom_pseudo() {}; +Atom_pseudo::~Atom_pseudo() {}; +pseudo::pseudo() {}; +pseudo::~pseudo() {}; +SepPot::SepPot() +{ +} +SepPot::~SepPot() +{ +} +Sep_Cell::Sep_Cell() noexcept +{ +} +Sep_Cell::~Sep_Cell() noexcept +{ +} int ModuleSymmetry::Symmetry::symm_flag = 0; -void ModuleSymmetry::Symmetry::symmetrize_mat3(ModuleBase::matrix& sigma, const Lattice& lat)const {}; -void ModuleSymmetry::Symmetry::symmetrize_vec3_nat(double* v)const {}; +void ModuleSymmetry::Symmetry::symmetrize_mat3(ModuleBase::matrix& sigma, const Lattice& lat) const {}; +void ModuleSymmetry::Symmetry::symmetrize_vec3_nat(double* v) const {}; diff --git a/source/source_relax/unittests/test_bfgs_basic.cpp b/source/source_relax/unittests/test_bfgs_basic.cpp index 3df6a3f5a0b..e849c1c8d21 100644 --- a/source/source_relax/unittests/test_bfgs_basic.cpp +++ b/source/source_relax/unittests/test_bfgs_basic.cpp @@ -1,9 +1,10 @@ +#include "source_io/module_parameter/parameter.h" +#include "source_relax/bfgs_basic.h" #include "source_relax/ions_move_basic.h" #include "source_relax/relax_data.h" + #include "gmock/gmock.h" -#include "source_io/module_parameter/parameter.h" #include "gtest/gtest.h" -#include "source_relax/bfgs_basic.h" /************************************************ * unit tests of class BFGS_Basic ***********************************************/ @@ -260,7 +261,9 @@ TEST_F(BFGSBasicTest, NewStepWarningQuit) std::vector etot_info(2, 0.0); testing::internal::CaptureStdout(); - EXPECT_EXIT(bfgs.new_step_for_testing(lat0, update_iter, ofs, etot_info, test_relax_method), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(bfgs.new_step_for_testing(lat0, update_iter, ofs, etot_info, test_relax_method), + ::testing::ExitedWithCode(1), + ""); std::string output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("bfgs_ndim > 1 not implemented yet")); } @@ -348,7 +351,9 @@ TEST_F(BFGSBasicTest, ComputeTrustRadiusWarningQuit) std::ofstream ofs("test_log.log"); testing::internal::CaptureStdout(); - EXPECT_EXIT(bfgs.compute_trust_radius_for_testing(ofs, etot_info, test_relax_method), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(bfgs.compute_trust_radius_for_testing(ofs, etot_info, test_relax_method), + ::testing::ExitedWithCode(1), + ""); std::string output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("bfgs history already reset at previous step, we got trapped!")); } \ No newline at end of file diff --git a/source/source_relax/unittests/test_cg_base.cpp b/source/source_relax/unittests/test_cg_base.cpp index 72223537d09..16ed352d98f 100644 --- a/source/source_relax/unittests/test_cg_base.cpp +++ b/source/source_relax/unittests/test_cg_base.cpp @@ -1,8 +1,7 @@ #include "../cg_base.h" -#include - #include +#include // Unit tests for CG_Base. normalize() is protected, so it is exercised // through a test subclass that re-exposes it (instead of a @@ -170,11 +169,11 @@ TEST(CGBase, BrentSameSignBranch) cg.Brent(fa, fb, fc, xa, xb, xc, best_x, xpt); const double dmove = (0.5 * 1.0 - 0.0 * 0.5) / (1.0 - 0.5); // = 1.0 - EXPECT_DOUBLE_EQ(best_x, dmove); // xpt was 0 + EXPECT_DOUBLE_EQ(best_x, dmove); // xpt was 0 EXPECT_DOUBLE_EQ(xpt, dmove); EXPECT_DOUBLE_EQ(xc, dmove); - EXPECT_DOUBLE_EQ(xb, 0.5); // old xc - EXPECT_DOUBLE_EQ(fb, 0.5); // old fc + EXPECT_DOUBLE_EQ(xb, 0.5); // old xc + EXPECT_DOUBLE_EQ(fb, 0.5); // old fc } TEST(CGBase, BrentOppositeSignQuadratic) @@ -194,7 +193,7 @@ TEST(CGBase, BrentOppositeSignQuadratic) cg.Brent(fa, fb, fc, xa, xb, xc, best_x, xpt); - EXPECT_DOUBLE_EQ(xpt, best_x); // xpt = dmove = best_x + old xpt(0) + EXPECT_DOUBLE_EQ(xpt, best_x); // xpt = dmove = best_x + old xpt(0) EXPECT_DOUBLE_EQ(xc, xpt); } diff --git a/source/source_relax/unittests/test_ions_move_basic.cpp b/source/source_relax/unittests/test_ions_move_basic.cpp index a82e80672da..9a1f1ad0371 100644 --- a/source/source_relax/unittests/test_ions_move_basic.cpp +++ b/source/source_relax/unittests/test_ions_move_basic.cpp @@ -1,8 +1,9 @@ -#include "gmock/gmock.h" -#include "gtest/gtest.h" +#include "for_test.h" #include "source_relax/ions_move_basic.h" #include "source_relax/relax_data.h" -#include "for_test.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" /************************************************ * unit tests of namespace Ions_Move_Basic @@ -85,7 +86,7 @@ TEST_F(IonsMoveBasicTest, MoveAtoms) ifs.close(); std::remove("test_move_atoms.log"); - EXPECT_THAT(output , ::testing::HasSubstr(expected_output)); + EXPECT_THAT(output, ::testing::HasSubstr(expected_output)); EXPECT_DOUBLE_EQ(pos[0], 0.0); EXPECT_DOUBLE_EQ(pos[1], 1.0); EXPECT_DOUBLE_EQ(pos[2], 2.0); @@ -102,7 +103,7 @@ TEST_F(IonsMoveBasicTest, CheckConvergedCase1) int update_iter = 1; const int test_relax_method = 1; const std::string out_level = "ie"; - const double force_thr = -1; // Input_para default; this test never set it + const double force_thr = -1; // Input_para default; this test never set it const double force_thr_ev = -1; // Input_para default; this test never set it std::vector etot_info(2, 0.0); for (int i = 0; i < Ions_Move_Basic::dim; ++i) @@ -113,7 +114,15 @@ TEST_F(IonsMoveBasicTest, CheckConvergedCase1) // Call the function being tested std::ofstream ofs("test_check_converged_case1.log"); testing::internal::CaptureStdout(); - bool converged = Ions_Move_Basic::check_converged(ucell, grad, update_iter, ofs, etot_info, force_thr, force_thr_ev, out_level, test_relax_method); + bool converged = Ions_Move_Basic::check_converged(ucell, + grad, + update_iter, + ofs, + etot_info, + force_thr, + force_thr_ev, + out_level, + test_relax_method); std::string std_outout = testing::internal::GetCapturedStdout(); ofs.close(); @@ -130,7 +139,7 @@ TEST_F(IonsMoveBasicTest, CheckConvergedCase1) "movement is possible.\n it may converged, otherwise no movement of atom is allowed.\n"; std::string expected_std = " ETOT DIFF (eV) : 0\n LARGEST GRAD (eV/Angstrom) : 0\n"; - EXPECT_THAT(ofs_output , ::testing::HasSubstr(expected_ofs)); + EXPECT_THAT(ofs_output, ::testing::HasSubstr(expected_ofs)); EXPECT_EQ(expected_std, std_outout); EXPECT_EQ(update_iter, 1); EXPECT_EQ(converged, true); @@ -153,7 +162,15 @@ TEST_F(IonsMoveBasicTest, CheckConvergedCase2) // Call the function being tested std::ofstream ofs("test_check_converged_case2.log"); testing::internal::CaptureStdout(); - bool converged = Ions_Move_Basic::check_converged(ucell, grad, update_iter, ofs, etot_info, force_thr, force_thr_ev, out_level, test_relax_method); + bool converged = Ions_Move_Basic::check_converged(ucell, + grad, + update_iter, + ofs, + etot_info, + force_thr, + force_thr_ev, + out_level, + test_relax_method); std::string std_outout = testing::internal::GetCapturedStdout(); ofs.close(); @@ -168,7 +185,7 @@ TEST_F(IonsMoveBasicTest, CheckConvergedCase2) " energy difference (ry) = 0\n largest gradient (ry/bohr) = 0.1\n"; std::string expected_std = " ETOT DIFF (eV) : 0\n LARGEST GRAD (eV/Angstrom) : 2.57111\n"; - EXPECT_THAT(ofs_output , ::testing::HasSubstr(expected_ofs)); + EXPECT_THAT(ofs_output, ::testing::HasSubstr(expected_ofs)); EXPECT_EQ(expected_std, std_outout); EXPECT_EQ(update_iter, 2); EXPECT_EQ(converged, true); @@ -191,7 +208,15 @@ TEST_F(IonsMoveBasicTest, CheckConvergedCase3) // Call the function being tested std::ofstream ofs("test_check_converged_case3.log"); testing::internal::CaptureStdout(); - bool converged = Ions_Move_Basic::check_converged(ucell, grad, update_iter, ofs, etot_info, force_thr, force_thr_ev, out_level, test_relax_method); + bool converged = Ions_Move_Basic::check_converged(ucell, + grad, + update_iter, + ofs, + etot_info, + force_thr, + force_thr_ev, + out_level, + test_relax_method); std::string std_outout = testing::internal::GetCapturedStdout(); ofs.close(); @@ -206,7 +231,7 @@ TEST_F(IonsMoveBasicTest, CheckConvergedCase3) " energy difference (ry) = 1\n largest gradient (ry/bohr) = 0.1\n"; std::string expected_std = " ETOT DIFF (eV) : 13.6057\n LARGEST GRAD (eV/Angstrom) : 2.57111\n"; - EXPECT_THAT(ofs_output , ::testing::HasSubstr(expected_ofs)); + EXPECT_THAT(ofs_output, ::testing::HasSubstr(expected_ofs)); EXPECT_EQ(expected_std, std_outout); EXPECT_EQ(update_iter, 1); EXPECT_EQ(converged, false); @@ -235,7 +260,7 @@ TEST_F(IonsMoveBasicTest, TerminateConverged) std::string expected_ofs = " end of geometry optimization\n istep = 2\n " " update iteration = 5\n"; - EXPECT_THAT(ofs_output , ::testing::HasSubstr(expected_ofs)); + EXPECT_THAT(ofs_output, ::testing::HasSubstr(expected_ofs)); } // Test the terminate() function when not converged @@ -257,9 +282,10 @@ TEST_F(IonsMoveBasicTest, TerminateNotConverged) ifs.close(); std::remove("test_terminate_not_converged.log"); - std::string expected_ofs = " the maximum number of steps has been reached.\n Relaxation is not converged after reaching relax_nmax!\n end of geometry optimization.\n"; + std::string expected_ofs = " the maximum number of steps has been reached.\n Relaxation is not converged after " + "reaching relax_nmax!\n end of geometry optimization.\n"; - EXPECT_THAT(ofs_output , ::testing::HasSubstr(expected_ofs)); + EXPECT_THAT(ofs_output, ::testing::HasSubstr(expected_ofs)); } // Test the setup_etot() function diff --git a/source/source_relax/unittests/test_ions_move_bfgs.cpp b/source/source_relax/unittests/test_ions_move_bfgs.cpp index cd634686a6c..0514e113bf8 100644 --- a/source/source_relax/unittests/test_ions_move_bfgs.cpp +++ b/source/source_relax/unittests/test_ions_move_bfgs.cpp @@ -1,10 +1,11 @@ -#include "source_relax/relax_criteria.h" #include "for_test.h" -#include "gtest/gtest.h" -#include "gmock/gmock.h" #include "source_io/module_parameter/parameter.h" #include "source_relax/ions_move_basic.h" #include "source_relax/ions_move_bfgs.h" +#include "source_relax/relax_criteria.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" /************************************************ * unit tests of class Ions_Move_BFGS @@ -103,8 +104,12 @@ TEST_F(IonsMoveBFGSTest, StartCase2) ucell.atoms[0].tau = std::vector>(2); ucell.atoms[0].taud = std::vector>(2); ucell.atoms[0].mbl = std::vector>(2, {1, 1, 1}); - ucell.atoms[0].tau[0].x = 0.0; ucell.atoms[0].tau[0].y = 0.0; ucell.atoms[0].tau[0].z = 0.0; - ucell.atoms[0].tau[1].x = 1.0; ucell.atoms[0].tau[1].y = 0.0; ucell.atoms[0].tau[1].z = 0.0; + ucell.atoms[0].tau[0].x = 0.0; + ucell.atoms[0].tau[0].y = 0.0; + ucell.atoms[0].tau[0].z = 0.0; + ucell.atoms[0].tau[1].x = 1.0; + ucell.atoms[0].tau[1].y = 0.0; + ucell.atoms[0].tau[1].z = 0.0; ucell.lat0 = 1.0; ucell.set_atom_flag = true; @@ -250,7 +255,13 @@ TEST_F(IonsMoveBFGSTest, BfgsRoutineCase1) // Call the function being tested std::ofstream ofs("test_bfgs_routine_case1.log"); testing::internal::CaptureStdout(); - bfgs.bfgs_routine_for_testing(lat0, istep, update_iter, ofs, etot_info, criteria.out_level, criteria.test_relax_method); + bfgs.bfgs_routine_for_testing(lat0, + istep, + update_iter, + ofs, + etot_info, + criteria.out_level, + criteria.test_relax_method); std::string std_outout = testing::internal::GetCapturedStdout(); ofs.close(); @@ -315,7 +326,13 @@ TEST_F(IonsMoveBFGSTest, BfgsRoutineCase2) // Call the function being tested std::ofstream ofs("test_bfgs_routine_case2.log"); testing::internal::CaptureStdout(); - bfgs.bfgs_routine_for_testing(lat0, istep, update_iter, ofs, etot_info, criteria.out_level, criteria.test_relax_method); + bfgs.bfgs_routine_for_testing(lat0, + istep, + update_iter, + ofs, + etot_info, + criteria.out_level, + criteria.test_relax_method); std::string std_outout = testing::internal::GetCapturedStdout(); ofs.close(); @@ -328,7 +345,7 @@ TEST_F(IonsMoveBFGSTest, BfgsRoutineCase2) std::string expected_ofs = " quadratic interpolation is impossible.\n"; std::string expected_std = ""; - EXPECT_THAT(ofs_output, ::testing::HasSubstr(expected_ofs)); + EXPECT_THAT(ofs_output, ::testing::HasSubstr(expected_ofs)); EXPECT_EQ(expected_std, std_outout); EXPECT_DOUBLE_EQ(Ions_Move_Basic::trust_radius, -0.5); @@ -374,7 +391,13 @@ TEST_F(IonsMoveBFGSTest, BfgsRoutineCase3) // Call the function being tested std::ofstream ofs("test_bfgs_routine_case3.log"); - bfgs.bfgs_routine_for_testing(lat0, istep, update_iter, ofs, etot_info, criteria.out_level, criteria.test_relax_method); + bfgs.bfgs_routine_for_testing(lat0, + istep, + update_iter, + ofs, + etot_info, + criteria.out_level, + criteria.test_relax_method); ofs.close(); // Check the results @@ -385,7 +408,7 @@ TEST_F(IonsMoveBFGSTest, BfgsRoutineCase3) std::string expected_ofs = " check the norm of new move 410 (Bohr)\n Uphill move : resetting bfgs history\n"; - EXPECT_THAT(ofs_output, ::testing::HasSubstr(expected_ofs)); + EXPECT_THAT(ofs_output, ::testing::HasSubstr(expected_ofs)); EXPECT_DOUBLE_EQ(Ions_Move_Basic::trust_radius, 0.2); EXPECT_DOUBLE_EQ(etot_info[0], 0.9); EXPECT_DOUBLE_EQ(bfgs.get_tr_min_hit(), false); @@ -432,7 +455,15 @@ TEST_F(IonsMoveBFGSTest, BfgsRoutineWarningQuit1) // Check the results std::ofstream ofs("test_bfgs_routine_warning_quit1.log"); testing::internal::CaptureStdout(); - EXPECT_EXIT(bfgs.bfgs_routine_for_testing(lat0, istep, update_iter, ofs, etot_info, criteria.out_level, criteria.test_relax_method), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(bfgs.bfgs_routine_for_testing(lat0, + istep, + update_iter, + ofs, + etot_info, + criteria.out_level, + criteria.test_relax_method), + ::testing::ExitedWithCode(1), + ""); std::string output = testing::internal::GetCapturedStdout(); ofs.close(); std::remove("test_bfgs_routine_warning_quit1.log"); @@ -456,7 +487,15 @@ TEST_F(IonsMoveBFGSTest, BfgsRoutineWarningQuit2) // Check the results std::ofstream ofs("test_bfgs_routine_warning_quit2.log"); testing::internal::CaptureStdout(); - EXPECT_EXIT(bfgs.bfgs_routine_for_testing(lat0, istep, update_iter, ofs, etot_info, criteria.out_level, criteria.test_relax_method), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(bfgs.bfgs_routine_for_testing(lat0, + istep, + update_iter, + ofs, + etot_info, + criteria.out_level, + criteria.test_relax_method), + ::testing::ExitedWithCode(1), + ""); std::string output = testing::internal::GetCapturedStdout(); ofs.close(); std::remove("test_bfgs_routine_warning_quit2.log"); diff --git a/source/source_relax/unittests/test_ions_move_bfgs2.cpp b/source/source_relax/unittests/test_ions_move_bfgs2.cpp index 8d176b780f7..c4f0022b37e 100644 --- a/source/source_relax/unittests/test_ions_move_bfgs2.cpp +++ b/source/source_relax/unittests/test_ions_move_bfgs2.cpp @@ -1,22 +1,21 @@ // source/source_relax/unittests/test_ions_move_bfgs2.cpp -#include "gtest/gtest.h" -#include "gmock/gmock.h" #include "for_test.h" - -#include "source_relax/ions_move_bfgs2.h" - #include "source_io/module_parameter/parameter.h" - #include "source_relax/ions_move_basic.h" // for Ions_Move_Basic static members +#include "source_relax/ions_move_bfgs2.h" #include "source_relax/relax_data.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" #include /************************************************ * unit tests for Ions_Move_BFGS2 (no MockUnitCell) ***********************************************/ -class BFGSTest : public ::testing::Test { -protected: +class BFGSTest : public ::testing::Test +{ + protected: Ions_Move_BFGS2 bfgs; void SetUp() override { @@ -167,13 +166,21 @@ TEST_F(BFGSTest, RelaxStepAutoInitialize) ucell.atoms[0].tau = std::vector>(2); ucell.atoms[0].taud = std::vector>(2); ucell.atoms[0].mbl = std::vector>(2, {1, 1, 1}); - ucell.atoms[0].tau[0].x = 0.0; ucell.atoms[0].tau[0].y = 0.0; ucell.atoms[0].tau[0].z = 0.0; - ucell.atoms[0].tau[1].x = 1.0; ucell.atoms[0].tau[1].y = 0.0; ucell.atoms[0].tau[1].z = 0.0; + ucell.atoms[0].tau[0].x = 0.0; + ucell.atoms[0].tau[0].y = 0.0; + ucell.atoms[0].tau[0].z = 0.0; + ucell.atoms[0].tau[1].x = 1.0; + ucell.atoms[0].tau[1].y = 0.0; + ucell.atoms[0].tau[1].z = 0.0; ucell.lat0 = 1.0; ModuleBase::matrix force(2, 3); - force(0, 0) = 0.1; force(0, 1) = 0.0; force(0, 2) = 0.0; - force(1, 0) = -0.1; force(1, 1) = 0.0; force(1, 2) = 0.0; + force(0, 0) = 0.1; + force(0, 1) = 0.0; + force(0, 2) = 0.0; + force(1, 0) = -0.1; + force(1, 1) = 0.0; + force(1, 2) = 0.0; std::ofstream ofs_running; @@ -198,10 +205,8 @@ TEST_F(BFGSTest, DetermineStepScaling) bfgs.allocate(size); std::vector steplength = {1.0, 0.1}; - std::vector> dpos = { - ModuleBase::Vector3(1.0, 1.0, 1.0), - ModuleBase::Vector3(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; call_DetermineStep(steplength, dpos, maxstep); @@ -262,15 +267,21 @@ TEST_F(BFGSTest, GetPosAndPostaud) ucell.atoms[0].mbl = std::vector>(2, {1, 1, 1}); // set coordinates - ucell.atoms[0].tau[0].x = 0.0; ucell.atoms[0].tau[0].y = 0.0; ucell.atoms[0].tau[0].z = 0.0; - ucell.atoms[0].tau[1].x = 1.0; ucell.atoms[0].tau[1].y = 0.0; ucell.atoms[0].tau[1].z = 0.0; - - // allocate mapping arrays + ucell.atoms[0].tau[0].x = 0.0; + ucell.atoms[0].tau[0].y = 0.0; + ucell.atoms[0].tau[0].z = 0.0; + ucell.atoms[0].tau[1].x = 1.0; + ucell.atoms[0].tau[1].y = 0.0; + ucell.atoms[0].tau[1].z = 0.0; + + // allocate mapping arrays ucell.iat2it.resize(ucell.nat); ucell.iat2ia.resize(ucell.nat); int k = 0; - for (int it = 0; it < ucell.ntype; ++it) { - for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { ucell.iat2it[k] = it; ucell.iat2ia[k] = ia; ++k; @@ -304,8 +315,10 @@ TEST_F(BFGSTest, CalculateLargestGrad) ucell.iat2it.resize(ucell.nat); ucell.iat2ia.resize(ucell.nat); int k = 0; - for (int it = 0; it < ucell.ntype; ++it) { - for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { ucell.iat2it[k] = it; ucell.iat2ia[k] = ia; ++k; @@ -314,10 +327,10 @@ TEST_F(BFGSTest, CalculateLargestGrad) // build force matrix: 2 atoms x 3 components ModuleBase::matrix force(2, 3); - force(0, 0) = -2.0; // this yields grad component = -(-2.0)*lat0 = 4.0 -> divided by lat0 => 2.0 + force(0, 0) = -2.0; // this yields grad component = -(-2.0)*lat0 = 4.0 -> divided by lat0 => 2.0 force(0, 1) = 0.0; force(0, 2) = 1.0; - force(1, 0) = 3.0; // this yields abs = 6.0 -> divided by lat0 => 3.0 (this should be largest) + force(1, 0) = 3.0; // this yields abs = 6.0 -> divided by lat0 => 3.0 (this should be largest) force(1, 1) = -1.0; force(1, 2) = 0.0; @@ -346,20 +359,30 @@ TEST_F(BFGSTest, RelaxStepBasic) ucell.iat2it.resize(ucell.nat); ucell.iat2ia.resize(ucell.nat); int k = 0; - for (int it = 0; it < ucell.ntype; ++it) { - for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { ucell.iat2it[k] = it; ucell.iat2ia[k] = ia; ++k; } } // Set initial positions - ucell.atoms[0].tau[0].x = 0.0; ucell.atoms[0].tau[0].y = 0.0; ucell.atoms[0].tau[0].z = 0.0; - ucell.atoms[0].tau[1].x = 1.0; ucell.atoms[0].tau[1].y = 0.0; ucell.atoms[0].tau[1].z = 0.0; + ucell.atoms[0].tau[0].x = 0.0; + ucell.atoms[0].tau[0].y = 0.0; + ucell.atoms[0].tau[0].z = 0.0; + ucell.atoms[0].tau[1].x = 1.0; + ucell.atoms[0].tau[1].y = 0.0; + ucell.atoms[0].tau[1].z = 0.0; // Setup force matrix ModuleBase::matrix force(2, 3); - force(0, 0) = 0.1; force(0, 1) = 0.0; force(0, 2) = 0.0; - force(1, 0) = -0.1; force(1, 1) = 0.0; force(1, 2) = 0.0; + force(0, 0) = 0.1; + force(0, 1) = 0.0; + force(0, 2) = 0.0; + force(1, 0) = -0.1; + force(1, 1) = 0.0; + force(1, 2) = 0.0; // Allocate and call relax_step bfgs.allocate(ucell.nat); diff --git a/source/source_relax/unittests/test_ions_move_cg.cpp b/source/source_relax/unittests/test_ions_move_cg.cpp index 5509bc776ed..91a11d63afa 100644 --- a/source/source_relax/unittests/test_ions_move_cg.cpp +++ b/source/source_relax/unittests/test_ions_move_cg.cpp @@ -1,10 +1,11 @@ -#include "source_relax/relax_criteria.h" -#include #include "for_test.h" -#include "gtest/gtest.h" -#include "gmock/gmock.h" #include "source_relax/ions_move_basic.h" #include "source_relax/ions_move_cg.h" +#include "source_relax/relax_criteria.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" +#include /************************************************ * unit tests of class Ions_Move_CG @@ -24,7 +25,7 @@ class IonsMoveCGTest : public ::testing::Test im_cg.allocate(Ions_Move_Basic::dim); criteria.force_thr = 0.001; - // ban the 'cout' + // ban the 'cout' // mohan add 2025-05-02 std::cout.rdbuf(NULL); } @@ -38,10 +39,10 @@ class IonsMoveCGTest : public ::testing::Test for (int it = 0; it < ucell.ntype; it++) { Atom* atom = &ucell.atoms[it]; - atom->label="test"; + atom->label = "test"; for (int ia = 0; ia < atom->na; ia++) { - atom->mag[ia]= 1; + atom->mag[ia] = 1; for (int ik = 0; ik < 3; ++ik) { atom->tau[ia][ik] = 1; @@ -142,7 +143,6 @@ TEST_F(IonsMoveCGTest, TestStartConverged) ifs.close(); std::remove("TestStartConverged.log"); - std::regex pattern(R"(==> .*::.*\t[\d\.]+ GB\t\d+ s\n )"); output = std::regex_replace(output, pattern, ""); EXPECT_THAT(output, testing::HasSubstr(expected_output)); diff --git a/source/source_relax/unittests/test_ions_move_lbfgs.cpp b/source/source_relax/unittests/test_ions_move_lbfgs.cpp index 93275d92b9e..416f2d18a2f 100644 --- a/source/source_relax/unittests/test_ions_move_lbfgs.cpp +++ b/source/source_relax/unittests/test_ions_move_lbfgs.cpp @@ -1,14 +1,14 @@ -#include -#include -#include - -#include "gmock/gmock.h" -#include "gtest/gtest.h" #include "source_cell/unitcell.h" #include "source_relax/ions_move_basic.h" #include "source_relax/ions_move_lbfgs.h" #include "source_relax/relax_criteria.h" +#include "gmock/gmock.h" +#include "gtest/gtest.h" +#include +#include +#include + /************************************************ * unit tests of class Ions_Move_LBFGS ***********************************************/ diff --git a/source/source_relax/unittests/test_ions_move_methods.cpp b/source/source_relax/unittests/test_ions_move_methods.cpp index 794d49c2073..d9b25d9971c 100644 --- a/source/source_relax/unittests/test_ions_move_methods.cpp +++ b/source/source_relax/unittests/test_ions_move_methods.cpp @@ -1,5 +1,7 @@ -#include "source_relax/relax_criteria.h" #include "for_test.h" +#include "source_relax/ions_move_methods.h" +#include "source_relax/relax_criteria.h" + #include "gmock/gmock.h" #include "gtest/gtest.h" #include @@ -8,7 +10,6 @@ #include #include #include -#include "source_relax/ions_move_methods.h" /************************************************ * unit tests of class Ions_Move_Methods ***********************************************/ @@ -24,68 +25,67 @@ * - Ions_Move_Methods::get_update_iter() */ - // Mock the remake_cell function from update_cell.h +// Mock the remake_cell function from update_cell.h namespace unitcell { - // Track if remake_cell was called and with what lattice - static bool remake_cell_called = false; - static std::string remake_cell_latName; - static ModuleBase::Matrix3 remake_cell_latvec; - - void remake_cell(Lattice& lat) - { - remake_cell_called = true; - remake_cell_latName = lat.latName; - remake_cell_latvec = lat.latvec; +// Track if remake_cell was called and with what lattice +static bool remake_cell_called = false; +static std::string remake_cell_latName; +static ModuleBase::Matrix3 remake_cell_latvec; - // Mock implementation: enforce simple cubic structure for "sc" - if (lat.latName == "sc") - { - double celldm = std::sqrt(lat.latvec.e11 * lat.latvec.e11 + - lat.latvec.e12 * lat.latvec.e12 + - lat.latvec.e13 * lat.latvec.e13); - lat.latvec.Zero(); - lat.latvec.e11 = celldm; - lat.latvec.e22 = celldm; - lat.latvec.e33 = celldm; - } - // Mock implementation: enforce FCC structure for "fcc" - else if (lat.latName == "fcc") - { - double celldm = std::sqrt(lat.latvec.e11 * lat.latvec.e11 + - lat.latvec.e12 * lat.latvec.e12 + - lat.latvec.e13 * lat.latvec.e13) / std::sqrt(2.0); - lat.latvec.e11 = -celldm; - lat.latvec.e12 = 0.0; - lat.latvec.e13 = celldm; - lat.latvec.e21 = 0.0; - lat.latvec.e22 = celldm; - lat.latvec.e23 = celldm; - lat.latvec.e31 = -celldm; - lat.latvec.e32 = celldm; - lat.latvec.e33 = 0.0; - } - } +void remake_cell(Lattice& lat) +{ + remake_cell_called = true; + remake_cell_latName = lat.latName; + remake_cell_latvec = lat.latvec; - void update_pos_tau(const Lattice&, const double*, const int, const int, Atom*) + // Mock implementation: enforce simple cubic structure for "sc" + if (lat.latName == "sc") { + double celldm = std::sqrt(lat.latvec.e11 * lat.latvec.e11 + lat.latvec.e12 * lat.latvec.e12 + + lat.latvec.e13 * lat.latvec.e13); + lat.latvec.Zero(); + lat.latvec.e11 = celldm; + lat.latvec.e22 = celldm; + lat.latvec.e33 = celldm; } - - // Helper function to reset mock state - void reset_remake_cell_mock() + // Mock implementation: enforce FCC structure for "fcc" + else if (lat.latName == "fcc") { - remake_cell_called = false; - remake_cell_latName = ""; - remake_cell_latvec.Zero(); + double celldm = std::sqrt(lat.latvec.e11 * lat.latvec.e11 + lat.latvec.e12 * lat.latvec.e12 + + lat.latvec.e13 * lat.latvec.e13) + / std::sqrt(2.0); + lat.latvec.e11 = -celldm; + lat.latvec.e12 = 0.0; + lat.latvec.e13 = celldm; + lat.latvec.e21 = 0.0; + lat.latvec.e22 = celldm; + lat.latvec.e23 = celldm; + lat.latvec.e31 = -celldm; + lat.latvec.e32 = celldm; + lat.latvec.e33 = 0.0; } +} - // Helper function to check if remake_cell was called - bool was_remake_cell_called() - { - return remake_cell_called; - } +void update_pos_tau(const Lattice&, const double*, const int, const int, Atom*) +{ } +// Helper function to reset mock state +void reset_remake_cell_mock() +{ + remake_cell_called = false; + remake_cell_latName = ""; + remake_cell_latvec.Zero(); +} + +// Helper function to check if remake_cell was called +bool was_remake_cell_called() +{ + return remake_cell_called; +} +} // namespace unitcell + // Define a fixture for the tests class IonsMoveMethodsTest : public ::testing::Test { @@ -191,8 +191,6 @@ TEST_F(IonsMoveMethodsTest, CalMovementWarningQuit) std::remove("log"); } - - // Test the get_ediff() function TEST_F(IonsMoveMethodsTest, GetEdiff) { diff --git a/source/source_relax/unittests/test_ions_move_sd.cpp b/source/source_relax/unittests/test_ions_move_sd.cpp index 507b3b216a9..59c93f8585f 100644 --- a/source/source_relax/unittests/test_ions_move_sd.cpp +++ b/source/source_relax/unittests/test_ions_move_sd.cpp @@ -1,10 +1,11 @@ -#include "source_relax/relax_criteria.h" -#include #include "for_test.h" -#include "gmock/gmock.h" -#include "gtest/gtest.h" #include "source_relax/ions_move_basic.h" #include "source_relax/ions_move_sd.h" +#include "source_relax/relax_criteria.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" +#include /************************************************ * unit tests of class Ions_Move_SD @@ -121,7 +122,7 @@ TEST_F(IonsMoveSDTest, TestStartNotConverged) { for (int ik = 0; ik < 3; ++ik) { - atom->tau[ia][ik] = (ik + 1)/3; + atom->tau[ia][ik] = (ik + 1) / 3; atom->mbl[ia][ik] = 1; } } diff --git a/source/source_relax/unittests/test_lat_change_method.cpp b/source/source_relax/unittests/test_lat_change_method.cpp index e0b8304bafb..4cd743152e7 100644 --- a/source/source_relax/unittests/test_lat_change_method.cpp +++ b/source/source_relax/unittests/test_lat_change_method.cpp @@ -1,9 +1,9 @@ -#include -#include "source_relax/lat_change_method.h" +#include "for_test.h" #include "mock_remake_cell.h" +#include "source_relax/lat_change_method.h" -#include "for_test.h" #include "gtest/gtest.h" +#include /************************************************ * unit tests of class Lattice_Change_Methods ***********************************************/ @@ -25,7 +25,12 @@ void Lattice_Change_CG::allocate(void) { } -bool Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_in, const double &etot_in, std::ofstream& ofs, std::vector& etot_info, const Relax_Criteria& criteria) +bool Lattice_Change_CG::start(UnitCell& ucell, + const ModuleBase::matrix& stress_in, + const double& etot_in, + std::ofstream& ofs, + std::vector& etot_info, + const Relax_Criteria& criteria) { return false; } diff --git a/source/source_relax/unittests/test_lattice_change_basic.cpp b/source/source_relax/unittests/test_lattice_change_basic.cpp index a77254da56d..4013d03c5d7 100644 --- a/source/source_relax/unittests/test_lattice_change_basic.cpp +++ b/source/source_relax/unittests/test_lattice_change_basic.cpp @@ -1,12 +1,12 @@ -#include +#include "for_test.h" +#include "mock_remake_cell.h" +#include "source_io/module_parameter/parameter.h" #include "source_relax/lattice_change_basic.h" #include "source_relax/relax_data.h" -#include "mock_remake_cell.h" -#include "for_test.h" -#include "gtest/gtest.h" #include "gmock/gmock.h" -#include "source_io/module_parameter/parameter.h" +#include "gtest/gtest.h" +#include /************************************************ * unit tests of namespace Lattice_Change_Basic @@ -219,7 +219,8 @@ TEST_F(LatticeChangeBasicTest, ChangeLattice) EXPECT_DOUBLE_EQ(ucell.invGGT.e33, 0.36); } -// Test for check_converged when ucell.lat_axis_free[0] == 1 && ucell.lat_axis_free[1] == 1 && ucell.lat_axis_free[2] == 1, but not converged +// Test for check_converged when ucell.lat_axis_free[0] == 1 && ucell.lat_axis_free[1] == 1 && ucell.lat_axis_free[2] == +// 1, but not converged TEST_F(LatticeChangeBasicTest, CheckConvergedCase1) { // Set up test data @@ -256,7 +257,8 @@ TEST_F(LatticeChangeBasicTest, CheckConvergedCase1) std::remove("test_check_converged_case1.log"); } -// Test for check_converged when ucell.lat_axis_free[0] == 1 && ucell.lat_axis_free[1] == 1 && ucell.lat_axis_free[2] == 1 && largest_grad == 0 +// Test for check_converged when ucell.lat_axis_free[0] == 1 && ucell.lat_axis_free[1] == 1 && ucell.lat_axis_free[2] == +// 1 && largest_grad == 0 TEST_F(LatticeChangeBasicTest, CheckConvergedCase2) { // Set up test data @@ -293,7 +295,8 @@ TEST_F(LatticeChangeBasicTest, CheckConvergedCase2) std::remove("test_check_converged_case2.log"); } -// Test for check_converged when ucell.lat_axis_free[0] == 1 && ucell.lat_axis_free[1] == 1 && ucell.lat_axis_free[2] == 1, and converged +// Test for check_converged when ucell.lat_axis_free[0] == 1 && ucell.lat_axis_free[1] == 1 && ucell.lat_axis_free[2] == +// 1, and converged TEST_F(LatticeChangeBasicTest, CheckConvergedCase3) { // Set up test data @@ -478,7 +481,6 @@ TEST_F(LatticeChangeBasicTest, TerminateNotConverged) std::remove("test_terminate_not_converged.log"); } - TEST_F(LatticeChangeBasicTest, SetupEtot) { Lattice_Change_Basic::stress_step = 2; @@ -587,7 +589,7 @@ TEST_F(LatticeChangeBasicTest, ChangeLatticeVolumeRescaling) EXPECT_NEAR(ucell.omega, 1000.0, 1e-8); // Check that lattice vectors were rescaled uniformly - double expected_scale = std::pow(1000.0 / 1331.0, 1.0/3.0); // (old_vol / new_vol)^(1/3) + double expected_scale = std::pow(1000.0 / 1331.0, 1.0 / 3.0); // (old_vol / new_vol)^(1/3) EXPECT_NEAR(ucell.latvec.e11, 1.1 * expected_scale, 1e-10); EXPECT_NEAR(ucell.latvec.e22, 1.1 * expected_scale, 1e-10); EXPECT_NEAR(ucell.latvec.e33, 1.1 * expected_scale, 1e-10); @@ -804,7 +806,8 @@ TEST_F(LatticeChangeBasicTest, ChangeLatticeFixedIbravFCC) lat[8] = 0.0; // Apply a small move - for (int i = 0; i < 9; i++) move[i] = 0.01 * ucell.lat0; + for (int i = 0; i < 9; i++) + move[i] = 0.01 * ucell.lat0; fixed_ibrav = true; Lattice_Change_Basic::fixed_axes = "None"; @@ -821,15 +824,12 @@ TEST_F(LatticeChangeBasicTest, ChangeLatticeFixedIbravFCC) // Check that lattice maintains FCC structure // For FCC: a1 = (-a, 0, a), a2 = (0, a, a), a3 = (-a, a, 0) // All should have same magnitude - double mag1 = std::sqrt(ucell.latvec.e11*ucell.latvec.e11 + - ucell.latvec.e12*ucell.latvec.e12 + - ucell.latvec.e13*ucell.latvec.e13); - double mag2 = std::sqrt(ucell.latvec.e21*ucell.latvec.e21 + - ucell.latvec.e22*ucell.latvec.e22 + - ucell.latvec.e23*ucell.latvec.e23); - double mag3 = std::sqrt(ucell.latvec.e31*ucell.latvec.e31 + - ucell.latvec.e32*ucell.latvec.e32 + - ucell.latvec.e33*ucell.latvec.e33); + double mag1 = std::sqrt(ucell.latvec.e11 * ucell.latvec.e11 + ucell.latvec.e12 * ucell.latvec.e12 + + ucell.latvec.e13 * ucell.latvec.e13); + double mag2 = std::sqrt(ucell.latvec.e21 * ucell.latvec.e21 + ucell.latvec.e22 * ucell.latvec.e22 + + ucell.latvec.e23 * ucell.latvec.e23); + double mag3 = std::sqrt(ucell.latvec.e31 * ucell.latvec.e31 + ucell.latvec.e32 * ucell.latvec.e32 + + ucell.latvec.e33 * ucell.latvec.e33); EXPECT_NEAR(mag1, mag2, 1e-10); EXPECT_NEAR(mag2, mag3, 1e-10); diff --git a/source/source_relax/unittests/test_lattice_change_cg.cpp b/source/source_relax/unittests/test_lattice_change_cg.cpp index fc073052cb9..b4070009e2c 100644 --- a/source/source_relax/unittests/test_lattice_change_cg.cpp +++ b/source/source_relax/unittests/test_lattice_change_cg.cpp @@ -1,9 +1,10 @@ -#include "source_relax/relax_criteria.h" #include "for_test.h" -#include "gtest/gtest.h" #include "mock_remake_cell.h" #include "source_relax/lattice_change_basic.h" #include "source_relax/lattice_change_cg.h" +#include "source_relax/relax_criteria.h" + +#include "gtest/gtest.h" /************************************************ * unit tests of class Lattice_Change_CG @@ -108,9 +109,9 @@ TEST_F(LatticeChangeCGTest, TestStartConverged) ofs.close(); // Check output - std::string expected_output - = " Largest stress is 0, movement is impossible.\n end of lattice optimization\n stress_step = 1\n " - " update iteration = 5\n"; + std::string expected_output = " Largest stress is 0, movement is impossible.\n end of lattice optimization\n " + " stress_step = 1\n " + " update iteration = 5\n"; std::ifstream ifs("test_lc_cg_start_converged.log"); std::string output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); @@ -336,4 +337,3 @@ TEST_F(LatticeChangeCGTest, TestStartNoTrial) ifs.close(); std::remove("test_lc_cg_start_notrial.log"); } - diff --git a/source/source_relax/unittests/test_line_search.cpp b/source/source_relax/unittests/test_line_search.cpp index 97382e64193..ec01a354ab3 100644 --- a/source/source_relax/unittests/test_line_search.cpp +++ b/source/source_relax/unittests/test_line_search.cpp @@ -1,10 +1,9 @@ #include "../line_search.h" -#include -#include - #include "gmock/gmock.h" #include "gtest/gtest.h" +#include +#include /************************************************ * unit test of class Line_Search ***********************************************/ diff --git a/source/source_relax/unittests/test_matrix_methods.cpp b/source/source_relax/unittests/test_matrix_methods.cpp index c40d86b32ef..72e754a365e 100644 --- a/source/source_relax/unittests/test_matrix_methods.cpp +++ b/source/source_relax/unittests/test_matrix_methods.cpp @@ -1,7 +1,6 @@ #include "../matrix_methods.h" #include - #include // Unit tests for the free matrix/vector helpers in matrix_methods.cpp. diff --git a/source/source_relax/unittests/test_relax_history.cpp b/source/source_relax/unittests/test_relax_history.cpp index c00b13ffaae..3496625bb20 100644 --- a/source/source_relax/unittests/test_relax_history.cpp +++ b/source/source_relax/unittests/test_relax_history.cpp @@ -33,9 +33,9 @@ TEST(RelaxHistory, LineBreakEveryFive) TEST(RelaxHistory, LongHistoryTruncated) { std::vector hist(150, 0.0); - hist.front() = 1.0; // first value kept - hist.back() = 2.0; // last value kept - hist[75] = 9.0; // middle value must be omitted + hist.front() = 1.0; // first value kept + hist.back() = 2.0; // last value kept + hist[75] = 9.0; // middle value must be omitted const std::string out = format_relax_history(hist); EXPECT_THAT(out, HasSubstr("... (omitted 130 step(s)) ...")); EXPECT_THAT(out, HasSubstr("1.000e+00")); diff --git a/source/source_relax/unittests/test_relax_nsync.cpp b/source/source_relax/unittests/test_relax_nsync.cpp index 7939822354e..7351c5eccf3 100644 --- a/source/source_relax/unittests/test_relax_nsync.cpp +++ b/source/source_relax/unittests/test_relax_nsync.cpp @@ -1,16 +1,16 @@ -#include -#include -#include -#include - -#include "gmock/gmock.h" -#include "gtest/gtest.h" #include "source_base/global_variable.h" #include "source_cell/unitcell.h" #include "source_io/module_parameter/parameter.h" #include "source_relax/lattice_change_basic.h" #include "source_relax/relax_nsync.h" +#include "gmock/gmock.h" +#include "gtest/gtest.h" +#include +#include +#include +#include + /************************************************ * unit tests of class IonCellOptimizer (relax_nsync.cpp) ***********************************************/ @@ -154,9 +154,9 @@ class IonCellOptimizerTest : public ::testing::Test inp.calculation = "relax"; inp.relax_method = "lbfgs"; inp.relax_nmax = 50; - inp.force_thr = 0.1; // Ry/Bohr - inp.force_thr_ev = 0.1 * 13.6058 / 0.529177; // eV/Angstrom - inp.stress_thr = 0.5; // kbar + inp.force_thr = 0.1; // Ry/Bohr + inp.force_thr_ev = 0.1 * 13.6058 / 0.529177; // eV/Angstrom + inp.stress_thr = 0.5; // kbar inp.cal_force = 1; inp.cal_stress = 1; inp.out_level = "m"; diff --git a/source/source_relax/unittests/test_relax_stru_io.cpp b/source/source_relax/unittests/test_relax_stru_io.cpp index 1d339a23c4c..1c553cadfe8 100644 --- a/source/source_relax/unittests/test_relax_stru_io.cpp +++ b/source/source_relax/unittests/test_relax_stru_io.cpp @@ -2,7 +2,6 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" - #include #include #include @@ -96,7 +95,7 @@ TEST(RelaxStruIO, HeaderWritesThreeStressRowsInKbar) // driver-loop control flow is covered by the relax sync/nsync suites. class HeaderBehaviorTest : public testing::Test { -protected: + protected: Input_para inp; ModuleBase::matrix stress; const double etot = -10.5; // Ry @@ -184,7 +183,9 @@ TEST_F(HeaderBehaviorTest, CalForceFalse) ASSERT_EQ(lines.size(), 7u); EXPECT_THAT(header, testing::Not(testing::HasSubstr("N/A"))); - EXPECT_THAT(lines[6], testing::HasSubstr("# NOTE: forces not computed (cal_force=0); per-atom f fields omitted intentionally")); + EXPECT_THAT( + lines[6], + testing::HasSubstr("# NOTE: forces not computed (cal_force=0); per-atom f fields omitted intentionally")); } TEST_F(HeaderBehaviorTest, CalForceStressBothFalse) @@ -197,7 +198,8 @@ TEST_F(HeaderBehaviorTest, CalForceStressBothFalse) ASSERT_EQ(lines.size(), 7u); EXPECT_EQ(count_substr(header, "# Stress (kbar): N/A N/A N/A"), 3); - EXPECT_THAT(lines[6], testing::HasSubstr("# NOTE: stress not computed (cal_stress=0); forces not computed (cal_force=0)")); + EXPECT_THAT(lines[6], + testing::HasSubstr("# NOTE: stress not computed (cal_stress=0); forces not computed (cal_force=0)")); } TEST_F(HeaderBehaviorTest, FinalStep) diff --git a/source/source_relax/unittests/test_relax_sync.cpp b/source/source_relax/unittests/test_relax_sync.cpp index e57b6e1fc67..dc27dd0157f 100644 --- a/source/source_relax/unittests/test_relax_sync.cpp +++ b/source/source_relax/unittests/test_relax_sync.cpp @@ -1,50 +1,38 @@ -#include "gtest/gtest.h" -#include "gmock/gmock.h" -#include #include "../relax_sync.h" +#include "relax_test.h" #include "source_cell/unitcell.h" #include "source_io/module_parameter/input_parameter.h" -#include "relax_test.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" #include +#include // Fixed mock SCF outputs (forces, stresses, energies) feeding relax_step, // previously read from support/*.txt. Inlining them keeps the test // independent of the working directory. Values are the first 3 steps of the // original files, so result_ref is unchanged. -namespace { +namespace +{ constexpr int kRelaxSteps = 3; constexpr int kRelaxNat = 5; -constexpr int kRelaxNforce = kRelaxSteps * kRelaxNat * 3; // 45 -constexpr int kRelaxNstress = kRelaxSteps * 3 * 3; // 27 +constexpr int kRelaxNforce = kRelaxSteps * kRelaxNat * 3; // 45 +constexpr int kRelaxNstress = kRelaxSteps * 3 * 3; // 27 const double kRelaxForces[kRelaxNforce] = { - 4.56e-08, -8.74e-08, 4.0671950e-03, - -5.59e-08, 1.479e-07, 1.81171376e-02, - 2.942e-07, -5.56e-08, -1.7979905e-03, - -1.691e-07, 1.612e-07, -1.7977699e-03, - -1.149e-07, -1.660e-07, -1.85885722e-02, - 6.73e-08, 5.90e-08, 4.8608620e-03, - -1.632e-07, -1.833e-07, -2.24081615e-02, - 4.93e-08, 1.246e-07, -1.8938742e-03, - 9.84e-08, 4.34e-08, -1.8939930e-03, - -5.19e-08, -4.38e-08, 2.13351667e-02, - 8.25e-08, 2.095e-07, 4.5297162e-03, - -5.19e-08, 3.57e-08, -3.6659292e-03, - 4.95e-08, -1.654e-07, -1.8863808e-03, - 4.14e-08, 2.440e-07, -1.8870376e-03, - -1.215e-07, -3.239e-07, 2.9096315e-03, + 4.56e-08, -8.74e-08, 4.0671950e-03, -5.59e-08, 1.479e-07, 1.81171376e-02, 2.942e-07, + -5.56e-08, -1.7979905e-03, -1.691e-07, 1.612e-07, -1.7977699e-03, -1.149e-07, -1.660e-07, + -1.85885722e-02, 6.73e-08, 5.90e-08, 4.8608620e-03, -1.632e-07, -1.833e-07, -2.24081615e-02, + 4.93e-08, 1.246e-07, -1.8938742e-03, 9.84e-08, 4.34e-08, -1.8939930e-03, -5.19e-08, + -4.38e-08, 2.13351667e-02, 8.25e-08, 2.095e-07, 4.5297162e-03, -5.19e-08, 3.57e-08, + -3.6659292e-03, 4.95e-08, -1.654e-07, -1.8863808e-03, 4.14e-08, 2.440e-07, -1.8870376e-03, + -1.215e-07, -3.239e-07, 2.9096315e-03, }; const double kRelaxStresses[kRelaxNstress] = { - 4.79987e-05, -2.0e-10, 6.0e-10, - -2.0e-10, 4.80038e-05, -4.0e-10, - 6.0e-10, -4.0e-10, -2.751687e-04, - 1.040725e-04, 2.0e-10, -1.6e-09, - 2.0e-10, 1.040682e-04, -7.0e-10, - -1.6e-09, -7.0e-10, 1.982021e-04, - 8.00586e-05, 1.0e-10, -1.0e-09, - 1.0e-10, 8.00526e-05, 2.8e-09, - -1.0e-09, 2.8e-09, -6.4010e-06, + 4.79987e-05, -2.0e-10, 6.0e-10, -2.0e-10, 4.80038e-05, -4.0e-10, 6.0e-10, -4.0e-10, -2.751687e-04, + 1.040725e-04, 2.0e-10, -1.6e-09, 2.0e-10, 1.040682e-04, -7.0e-10, -1.6e-09, -7.0e-10, 1.982021e-04, + 8.00586e-05, 1.0e-10, -1.0e-09, 1.0e-10, 8.00526e-05, 2.8e-09, -1.0e-09, 2.8e-09, -6.4010e-06, }; const double kRelaxEnergies[kRelaxSteps] = { @@ -56,342 +44,431 @@ const double kRelaxEnergies[kRelaxSteps] = { class Test_SETGRAD : public testing::Test { - protected: - Relax rl; - std::vector result; - Input_para inp; - UnitCell ucell; - std::ofstream ofs; - - void SetUp() - { - inp.force_thr = 0.001; - inp.calculation = "cell-relax"; - - ModuleBase::matrix force_in, stress_in; - int nat = 3; - - force_in.create(nat,3); - stress_in.create(3,3); - - force_in(0,0) = 0.1; force_in(0,1) = 0.1; force_in(0,2)= 0.1; - force_in(1,0) = 0; force_in(1,1) = 0.1; force_in(1,2)= 0.1; - force_in(2,0) = 0; force_in(2,1) = 0; force_in(2,2)= 0.1; - - stress_in(0,0) = 1; stress_in(0,1) = 1; stress_in(0,2)= 1; - stress_in(1,0) = 0; stress_in(1,1) = 1; stress_in(1,2)= 1; - stress_in(2,0) = 0; stress_in(2,1) = 0; stress_in(2,2)= 1; - - ucell.ntype = 1; - ucell.nat = nat; - ucell.atoms = new Atom[1]; - ucell.atoms[0].na = nat; - ucell.omega = 1.0; - ucell.lat0 = 1.0; - - ucell.iat2it.resize(nat); - ucell.iat2ia.resize(nat); - ucell.atoms[0].mbl.resize(nat); - ucell.atoms[0].taud.resize(nat); - ucell.atoms[0].tau.resize(nat); - ucell.atoms[0].dis.resize(nat); - ucell.atoms[0].mag.resize(nat); - ucell.atoms[0].vel.resize(nat); - - ucell.iat2it[0] = 0; - ucell.iat2it[1] = 0; - ucell.iat2it[2] = 0; - - ucell.iat2ia[0] = 0; - ucell.iat2ia[1] = 1; - ucell.iat2ia[2] = 2; - - ucell.atoms[0].mbl[0].x = 0; - ucell.atoms[0].mbl[0].y = 0; - ucell.atoms[0].mbl[0].z = 1; - - ucell.atoms[0].mbl[1].x = 0; - ucell.atoms[0].mbl[1].y = 1; - ucell.atoms[0].mbl[1].z = 0; - - ucell.atoms[0].mbl[2].x = 1; - ucell.atoms[0].mbl[2].y = 0; - ucell.atoms[0].mbl[2].z = 0; - - ucell.atoms[0].taud[0] = 0.0; - ucell.atoms[0].taud[1] = 0.0; - ucell.atoms[0].taud[2] = 0.0; - - ucell.atoms[0].tau.resize(nat); - - ucell.lat_axis_free[0] = 1; - ucell.lat_axis_free[1] = 1; - ucell.lat_axis_free[2] = 1; - - rl.init_relax(nat, inp); - rl.relax_step(ucell,force_in,stress_in,0.0, ofs); - - for(int i=0;i<3;i++) - { - result.push_back(ucell.atoms[0].taud[i].x); - result.push_back(ucell.atoms[0].taud[i].y); - result.push_back(ucell.atoms[0].taud[i].z); - } - push_result(); - - //reset lattice vector - ucell.latvec.Identity(); - inp.fixed_axes = "shape"; - rl.init_relax(nat, inp); - rl.relax_step(ucell,force_in,stress_in,0.0, ofs); - push_result(); - - //reset lattice vector - ucell.latvec.Identity(); - inp.fixed_axes = "volume"; - rl.init_relax(nat, inp); - rl.relax_step(ucell,force_in,stress_in,0.0, ofs); - push_result(); - - //reset lattice vector - ucell.latvec.Identity(); - inp.fixed_axes = "a"; //anything other than "None" - inp.fixed_ibrav = true; - ucell.latName = "sc"; - ucell.lat_axis_free[0] = 0; - ucell.lat_axis_free[1] = 0; - ucell.lat_axis_free[2] = 0; - rl.init_relax(nat, inp); - rl.relax_step(ucell,force_in,stress_in,0.0, ofs); - - push_result(); - } + protected: + Relax rl; + std::vector result; + Input_para inp; + UnitCell ucell; + std::ofstream ofs; - void push_result() + void SetUp() + { + inp.force_thr = 0.001; + inp.calculation = "cell-relax"; + + ModuleBase::matrix force_in, stress_in; + int nat = 3; + + force_in.create(nat, 3); + stress_in.create(3, 3); + + force_in(0, 0) = 0.1; + force_in(0, 1) = 0.1; + force_in(0, 2) = 0.1; + force_in(1, 0) = 0; + force_in(1, 1) = 0.1; + force_in(1, 2) = 0.1; + force_in(2, 0) = 0; + force_in(2, 1) = 0; + force_in(2, 2) = 0.1; + + stress_in(0, 0) = 1; + stress_in(0, 1) = 1; + stress_in(0, 2) = 1; + stress_in(1, 0) = 0; + stress_in(1, 1) = 1; + stress_in(1, 2) = 1; + stress_in(2, 0) = 0; + stress_in(2, 1) = 0; + stress_in(2, 2) = 1; + + ucell.ntype = 1; + ucell.nat = nat; + ucell.atoms = new Atom[1]; + ucell.atoms[0].na = nat; + ucell.omega = 1.0; + ucell.lat0 = 1.0; + + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); + ucell.atoms[0].mbl.resize(nat); + ucell.atoms[0].taud.resize(nat); + ucell.atoms[0].tau.resize(nat); + ucell.atoms[0].dis.resize(nat); + ucell.atoms[0].mag.resize(nat); + ucell.atoms[0].vel.resize(nat); + + ucell.iat2it[0] = 0; + ucell.iat2it[1] = 0; + ucell.iat2it[2] = 0; + + ucell.iat2ia[0] = 0; + ucell.iat2ia[1] = 1; + ucell.iat2ia[2] = 2; + + ucell.atoms[0].mbl[0].x = 0; + ucell.atoms[0].mbl[0].y = 0; + ucell.atoms[0].mbl[0].z = 1; + + ucell.atoms[0].mbl[1].x = 0; + ucell.atoms[0].mbl[1].y = 1; + ucell.atoms[0].mbl[1].z = 0; + + ucell.atoms[0].mbl[2].x = 1; + ucell.atoms[0].mbl[2].y = 0; + ucell.atoms[0].mbl[2].z = 0; + + ucell.atoms[0].taud[0] = 0.0; + ucell.atoms[0].taud[1] = 0.0; + ucell.atoms[0].taud[2] = 0.0; + + ucell.atoms[0].tau.resize(nat); + + ucell.lat_axis_free[0] = 1; + ucell.lat_axis_free[1] = 1; + ucell.lat_axis_free[2] = 1; + + rl.init_relax(nat, inp); + rl.relax_step(ucell, force_in, stress_in, 0.0, ofs); + + for (int i = 0; i < 3; i++) { - result.push_back(ucell.latvec.e11); - result.push_back(ucell.latvec.e12); - result.push_back(ucell.latvec.e13); - result.push_back(ucell.latvec.e21); - result.push_back(ucell.latvec.e22); - result.push_back(ucell.latvec.e23); - result.push_back(ucell.latvec.e31); - result.push_back(ucell.latvec.e32); - result.push_back(ucell.latvec.e33); + result.push_back(ucell.atoms[0].taud[i].x); + result.push_back(ucell.atoms[0].taud[i].y); + result.push_back(ucell.atoms[0].taud[i].z); } + push_result(); + + // reset lattice vector + ucell.latvec.Identity(); + inp.fixed_axes = "shape"; + rl.init_relax(nat, inp); + rl.relax_step(ucell, force_in, stress_in, 0.0, ofs); + push_result(); + + // reset lattice vector + ucell.latvec.Identity(); + inp.fixed_axes = "volume"; + rl.init_relax(nat, inp); + rl.relax_step(ucell, force_in, stress_in, 0.0, ofs); + push_result(); + + // reset lattice vector + ucell.latvec.Identity(); + inp.fixed_axes = "a"; // anything other than "None" + inp.fixed_ibrav = true; + ucell.latName = "sc"; + ucell.lat_axis_free[0] = 0; + ucell.lat_axis_free[1] = 0; + ucell.lat_axis_free[2] = 0; + rl.init_relax(nat, inp); + rl.relax_step(ucell, force_in, stress_in, 0.0, ofs); + + push_result(); + } + void push_result() + { + result.push_back(ucell.latvec.e11); + result.push_back(ucell.latvec.e12); + result.push_back(ucell.latvec.e13); + result.push_back(ucell.latvec.e21); + result.push_back(ucell.latvec.e22); + result.push_back(ucell.latvec.e23); + result.push_back(ucell.latvec.e31); + result.push_back(ucell.latvec.e32); + result.push_back(ucell.latvec.e33); + } }; TEST_F(Test_SETGRAD, relax_new) { - std::vector result_ref = - { - 0, 0, 0.24293434145, - 0, 0.242934341453, 0, - 0, 0, 0, - //paramter for taud - 1.2267616333,0.2267616333,0.22676163333, - 0, 1.2267616333 ,0.2267616333, - 0, 0, 1.22676163333, + std::vector result_ref = { + 0, + 0, + 0.24293434145, + 0, + 0.242934341453, + 0, + 0, + 0, + 0, + // paramter for taud + 1.2267616333, + 0.2267616333, + 0.22676163333, + 0, + 1.2267616333, + 0.2267616333, + 0, + 0, + 1.22676163333, // paramter for fisrt time after relaxation - 1.3677603495, 0, 0, - 0, 1.36776034956, 0, - 0, 0, 1.36776034956, + 1.3677603495, + 0, + 0, + 0, + 1.36776034956, + 0, + 0, + 0, + 1.36776034956, // paramter for second time after relaxation - 1.3677603495 ,0.3633367476,0.36333674766, - 0, 1.3677603495 ,0.36333674766, - 0, 0, 1.3677603495 , + 1.3677603495, + 0.3633367476, + 0.36333674766, + 0, + 1.3677603495, + 0.36333674766, + 0, + 0, + 1.3677603495, // paramter for third time after relaxation - 1,0,0,0,1,0,0,0,1 + 1, + 0, + 0, + 0, + 1, + 0, + 0, + 0, + 1 // paramter for fourth time after relaxation }; - for(int i=0;i result; - Input_para inp; - UnitCell ucell; - std::ofstream ofs; - - void SetUp() - { - int nstep = 3; - int nat = 5; - inp.calculation = "cell-relax"; - inp.force_thr = 0.001; - inp.fixed_axes = "a"; - inp.stress_thr = 0.01; - inp.fixed_ibrav = false; + protected: + Relax rl; + std::vector result; + Input_para inp; + UnitCell ucell; + std::ofstream ofs; + + void SetUp() + { + int nstep = 3; + int nat = 5; + inp.calculation = "cell-relax"; + inp.force_thr = 0.001; + inp.fixed_axes = "a"; + inp.stress_thr = 0.01; + inp.fixed_ibrav = false; - this->setup_cell(); + this->setup_cell(); - ModuleBase::matrix force_in, stress_in; - force_in.create(nat,3); - stress_in.create(3,3); + ModuleBase::matrix force_in, stress_in; + force_in.create(nat, 3); + stress_in.create(3, 3); - rl.init_relax(nat, inp); + rl.init_relax(nat, inp); - for(int istep=0;istep result_ref= - { - 0.5000000586,0.4999998876,0.009364595811, - 0.9999999281,1.901333279e-07,0.5476035454, - 3.782097706e-07,0.4999999285,0.4770375874, - 0.4999997826,2.072311863e-07,0.477037871, - 0.9999998523,0.9999997866,0.9349574003, + std::vector result_ref = { + 0.5000000586, + 0.4999998876, + 0.009364595811, + 0.9999999281, + 1.901333279e-07, + 0.5476035454, + 3.782097706e-07, + 0.4999999285, + 0.4770375874, + 0.4999997826, + 2.072311863e-07, + 0.477037871, + 0.9999998523, + 0.9999997866, + 0.9349574003, // paramter for taud after first relaxation - 4.006349654,-1.93128788e-07,5.793863639e-07, - -1.93128788e-07,4.006354579,-3.86257576e-07, - 6.757962549e-07,-4.505308366e-07,3.966870038, + 4.006349654, + -1.93128788e-07, + 5.793863639e-07, + -1.93128788e-07, + 4.006354579, + -3.86257576e-07, + 6.757962549e-07, + -4.505308366e-07, + 3.966870038, // paramter for latvec after first relaxation - 0.5000000566,0.4999998916,0.009177239183, - 0.9999999308,1.832935626e-07,0.5467689737, - 3.647769323e-07,0.4999999311,0.4771204124, - 0.4999997903,1.998879545e-07,0.4771206859, - 0.9999998574,0.9999997943,0.9358136888, + 0.5000000566, + 0.4999998916, + 0.009177239183, + 0.9999999308, + 1.832935626e-07, + 0.5467689737, + 3.647769323e-07, + 0.4999999311, + 0.4771204124, + 0.4999997903, + 1.998879545e-07, + 0.4771206859, + 0.9999998574, + 0.9999997943, + 0.9358136888, // paramter for taud after second relaxation - 3.999761277,-1.656764727e-07,4.97029418e-07, - -1.656764727e-07,3.999765501,-3.313529453e-07, - 5.797351131e-07,-3.864900754e-07,4.010925071, + 3.999761277, + -1.656764727e-07, + 4.97029418e-07, + -1.656764727e-07, + 3.999765501, + -3.313529453e-07, + 5.797351131e-07, + -3.864900754e-07, + 4.010925071, // paramter for latvec after second relaxation - 0.500000082,0.4999999574,0.01057784352, - 0.9999999149,1.939640249e-07,0.5455830599, - 3.795967781e-07,0.4999998795,0.4765373919, - 0.4999998037,2.756298268e-07,0.4765374602, - 0.9999998196,0.9999996936,0.9367652445, + 0.500000082, + 0.4999999574, + 0.01057784352, + 0.9999999149, + 1.939640249e-07, + 0.5455830599, + 3.795967781e-07, + 0.4999998795, + 0.4765373919, + 0.4999998037, + 2.756298268e-07, + 0.4765374602, + 0.9999998196, + 0.9999996936, + 0.9367652445, // paramter for taud after third relaxation - 4.017733155,-1.420363309e-07,2.637046077e-07, - -1.420364243e-07,4.017735987,3.126225134e-07, - 3.479123171e-07,2.578467568e-07,4.011674933 + 4.017733155, + -1.420363309e-07, + 2.637046077e-07, + -1.420364243e-07, + 4.017735987, + 3.126225134e-07, + 3.479123171e-07, + 2.578467568e-07, + 4.011674933 // paramter for latvec after third relaxation }; - for(int i=0;i (type0, ia0); 1 -> (type1, ia0); 2 -> (type1, ia1) ucell.iat2it.resize(3); ucell.iat2ia.resize(3); - ucell.iat2it[0] = 0; ucell.iat2it[1] = 1; ucell.iat2it[2] = 1; - ucell.iat2ia[0] = 0; ucell.iat2ia[1] = 0; ucell.iat2ia[2] = 1; + ucell.iat2it[0] = 0; + ucell.iat2it[1] = 1; + ucell.iat2it[2] = 1; + ucell.iat2ia[0] = 0; + ucell.iat2ia[1] = 0; + ucell.iat2ia[2] = 1; // mbl flags, laid out per atom ucell.atoms[0].mbl[0] = {mbl_flat[0], mbl_flat[1], mbl_flat[2]}; @@ -655,22 +736,24 @@ void free_two_type_cell(UnitCell& ucell) // lattice scales isotropically (off-diagonals stay zero, diagonal scales equal). TEST(RelaxSyncCellMove, FixedShapeIsotropicScale) { - const int mbl[9] = {1,1,1, 1,1,1, 1,1,1}; + const int mbl[9] = {1, 1, 1, 1, 1, 1, 1, 1, 1}; UnitCell ucell; make_two_type_cell(ucell, mbl); Input_para inp; inp.calculation = "cell-relax"; - inp.force_thr = 0.0; // force never converged -> take the CG move + inp.force_thr = 0.0; // force never converged -> take the CG move inp.force_thr_ev = 0.0; - inp.stress_thr = 0.0; // stress never converged either + inp.stress_thr = 0.0; // stress never converged either inp.fixed_axes = "shape"; inp.fixed_ibrav = false; ModuleBase::matrix force_in(3, 3); - force_in(0, 0) = 0.01; // nonzero -> not converged -> performs the move + force_in(0, 0) = 0.01; // nonzero -> not converged -> performs the move ModuleBase::matrix stress_in(3, 3); - stress_in(0, 0) = 1.0; stress_in(1, 1) = 2.0; stress_in(2, 2) = 3.0; + stress_in(0, 0) = 1.0; + stress_in(1, 1) = 2.0; + stress_in(2, 2) = 3.0; stress_in(0, 1) = 0.5; // anisotropic + off-diagonal parts must be dropped Relax rl; @@ -712,7 +795,7 @@ TEST(RelaxSyncMaxStep, GeometryMovedAndNotConverged) Input_para inp; inp.calculation = "relax"; inp.relax_method = "cg"; - inp.force_thr = 0.0; // never converged -> always take the move + inp.force_thr = 0.0; // never converged -> always take the move inp.force_thr_ev = 0.0; inp.fixed_axes = "None"; inp.fixed_ibrav = false; @@ -741,7 +824,7 @@ TEST(RelaxSyncMaxStep, GeometryMovedAndNotConverged) ModuleBase::matrix force_in(nat, 3); ModuleBase::matrix stress_in(3, 3); - force_in(0, 0) = 0.01; // nonzero -> not converged -> move + force_in(0, 0) = 0.01; // nonzero -> not converged -> move Relax rl; rl.init_relax(nat, inp); @@ -762,7 +845,7 @@ TEST(RelaxSyncMaxStep, GeometryMovedAndNotConverged) // back to the original omega after each move (relax_sync.cpp Step 1). TEST(RelaxSyncCellMove, FixedVolumePreservesOmega) { - const int mbl[9] = {1,1,1, 1,1,1, 1,1,1}; + const int mbl[9] = {1, 1, 1, 1, 1, 1, 1, 1, 1}; UnitCell ucell; make_two_type_cell(ucell, mbl); const double omega0 = ucell.omega; @@ -776,9 +859,11 @@ TEST(RelaxSyncCellMove, FixedVolumePreservesOmega) inp.fixed_ibrav = false; ModuleBase::matrix force_in(3, 3); - force_in(0, 0) = 0.01; // nonzero -> not converged + force_in(0, 0) = 0.01; // nonzero -> not converged ModuleBase::matrix stress_in(3, 3); - stress_in(0, 0) = 1.0; stress_in(1, 1) = -0.5; stress_in(0, 1) = 0.3; + stress_in(0, 0) = 1.0; + stress_in(1, 1) = -0.5; + stress_in(0, 1) = 0.3; Relax rl; rl.init_relax(3, inp); @@ -797,7 +882,7 @@ TEST(RelaxSyncCellMove, FixedVolumePreservesOmega) // must remain exactly at their initial values (relax_sync.cpp 572-589). TEST(RelaxSyncCellMove, PartialAxisFreeKeepsLockedRows) { - const int mbl[9] = {1,1,1, 1,1,1, 1,1,1}; + const int mbl[9] = {1, 1, 1, 1, 1, 1, 1, 1, 1}; UnitCell ucell; make_two_type_cell(ucell, mbl); ucell.lat_axis_free[0] = 1; @@ -813,10 +898,13 @@ TEST(RelaxSyncCellMove, PartialAxisFreeKeepsLockedRows) inp.fixed_ibrav = false; ModuleBase::matrix force_in(3, 3); - force_in(0, 0) = 0.01; // nonzero -> not converged + force_in(0, 0) = 0.01; // nonzero -> not converged ModuleBase::matrix stress_in(3, 3); - stress_in(0, 0) = 1.0; stress_in(1, 1) = 1.0; stress_in(2, 2) = 1.0; - stress_in(0, 1) = 0.2; stress_in(1, 0) = 0.2; + stress_in(0, 0) = 1.0; + stress_in(1, 1) = 1.0; + stress_in(2, 2) = 1.0; + stress_in(0, 1) = 0.2; + stress_in(1, 0) = 0.2; Relax rl; rl.init_relax(3, inp); @@ -842,12 +930,12 @@ TEST(RelaxSyncCellMove, PartialAxisFreeKeepsLockedRows) TEST(RelaxSyncCellMove, MixedMblMovesOnlyFreeComponents) { // atom0: free x only; atom1(type1,ia0): free y only; atom2(type1,ia1): free z only - const int mbl[9] = {1,0,0, 0,1,0, 0,0,1}; + const int mbl[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; UnitCell ucell; make_two_type_cell(ucell, mbl); Input_para inp; - inp.calculation = "relax"; // ions only, keep lattice fixed for clarity + inp.calculation = "relax"; // ions only, keep lattice fixed for clarity inp.force_thr = 0.0; inp.force_thr_ev = 0.0; inp.fixed_axes = "None"; diff --git a/source/source_relax/unittests/test_socket_driver.cpp b/source/source_relax/unittests/test_socket_driver.cpp index 95b338cd388..50e5fe25597 100644 --- a/source/source_relax/unittests/test_socket_driver.cpp +++ b/source/source_relax/unittests/test_socket_driver.cpp @@ -1,13 +1,12 @@ -#include "source_relax/socket_driver.h" - -#include "gmock/gmock.h" -#include "gtest/gtest.h" +#include "for_test.h" #include "mpi.h" #include "source_cell/unitcell.h" #include "source_esolver/esolver.h" #include "source_io/module_parameter/input_parameter.h" -#include "for_test.h" +#include "source_relax/socket_driver.h" +#include "gmock/gmock.h" +#include "gtest/gtest.h" #include #include #include @@ -277,9 +276,7 @@ void initialize_one_atom_cell(UnitCell& ucell) void send_positions(const int fd, const double x) { - const double identity[9] = {1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0}; + const double identity[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; const std::int32_t nat = 1; const double position[3] = {x, 0.0, 0.0}; send_header(fd, "POSDATA"); @@ -343,8 +340,7 @@ ForceResponse read_force_response(const int fd) } response.energy_hartree = read_value(fd); response.nat = read_value(fd); - response.forces_hartree_per_bohr - = read_doubles(fd, static_cast(3 * response.nat)); + response.forces_hartree_per_bohr = read_doubles(fd, static_cast(3 * response.nat)); response.virial_wire_hartree = read_doubles(fd, 9); const std::int32_t extra_bytes = read_value(fd); if (extra_bytes < 0) @@ -476,12 +472,10 @@ DriverResult run_driver_frame(const bool converged, TEST(SocketDriverTest, NonconvergedFrameIsPublishedWithMetadata) { ForceResponse response; - const DriverResult result = run_driver_frame( - false, true, false, - [&](const int fd) { - send_header(fd, "GETFORCE"); - response = read_force_response(fd); - }); + const DriverResult result = run_driver_frame(false, true, false, [&](const int fd) { + send_header(fd, "GETFORCE"); + response = read_force_response(fd); + }); EXPECT_EQ("FORCEREADY", response.header); EXPECT_EQ(0, result.exit_code); @@ -491,35 +485,27 @@ TEST(SocketDriverTest, NonconvergedFrameIsPublishedWithMetadata) TEST(SocketDriverTest, EnergyOnlyFrameMarksForceAndStressAbsent) { ForceResponse response; - const DriverResult result = run_driver_frame( - true, false, false, - [&](const int fd) { - send_header(fd, "GETFORCE"); - response = read_force_response(fd); - }); + const DriverResult result = run_driver_frame(true, false, false, [&](const int fd) { + send_header(fd, "GETFORCE"); + response = read_force_response(fd); + }); EXPECT_EQ("FORCEREADY", response.header); EXPECT_EQ(0, result.exit_code); EXPECT_THAT(response.extra, testing::HasSubstr("\"present\":[\"energy\"]")); EXPECT_THAT(response.extra, testing::Not(testing::HasSubstr("\"forces\""))); EXPECT_THAT(response.extra, testing::Not(testing::HasSubstr("\"stress\""))); - EXPECT_THAT(response.forces_hartree_per_bohr, - testing::ElementsAre(0.0, 0.0, 0.0)); - EXPECT_THAT(response.virial_wire_hartree, - testing::ElementsAre(0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0)); + EXPECT_THAT(response.forces_hartree_per_bohr, testing::ElementsAre(0.0, 0.0, 0.0)); + EXPECT_THAT(response.virial_wire_hartree, testing::ElementsAre(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)); } TEST(SocketDriverTest, EnergyAndStressFrameDoesNotAdvertiseForce) { ForceResponse response; - const DriverResult result = run_driver_frame( - true, false, true, - [&](const int fd) { - send_header(fd, "GETFORCE"); - response = read_force_response(fd); - }); + const DriverResult result = run_driver_frame(true, false, true, [&](const int fd) { + send_header(fd, "GETFORCE"); + response = read_force_response(fd); + }); EXPECT_EQ("FORCEREADY", response.header); EXPECT_EQ(0, result.exit_code); @@ -531,42 +517,32 @@ TEST(SocketDriverTest, EnergyAndStressFrameDoesNotAdvertiseForce) TEST(SocketDriverTest, EnergyAndForceFrameAdvertisesOnlyForce) { ForceResponse response; - const DriverResult result = run_driver_frame( - true, true, false, - [&](const int fd) { - send_header(fd, "GETFORCE"); - response = read_force_response(fd); - }); + const DriverResult result = run_driver_frame(true, true, false, [&](const int fd) { + send_header(fd, "GETFORCE"); + response = read_force_response(fd); + }); EXPECT_EQ("FORCEREADY", response.header); EXPECT_EQ(0, result.exit_code); EXPECT_THAT(response.extra, testing::HasSubstr("\"present\":[\"energy\",\"forces\"]")); EXPECT_THAT(response.extra, testing::Not(testing::HasSubstr("\"stress\""))); - EXPECT_THAT(response.forces_hartree_per_bohr, - testing::ElementsAre(2.0, 0.0, 0.0)); - EXPECT_THAT(response.virial_wire_hartree, - testing::ElementsAre(0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0)); + EXPECT_THAT(response.forces_hartree_per_bohr, testing::ElementsAre(2.0, 0.0, 0.0)); + EXPECT_THAT(response.virial_wire_hartree, testing::ElementsAre(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)); } TEST(SocketDriverTest, EnergyForceAndStressFrameAdvertisesBothDerivatives) { ForceResponse response; - const DriverResult result = run_driver_frame( - true, true, true, - [&](const int fd) { - send_header(fd, "GETFORCE"); - response = read_force_response(fd); - }); + const DriverResult result = run_driver_frame(true, true, true, [&](const int fd) { + send_header(fd, "GETFORCE"); + response = read_force_response(fd); + }); EXPECT_EQ("FORCEREADY", response.header); EXPECT_EQ(0, result.exit_code); - EXPECT_THAT(response.extra, - testing::HasSubstr("\"present\":[\"energy\",\"forces\",\"stress\"]")); + EXPECT_THAT(response.extra, testing::HasSubstr("\"present\":[\"energy\",\"forces\",\"stress\"]")); EXPECT_EQ(3u, response.forces_hartree_per_bohr.size()); - EXPECT_THAT(response.forces_hartree_per_bohr, - testing::ElementsAre(2.0, 0.0, 0.0)); + EXPECT_THAT(response.forces_hartree_per_bohr, testing::ElementsAre(2.0, 0.0, 0.0)); EXPECT_NE(0.0, response.virial_wire_hartree[0]); } diff --git a/source/source_relax/unittests/test_socket_frame.cpp b/source/source_relax/unittests/test_socket_frame.cpp index e420cab824f..9d9dd684514 100644 --- a/source/source_relax/unittests/test_socket_frame.cpp +++ b/source/source_relax/unittests/test_socket_frame.cpp @@ -1,7 +1,6 @@ #include "../socket_frame.h" #include "gtest/gtest.h" - #include #include #include @@ -11,13 +10,13 @@ namespace { using SocketFrame::CellValidation; -using SocketFrame::Matrix9; -using SocketFrame::VirialConversion; using SocketFrame::checked_position_count; using SocketFrame::make_ipi_virial; +using SocketFrame::Matrix9; using SocketFrame::transpose_matrix9; using SocketFrame::validate_ipi_cell; using SocketFrame::validate_positions; +using SocketFrame::VirialConversion; const double EPSILON = std::numeric_limits::epsilon(); @@ -54,9 +53,7 @@ TEST(SocketFrameTest, VirialUsesPositiveHalfVolumeAndWireTranspose) TEST(SocketFrameTest, RightHandedTriclinicCellReturnsKnownInverse) { const Matrix9 cell = {{2.0, 1.0, 0.0, 0.0, 3.0, 1.0, 0.0, 0.0, 4.0}}; - const Matrix9 inverse = {{0.5, -1.0 / 6.0, 1.0 / 24.0, - 0.0, 1.0 / 3.0, -1.0 / 12.0, - 0.0, 0.0, 0.25}}; + const Matrix9 inverse = {{0.5, -1.0 / 6.0, 1.0 / 24.0, 0.0, 1.0 / 3.0, -1.0 / 12.0, 0.0, 0.0, 0.25}}; const CellValidation out = validate_with_driver_thresholds(cell, inverse); @@ -72,41 +69,46 @@ TEST(SocketFrameTest, AseTriclinicInverseWireLayoutIsAcceptedAndRecomputedFromCe // and sends pinv(A) * Bohr as the inverse field. For this nonsingular // cell that received field is inv(H)^T, not inv(H). const double bohr_angstrom = 0.5291772105638411; - const Matrix9 cell_wire = {{5.0 / bohr_angstrom, 0.5 / bohr_angstrom, 0.25 / bohr_angstrom, - 0.0, 4.0 / bohr_angstrom, 0.75 / bohr_angstrom, - 0.0, 0.0, 3.0 / bohr_angstrom}}; - const Matrix9 ase_inverse_wire = {{bohr_angstrom / 5.0, 0.0, 0.0, - -bohr_angstrom / 40.0, bohr_angstrom / 4.0, 0.0, - -bohr_angstrom / 96.0, -bohr_angstrom / 16.0, + const Matrix9 cell_wire = {{5.0 / bohr_angstrom, + 0.5 / bohr_angstrom, + 0.25 / bohr_angstrom, + 0.0, + 4.0 / bohr_angstrom, + 0.75 / bohr_angstrom, + 0.0, + 0.0, + 3.0 / bohr_angstrom}}; + const Matrix9 ase_inverse_wire = {{bohr_angstrom / 5.0, + 0.0, + 0.0, + -bohr_angstrom / 40.0, + bohr_angstrom / 4.0, + 0.0, + -bohr_angstrom / 96.0, + -bohr_angstrom / 16.0, bohr_angstrom / 3.0}}; const Matrix9 inverse_computed_from_cell = {{bohr_angstrom / 5.0, - -bohr_angstrom / 40.0, - -bohr_angstrom / 96.0, - 0.0, - bohr_angstrom / 4.0, - -bohr_angstrom / 16.0, - 0.0, - 0.0, - bohr_angstrom / 3.0}}; + -bohr_angstrom / 40.0, + -bohr_angstrom / 96.0, + 0.0, + bohr_angstrom / 4.0, + -bohr_angstrom / 16.0, + 0.0, + 0.0, + bohr_angstrom / 3.0}}; const CellValidation out = validate_with_driver_thresholds(cell_wire, ase_inverse_wire); ASSERT_TRUE(out.ok) << out.message; EXPECT_NEAR(0.0, out.inverse_residual, 16.0 * EPSILON); - expect_matrix_near(inverse_computed_from_cell, - out.computed_inverse_wire_bohr_inv, - 16.0 * EPSILON); + expect_matrix_near(inverse_computed_from_cell, out.computed_inverse_wire_bohr_inv, 16.0 * EPSILON); } TEST(SocketFrameTest, RotatedDiagonalTracksRightSingularVectorsAndInverseOrder) { // Hand-multiplied U diag(5, 2, 0.5) V^T, with rational plane rotations. - const Matrix9 cell = {{4.0, -0.72, -0.96, - 3.0, 0.96, 1.28, - 0.0, -0.4, 0.3}}; - const Matrix9 inverse = {{0.16, 0.12, 0.0, - -0.18, 0.24, -1.6, - -0.24, 0.32, 1.2}}; + const Matrix9 cell = {{4.0, -0.72, -0.96, 3.0, 0.96, 1.28, 0.0, -0.4, 0.3}}; + const Matrix9 inverse = {{0.16, 0.12, 0.0, -0.18, 0.24, -1.6, -0.24, 0.32, 1.2}}; const CellValidation out = validate_with_driver_thresholds(cell, inverse); @@ -121,9 +123,7 @@ TEST(SocketFrameTest, InconsistentReceivedInverseIsRejected) const Matrix9 cell = {{2.0, 1.0, 0.0, 0.0, 3.0, 1.0, 0.0, 0.0, 4.0}}; // This is neither inv(cell) nor inv(cell)^T, so both supported wire // layouts must reject it. - const Matrix9 wrong_inverse = {{0.6, -1.0 / 6.0, 1.0 / 24.0, - 0.0, 1.0 / 3.0, -1.0 / 12.0, - 0.0, 0.0, 0.25}}; + const Matrix9 wrong_inverse = {{0.6, -1.0 / 6.0, 1.0 / 24.0, 0.0, 1.0 / 3.0, -1.0 / 12.0, 0.0, 0.0, 0.25}}; const CellValidation out = validate_with_driver_thresholds(cell, wrong_inverse); @@ -134,19 +134,13 @@ TEST(SocketFrameTest, InconsistentReceivedInverseIsRejected) TEST(SocketFrameTest, ReceivedInverseResidualUsesConditionScaledRelativeTolerance) { - const Matrix9 cell = {{1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0e-6}}; - Matrix9 accepted_inverse = {{1.0 + 1.0e-8, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0e6}}; + const Matrix9 cell = {{1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0e-6}}; + Matrix9 accepted_inverse = {{1.0 + 1.0e-8, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0e6}}; Matrix9 rejected_inverse = accepted_inverse; rejected_inverse[0] = 1.0 + 2.0e-8; - const CellValidation accepted - = validate_ipi_cell(cell, accepted_inverse, 1.0e12, 0.0, 64.0); - const CellValidation rejected - = validate_ipi_cell(cell, rejected_inverse, 1.0e12, 0.0, 64.0); + const CellValidation accepted = validate_ipi_cell(cell, accepted_inverse, 1.0e12, 0.0, 64.0); + const CellValidation rejected = validate_ipi_cell(cell, rejected_inverse, 1.0e12, 0.0, 64.0); ASSERT_TRUE(accepted.ok) << accepted.message; EXPECT_DOUBLE_EQ(1.0e6, accepted.condition_number_2); @@ -168,12 +162,8 @@ TEST(SocketFrameTest, NegativeAndZeroDeterminantsAreRejected) TEST(SocketFrameTest, NonrepresentablePositiveCellVolumeIsRejected) { - const Matrix9 huge_cell = {{1.0e200, 0.0, 0.0, - 0.0, 1.0e200, 0.0, - 0.0, 0.0, 1.0e200}}; - const Matrix9 tiny_inverse = {{1.0e-200, 0.0, 0.0, - 0.0, 1.0e-200, 0.0, - 0.0, 0.0, 1.0e-200}}; + const Matrix9 huge_cell = {{1.0e200, 0.0, 0.0, 0.0, 1.0e200, 0.0, 0.0, 0.0, 1.0e200}}; + const Matrix9 tiny_inverse = {{1.0e-200, 0.0, 0.0, 0.0, 1.0e-200, 0.0, 0.0, 0.0, 1.0e-200}}; const CellValidation out = validate_with_driver_thresholds(huge_cell, tiny_inverse); @@ -183,12 +173,8 @@ TEST(SocketFrameTest, NonrepresentablePositiveCellVolumeIsRejected) TEST(SocketFrameTest, UnderflowedCellVolumeIsRejectedAsZeroDeterminant) { - const Matrix9 tiny_cell = {{1.0e-200, 0.0, 0.0, - 0.0, 1.0e-200, 0.0, - 0.0, 0.0, 1.0e-200}}; - const Matrix9 huge_inverse = {{1.0e200, 0.0, 0.0, - 0.0, 1.0e200, 0.0, - 0.0, 0.0, 1.0e200}}; + const Matrix9 tiny_cell = {{1.0e-200, 0.0, 0.0, 0.0, 1.0e-200, 0.0, 0.0, 0.0, 1.0e-200}}; + const Matrix9 huge_inverse = {{1.0e200, 0.0, 0.0, 0.0, 1.0e200, 0.0, 0.0, 0.0, 1.0e200}}; const CellValidation out = validate_with_driver_thresholds(tiny_cell, huge_inverse); @@ -257,10 +243,7 @@ TEST(SocketFrameTest, PositionCountRejectsMismatchBeforeDerivingAllocationSize) std::size_t coordinate_count = 123; std::string message; - EXPECT_FALSE(checked_position_count(std::numeric_limits::max(), - 1, - coordinate_count, - message)); + EXPECT_FALSE(checked_position_count(std::numeric_limits::max(), 1, coordinate_count, message)); EXPECT_EQ(123u, coordinate_count); EXPECT_NE(std::string::npos, message.find("match")); } @@ -289,12 +272,8 @@ TEST(SocketFrameTest, PositionsRequireExactSizeAndFiniteCoordinates) TEST(SocketFrameTest, SmallStressAsymmetryIsAveragedBeforeConversion) { - const Matrix9 stress = {{1.0, 2.1, 3.2, - 1.9, 5.0, 6.3, - 2.8, 5.7, 9.0}}; - const Matrix9 expected = {{1.0, 2.0, 3.0, - 2.0, 5.0, 6.0, - 3.0, 6.0, 9.0}}; + const Matrix9 stress = {{1.0, 2.1, 3.2, 1.9, 5.0, 6.3, 2.8, 5.7, 9.0}}; + const Matrix9 expected = {{1.0, 2.0, 3.0, 2.0, 5.0, 6.0, 3.0, 6.0, 9.0}}; const VirialConversion out = make_ipi_virial(stress, 2.0, 0.61, 0.0); @@ -305,9 +284,7 @@ TEST(SocketFrameTest, SmallStressAsymmetryIsAveragedBeforeConversion) TEST(SocketFrameTest, ExcessiveStressAsymmetryIsRejected) { - const Matrix9 stress = {{1.0, 2.1, 3.2, - 1.9, 5.0, 6.3, - 2.8, 5.7, 9.0}}; + const Matrix9 stress = {{1.0, 2.1, 3.2, 1.9, 5.0, 6.3, 2.8, 5.7, 9.0}}; const VirialConversion out = make_ipi_virial(stress, 2.0, 0.59, 0.0); @@ -318,20 +295,14 @@ TEST(SocketFrameTest, ExcessiveStressAsymmetryIsRejected) TEST(SocketFrameTest, StressAsymmetryUsesAbsolutePlusRelativeTolerance) { - const Matrix9 accepted_stress = {{10.0, 2.0 + 4.0e-8, 3.0, - 2.0 - 4.0e-8, 5.0, 6.0, - 3.0, 6.0, 9.0}}; + const Matrix9 accepted_stress = {{10.0, 2.0 + 4.0e-8, 3.0, 2.0 - 4.0e-8, 5.0, 6.0, 3.0, 6.0, 9.0}}; Matrix9 rejected_stress = accepted_stress; rejected_stress[1] = 2.0 + 6.0e-8; rejected_stress[3] = 2.0 - 6.0e-8; - const Matrix9 expected = {{10.0, 2.0, 3.0, - 2.0, 5.0, 6.0, - 3.0, 6.0, 9.0}}; + const Matrix9 expected = {{10.0, 2.0, 3.0, 2.0, 5.0, 6.0, 3.0, 6.0, 9.0}}; - const VirialConversion accepted - = make_ipi_virial(accepted_stress, 2.0, 1.0e-10, 1.0e-8); - const VirialConversion rejected - = make_ipi_virial(rejected_stress, 2.0, 1.0e-10, 1.0e-8); + const VirialConversion accepted = make_ipi_virial(accepted_stress, 2.0, 1.0e-10, 1.0e-8); + const VirialConversion rejected = make_ipi_virial(rejected_stress, 2.0, 1.0e-10, 1.0e-8); ASSERT_TRUE(accepted.ok) << accepted.message; expect_matrix_near(expected, accepted.wire_virial_hartree, 4.0 * EPSILON); @@ -347,19 +318,13 @@ TEST(SocketFrameTest, NonpositiveOrNonfiniteVolumeIsRejected) EXPECT_FALSE(make_ipi_virial(zero_stress, 0.0, 1.0e-10, 1.0e-8).ok); EXPECT_FALSE(make_ipi_virial(zero_stress, -1.0, 1.0e-10, 1.0e-8).ok); - EXPECT_FALSE(make_ipi_virial(zero_stress, - std::numeric_limits::infinity(), - 1.0e-10, - 1.0e-8) - .ok); + EXPECT_FALSE(make_ipi_virial(zero_stress, std::numeric_limits::infinity(), 1.0e-10, 1.0e-8).ok); } TEST(SocketFrameTest, FiniteStressAndVolumeRejectConvertedVirialOverflow) { const double largest_finite = std::numeric_limits::max(); - const Matrix9 stress = {{largest_finite, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0}}; + const Matrix9 stress = {{largest_finite, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}}; const VirialConversion out = make_ipi_virial(stress, 4.0, 1.0e-10, 1.0e-8); diff --git a/source/source_relax/unittests/test_socket_ipi.cpp b/source/source_relax/unittests/test_socket_ipi.cpp index 2f02e1cd55a..6b4420c2c70 100644 --- a/source/source_relax/unittests/test_socket_ipi.cpp +++ b/source/source_relax/unittests/test_socket_ipi.cpp @@ -1,11 +1,10 @@ #include "../socket_ipi.h" #include "gtest/gtest.h" - #include #include -#include #include +#include #include #include #include