diff --git a/.clang-format b/.clang-format index 8970df99c2..5ea3e06441 100644 --- a/.clang-format +++ b/.clang-format @@ -1,49 +1,63 @@ --- Language: Cpp -BasedOnStyle: Microsoft +BasedOnStyle: LLVM -AlwaysBreakTemplateDeclarations: Yes +ColumnLimit: 120 -AllowAllArgumentsOnNextLine: false -AllowAllParametersOfDeclarationOnNextLine: false +#AlwaysBreakTemplateDeclarations: Yes +# +#AllowAllArgumentsOnNextLine: false +#AllowAllParametersOfDeclarationOnNextLine: false BinPackArguments: false BinPackParameters: false - -BreakBeforeBinaryOperators: All -BreakBeforeTernaryOperators: true - -Cpp11BracedListStyle: true -FixNamespaceComments: true - -# InsertBraces: true # only for clang-format version15.0.0 or later - -# About include -IncludeBlocks: Regroup -IncludeIsMainRegex: '([-_](test|unittest))?$' -IncludeIsMainSourceRegex: '' - +PackArguments: + BinPack: UseBreakAfter + BreakAfter: 3 +PackParameters: + BinPack: UseBreakAfter + BreakAfter: 3 +# +#BreakBeforeBinaryOperators: None +#BreakBeforeTernaryOperators: true +BreakBeforeBraces: Allman +BreakAfterOpenBracketFunction: true +AlignAfterOpenBracket: DontAlign +# +#Cpp11BracedListStyle: true +#FixNamespaceComments: true +# +InsertBraces: true # only for clang-format version15.0.0 or later +# +## About include +#IncludeBlocks: Regroup +#IncludeIsMainRegex: '([-_](test|unittest))?$' +#IncludeIsMainSourceRegex: '' +# IndentWidth: 4 -IndentWrappedFunctionNames: true - -# About Point -DerivePointerAlignment: false -PointerAlignment: Left - -ReflowComments: true -SortIncludes: true -SortUsingDeclarations: true - -# About space -SpaceAfterCStyleCast: false -SpaceAfterLogicalNot: false -SpaceBeforeRangeBasedForLoopColon: false -SpaceInEmptyParentheses: false -SpacesInAngles: false -SpacesInCStyleCastParentheses: false -SpacesInContainerLiterals: false -SpacesInParentheses: false -SpacesInSquareBrackets: false +#IndentWrappedFunctionNames: true +# +## About Point +#DerivePointerAlignment: false +#PointerAlignment: Left +PenaltyBreakBeforeMemberAccess: 1000 +PenaltyBreakScopeResolution: 1000 +# +#ReflowComments: true +#SortIncludes: true +#SortUsingDeclarations: true +# +## About space +#SpaceAfterCStyleCast: false +#SpaceAfterLogicalNot: false +#SpaceBeforeRangeBasedForLoopColon: false +#SpaceInEmptyParentheses: false +#SpacesInAngles: false +#SpacesInCStyleCastParentheses: false +#SpacesInContainerLiterals: false +#SpacesInParentheses: false +#SpacesInSquareBrackets: false +# Standard: c++11 TabWidth: 4 UseTab: Never diff --git a/source/source_md/fire.cpp b/source/source_md/fire.cpp index cb4611c6e9..dbb7f0f639 100644 --- a/source/source_md/fire.cpp +++ b/source/source_md/fire.cpp @@ -6,13 +6,19 @@ #endif #include "source_base/timer.h" -FIRE::FIRE(const MD_para& mdp_in, - const bool cal_stress_in, - const bool init_vel, - const int my_rank_in, - const double force_thr_in, - MDCell& mdcell_in) - : MD_base(mdp_in, cal_stress_in, init_vel, my_rank_in, mdcell_in) +FIRE::FIRE( + const MD_para &mdp_in, + const bool cal_stress_in, + const bool init_vel, + const int my_rank_in, + const double force_thr_in, + MDCell &mdcell_in) + : MD_base( + mdp_in, + cal_stress_in, + init_vel, + my_rank_in, + mdcell_in) { force_thr = force_thr_in; dt_max = -1.0; @@ -32,11 +38,9 @@ FIRE::FIRE(const MD_para& mdp_in, // force_thr = 1e-3; } -FIRE::~FIRE() -{ -} +FIRE::~FIRE() {} -void FIRE::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp) +void FIRE::setup(ModuleESolver::ESolver *p_esolver, const std::string &global_readin_dir, DomainDecomposition &decomp) { ModuleBase::TITLE("FIRE", "setup"); ModuleBase::timer::start("FIRE", "setup"); @@ -50,7 +54,7 @@ void FIRE::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_re return; } -void FIRE::first_half(std::ofstream& ofs) +void FIRE::first_half(std::ofstream &ofs) { ModuleBase::TITLE("FIRE", "first_half"); ModuleBase::timer::start("FIRE", "first_half"); @@ -66,7 +70,6 @@ void FIRE::first_half(std::ofstream& ofs) return; } - void FIRE::second_half(void) { ModuleBase::TITLE("FIRE", "second_half"); @@ -81,8 +84,7 @@ void FIRE::second_half(void) return; } - -void FIRE::print_md(std::ofstream& ofs, const bool& cal_stress) +void FIRE::print_md(std::ofstream &ofs, const bool &cal_stress) { MD_base::print_md(ofs, cal_stress); @@ -94,8 +96,7 @@ void FIRE::print_md(std::ofstream& ofs, const bool& cal_stress) return; } - -void FIRE::write_restart(const std::string& global_out_dir) +void FIRE::write_restart(const std::string &global_out_dir) { if (!my_rank) { @@ -118,8 +119,7 @@ void FIRE::write_restart(const std::string& global_out_dir) return; } - -void FIRE::restart(const std::string& global_readin_dir) +void FIRE::restart(const std::string &global_readin_dir) { bool ok = true; @@ -137,16 +137,21 @@ void FIRE::restart(const std::string& global_readin_dir) if (ok) { file >> step_rst_ >> md_tfirst >> alpha >> negative_count >> dt_max >> md_dt; - if(!file) + if (!file) { ok = false; - } + } file.close(); } } #ifdef __MPI - MPI_Bcast(&ok, 1, MPI_C_BOOL, 0, mdcell.communicator()); + MPI_Bcast( + &ok, + 1, + MPI_C_BOOL, + 0, + mdcell.communicator()); #endif if (!ok) @@ -155,12 +160,42 @@ void FIRE::restart(const std::string& global_readin_dir) } #ifdef __MPI - MPI_Bcast(&step_rst_, 1, MPI_INT, 0, mdcell.communicator()); - MPI_Bcast(&md_tfirst, 1, MPI_DOUBLE, 0, mdcell.communicator()); - MPI_Bcast(&alpha, 1, MPI_DOUBLE, 0, mdcell.communicator()); - MPI_Bcast(&negative_count, 1, MPI_INT, 0, mdcell.communicator()); - MPI_Bcast(&dt_max, 1, MPI_DOUBLE, 0, mdcell.communicator()); - MPI_Bcast(&md_dt, 1, MPI_DOUBLE, 0, mdcell.communicator()); + MPI_Bcast( + &step_rst_, + 1, + MPI_INT, + 0, + mdcell.communicator()); + MPI_Bcast( + &md_tfirst, + 1, + MPI_DOUBLE, + 0, + mdcell.communicator()); + MPI_Bcast( + &alpha, + 1, + MPI_DOUBLE, + 0, + mdcell.communicator()); + MPI_Bcast( + &negative_count, + 1, + MPI_INT, + 0, + mdcell.communicator()); + MPI_Bcast( + &dt_max, + 1, + MPI_DOUBLE, + 0, + mdcell.communicator()); + MPI_Bcast( + &md_dt, + 1, + MPI_DOUBLE, + 0, + mdcell.communicator()); #endif return; @@ -172,7 +207,7 @@ void FIRE::check_force(void) std::int64_t movable_dof = 0; - for (const LocalAtom& atom : mdcell.owned_atoms()) + for (const LocalAtom &atom : mdcell.owned_atoms()) { for (int j = 0; j < 3; ++j) { @@ -197,13 +232,25 @@ void FIRE::check_force(void) } } - #ifdef __MPI +#ifdef __MPI if (mdcell.mpi_size() > 1) { double global_max = 0.0; std::int64_t global_movable_dof = 0; - MPI_Allreduce(&max, &global_max, 1, MPI_DOUBLE, MPI_MAX, mdcell.communicator()); - MPI_Allreduce(&movable_dof, &global_movable_dof, 1, MPI_INT64_T, MPI_SUM, mdcell.communicator()); + MPI_Allreduce( + &max, + &global_max, + 1, + MPI_DOUBLE, + MPI_MAX, + mdcell.communicator()); + MPI_Allreduce( + &movable_dof, + &global_movable_dof, + 1, + MPI_INT64_T, + MPI_SUM, + mdcell.communicator()); max = global_max; movable_dof = global_movable_dof; } @@ -224,7 +271,6 @@ void FIRE::check_force(void) return; } - void FIRE::check_fire(void) { double P = 0.0; @@ -242,7 +288,7 @@ void FIRE::check_fire(void) // Compute P, |F| and |v| only on movable degrees of freedom. // Fixed atoms/directions may have non-zero raw forces, but they should not // affect the FIRE velocity projection or adaptive time-step control. - for (LocalAtom& atom : mdcell.owned_atoms()) + for (LocalAtom &atom : mdcell.owned_atoms()) { for (int j = 0; j < 3; ++j) { @@ -261,14 +307,26 @@ void FIRE::check_fire(void) } } - #ifdef __MPI +#ifdef __MPI if (mdcell.mpi_size() > 1) { double local_values[3] = {P, sumforce, normvel}; double global_values[3] = {0.0, 0.0, 0.0}; std::int64_t global_movable_dof = 0; - MPI_Allreduce(local_values, global_values, 3, MPI_DOUBLE, MPI_SUM, mdcell.communicator()); - MPI_Allreduce(&movable_dof, &global_movable_dof, 1, MPI_INT64_T, MPI_SUM, mdcell.communicator()); + MPI_Allreduce( + local_values, + global_values, + 3, + MPI_DOUBLE, + MPI_SUM, + mdcell.communicator()); + MPI_Allreduce( + &movable_dof, + &global_movable_dof, + 1, + MPI_INT64_T, + MPI_SUM, + mdcell.communicator()); P = global_values[0]; sumforce = global_values[1]; normvel = global_values[2]; @@ -289,7 +347,7 @@ void FIRE::check_fire(void) // Avoid 0/0. In a truly converged case check_force() should stop the run. if (sumforce > 0.0 && normvel > 0.0) { - for (LocalAtom& atom : mdcell.owned_atoms()) + for (LocalAtom &atom : mdcell.owned_atoms()) { for (int j = 0; j < 3; ++j) { @@ -299,8 +357,7 @@ void FIRE::check_fire(void) continue; } - atom.vel[j] = (1.0 - alpha) * atom.vel[j] - + alpha * atom.force[j] / sumforce * normvel; + atom.vel[j] = (1.0 - alpha) * atom.vel[j] + alpha * atom.force[j] / sumforce * normvel; } } } @@ -319,7 +376,7 @@ void FIRE::check_fire(void) md_dt *= fdec; negative_count = 0; - for (LocalAtom& atom : mdcell.owned_atoms()) + for (LocalAtom &atom : mdcell.owned_atoms()) { for (int j = 0; j < 3; ++j) { diff --git a/source/source_md/fire.h b/source/source_md/fire.h index 5b84c3be9c..869b63dd1c 100644 --- a/source/source_md/fire.h +++ b/source/source_md/fire.h @@ -14,46 +14,37 @@ class FIRE : public MD_base { public: /// @param force_thr_in force convergence threshold used by the FIRE descent - FIRE(const MD_para& mdp_in, - const bool cal_stress_in, - const bool init_vel, - const int my_rank_in, - const double force_thr_in, - MDCell& mdcell_in); + FIRE( + const MD_para &mdp_in, + const bool cal_stress_in, + const bool init_vel, + const int my_rank_in, + const double force_thr_in, + MDCell &mdcell_in); ~FIRE(); /// @brief mixing coefficient of the FIRE velocity update - double get_alpha() const - { - return alpha; - } + double get_alpha() const { return alpha; } /// @brief largest time step the adaptive scheme is allowed to reach - double get_dt_max() const - { - return dt_max; - } + double get_dt_max() const { return dt_max; } /// @brief number of consecutive steps with negative power - int get_negative_count() const - { - return negative_count; - } + int get_negative_count() const { return negative_count; } private: + void setup(ModuleESolver::ESolver *p_esolver, const std::string &global_readin_dir, DomainDecomposition &decomp); - void setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp); - - void first_half(std::ofstream& ofs); + void first_half(std::ofstream &ofs); void second_half(); - void print_md(std::ofstream& ofs, const bool& cal_stress); + void print_md(std::ofstream &ofs, const bool &cal_stress); - void restart(const std::string& global_readin_dir); + void restart(const std::string &global_readin_dir); - void write_restart(const std::string& global_out_dir); + void write_restart(const std::string &global_out_dir); /** * @brief check the atomic forces converged or not diff --git a/source/source_md/langevin.cpp b/source/source_md/langevin.cpp index 204e946bbc..b8694cb9b5 100644 --- a/source/source_md/langevin.cpp +++ b/source/source_md/langevin.cpp @@ -1,26 +1,32 @@ #include "langevin.h" #include "md_func.h" -#include "source_cell/unitcell.h" #include "source_base/timer.h" +#include "source_cell/unitcell.h" -Langevin::Langevin(const MD_para& mdp_in, - const bool cal_stress_in, - const bool init_vel, - const int my_rank_in, - MDCell& mdcell_in) - : MD_base(mdp_in, cal_stress_in, init_vel, my_rank_in, mdcell_in) +Langevin::Langevin( + const MD_para &mdp_in, + const bool cal_stress_in, + const bool init_vel, + const int my_rank_in, + MDCell &mdcell_in) + : MD_base( + mdp_in, + cal_stress_in, + init_vel, + my_rank_in, + mdcell_in) { /// convert to a.u. unit - assert(ModuleBase::AU_to_FS!=0.0); + assert(ModuleBase::AU_to_FS != 0.0); md_damp = mdp.md_damp / ModuleBase::AU_to_FS; total_force.resize(static_cast(mdcell.owned_atoms().size())); } - -void Langevin::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp) +void Langevin::setup( + ModuleESolver::ESolver *p_esolver, const std::string &global_readin_dir, DomainDecomposition &decomp) { ModuleBase::TITLE("Langevin", "setup"); ModuleBase::timer::start("Langevin", "setup"); @@ -33,18 +39,20 @@ void Langevin::setup(ModuleESolver::ESolver* p_esolver, const std::string& globa return; } - -void Langevin::first_half(std::ofstream& ofs) +void Langevin::first_half(std::ofstream &ofs) { ModuleBase::TITLE("Langevin", "first_half"); ModuleBase::timer::start("Langevin", "first_half"); for (int i = 0; i < mdcell.owned_atoms().size(); ++i) { - LocalAtom& atom = mdcell.owned_atoms()[static_cast(i)]; + LocalAtom &atom = mdcell.owned_atoms()[static_cast(i)]; for (int k = 0; k < 3; ++k) { - if (atom.mbl[k]) atom.vel[k] += 0.5 * total_force[i][k] * md_dt / atom.mass; + if (atom.mbl[k]) + { + atom.vel[k] += 0.5 * total_force[i][k] * md_dt / atom.mass; + } } } MD_base::update_pos(); @@ -53,7 +61,6 @@ void Langevin::first_half(std::ofstream& ofs) return; } - void Langevin::second_half() { ModuleBase::TITLE("Langevin", "second_half"); @@ -62,10 +69,13 @@ void Langevin::second_half() post_force(); for (int i = 0; i < mdcell.owned_atoms().size(); ++i) { - LocalAtom& atom = mdcell.owned_atoms()[static_cast(i)]; + LocalAtom &atom = mdcell.owned_atoms()[static_cast(i)]; for (int k = 0; k < 3; ++k) { - if (atom.mbl[k]) atom.vel[k] += 0.5 * total_force[i][k] * md_dt / atom.mass; + if (atom.mbl[k]) + { + atom.vel[k] += 0.5 * total_force[i][k] * md_dt / atom.mass; + } } } @@ -73,31 +83,31 @@ void Langevin::second_half() return; } - -void Langevin::print_md(std::ofstream& ofs, const bool& cal_stress) +void Langevin::print_md(std::ofstream &ofs, const bool &cal_stress) { MD_base::print_md(ofs, cal_stress); return; } - -void Langevin::write_restart(const std::string& global_out_dir) +void Langevin::write_restart(const std::string &global_out_dir) { MD_base::write_restart(global_out_dir); return; } - -void Langevin::restart(const std::string& global_readin_dir) +void Langevin::restart(const std::string &global_readin_dir) { MD_base::restart(global_readin_dir); return; } - void Langevin::post_force() { - double t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); + double t_target = MD_func::target_temp( + step_ + step_rst_, + mdp.md_nstep, + md_tfirst, + md_tlast); total_force.resize(static_cast(mdcell.owned_atoms().size())); for (int i = 0; i < mdcell.owned_atoms().size(); ++i) @@ -107,12 +117,11 @@ void Langevin::post_force() { random_value[k] = static_cast(std::rand()) / RAND_MAX - 0.5; } - const LocalAtom& atom = mdcell.owned_atoms()[static_cast(i)]; + const LocalAtom &atom = mdcell.owned_atoms()[static_cast(i)]; ModuleBase::Vector3 fictitious_force = -atom.mass * atom.vel / md_damp; for (int j = 0; j < 3; ++j) { - fictitious_force[j] += sqrt(24.0 * t_target * atom.mass / md_damp / md_dt) - * random_value[j]; + fictitious_force[j] += sqrt(24.0 * t_target * atom.mass / md_damp / md_dt) * random_value[j]; } total_force[i] = atom.force + fictitious_force; } diff --git a/source/source_md/langevin.h b/source/source_md/langevin.h index 58ad4b63a6..dc7fa31480 100644 --- a/source/source_md/langevin.h +++ b/source/source_md/langevin.h @@ -16,24 +16,25 @@ class Langevin : public MD_base { public: - Langevin(const MD_para& mdp_in, - const bool cal_stress_in, - const bool init_vel, - const int my_rank_in, - MDCell& mdcell_in); + Langevin( + const MD_para &mdp_in, + const bool cal_stress_in, + const bool init_vel, + const int my_rank_in, + MDCell &mdcell_in); private: - void setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp); + void setup(ModuleESolver::ESolver *p_esolver, const std::string &global_readin_dir, DomainDecomposition &decomp); - void first_half(std::ofstream& ofs); + void first_half(std::ofstream &ofs); void second_half(); - void print_md(std::ofstream& ofs, const bool& cal_stress); + void print_md(std::ofstream &ofs, const bool &cal_stress); - void write_restart(const std::string& global_out_dir); + void write_restart(const std::string &global_out_dir); - void restart(const std::string& global_readin_dir); + void restart(const std::string &global_readin_dir); /** * @brief calculate fictitious forces @@ -41,8 +42,8 @@ class Langevin : public MD_base */ void post_force(); - std::vector > total_force; ///< total force = true force + Langevin fictitious_force - double md_damp; ///< damping factor + std::vector> total_force; ///< total force = true force + Langevin fictitious_force + double md_damp; ///< damping factor }; #endif diff --git a/source/source_md/md_base.cpp b/source/source_md/md_base.cpp index 6146e25117..cf08f8fb61 100644 --- a/source/source_md/md_base.cpp +++ b/source/source_md/md_base.cpp @@ -1,19 +1,22 @@ #include "md_base.h" + #include "md_func.h" #include "source_cell/unitcell.h" #ifdef __MPI #include "mpi.h" #endif #include "source_io/module_output/print_info.h" + #include #include -MD_base::MD_base(const MD_para& mdp_in, - const bool cal_stress_in, - const bool init_vel, - const int my_rank_in, - MDCell& mdcell_in) -: mdp(mdp_in), mdcell(mdcell_in) +MD_base::MD_base( + const MD_para &mdp_in, + const bool cal_stress_in, + const bool init_vel, + const int my_rank_in, + MDCell &mdcell_in) + : mdp(mdp_in), mdcell(mdcell_in) { #ifdef __MPI my_rank = mdcell.mpi_rank(); @@ -28,8 +31,8 @@ MD_base::MD_base(const MD_para& mdp_in, virial.create(3, 3); stress.create(3, 3); - assert(ModuleBase::AU_to_FS!=0.0); - assert(ModuleBase::Hartree_to_K!=0.0); + assert(ModuleBase::AU_to_FS != 0.0); + assert(ModuleBase::Hartree_to_K != 0.0); /// convert to a.u. unit md_dt = mdp.md_dt / ModuleBase::AU_to_FS; @@ -39,15 +42,19 @@ MD_base::MD_base(const MD_para& mdp_in, step_ = 0; step_rst_ = 0; - MD_func::init_vel(mdcell, init_vel, mdp.md_restart, md_tfirst, frozen_freedom_); + MD_func::init_vel( + mdcell, + init_vel, + mdp.md_restart, + md_tfirst, + frozen_freedom_); t_current = MD_func::current_temp(kinetic, mdcell, frozen_freedom_); } - MD_base::~MD_base() {} - -void MD_base::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp) +void MD_base::setup( + ModuleESolver::ESolver *p_esolver, const std::string &global_readin_dir, DomainDecomposition &decomp) { if (mdp.md_restart) { @@ -59,10 +66,22 @@ void MD_base::setup(ModuleESolver::ESolver* p_esolver, const std::string& global const int force_step = 0; const int istep_print = step_ + step_rst_ + 1; - ModuleIO::print_screen(stress_step, force_step, istep_print); - - MD_func::force_virial(p_esolver, step_, mdcell, decomp, potential, cal_stress, virial, mdp.md_out_force); - MD_func::compute_stress(mdcell, cal_stress, virial, stress); + ModuleIO::print_screen(stress_step, force_step, istep_print); + + MD_func::force_virial( + p_esolver, + step_, + mdcell, + decomp, + potential, + cal_stress, + virial, + mdp.md_out_force); + MD_func::compute_stress( + mdcell, + cal_stress, + virial, + stress); if (mdcell.has_backing_unitcell()) { mdcell.backing_unitcell().ionic_position_updated = true; @@ -71,8 +90,7 @@ void MD_base::setup(ModuleESolver::ESolver* p_esolver, const std::string& global return; } - -void MD_base::first_half(std::ofstream& ofs) +void MD_base::first_half(std::ofstream &ofs) { update_vel(); update_pos(); @@ -80,7 +98,6 @@ void MD_base::first_half(std::ofstream& ofs) return; } - void MD_base::second_half() { update_vel(); @@ -88,13 +105,12 @@ void MD_base::second_half() return; } - void MD_base::update_pos() { - std::vector& atoms = mdcell.owned_atoms(); + std::vector &atoms = mdcell.owned_atoms(); for (std::size_t i = 0; i < atoms.size(); ++i) { - LocalAtom& atom = atoms[i]; + LocalAtom &atom = atoms[i]; ModuleBase::Vector3 pos; for (int k = 0; k < 3; ++k) { @@ -124,13 +140,12 @@ void MD_base::update_pos() return; } - void MD_base::update_vel() { - std::vector& atoms = mdcell.owned_atoms(); + std::vector &atoms = mdcell.owned_atoms(); for (std::size_t i = 0; i < atoms.size(); ++i) { - LocalAtom& atom = atoms[i]; + LocalAtom &atom = atoms[i]; for (int k = 0; k < 3; ++k) { if (atom.mbl[k]) @@ -142,17 +157,16 @@ void MD_base::update_vel() return; } - -void MD_base::print_md(std::ofstream& ofs, const bool& cal_stress) +void MD_base::print_md(std::ofstream &ofs, const bool &cal_stress) { t_current = MD_func::current_temp(kinetic, mdcell, frozen_freedom_); - if (my_rank!=0) + if (my_rank != 0) { return; } - assert(ModuleBase::BOHR_RADIUS_SI>0.0); + assert(ModuleBase::BOHR_RADIUS_SI > 0.0); const double unit_transform = ModuleBase::HARTREE_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; double press = 0.0; @@ -196,7 +210,7 @@ void MD_base::print_md(std::ofstream& ofs, const bool& cal_stress) if (cal_stress) { MD_func::print_stress(ofs, virial, stress); - ofs << std::endl; + ofs << std::endl; } ofs << " ------------------------------------------------------------------------------------------------" @@ -228,8 +242,7 @@ void MD_base::print_md(std::ofstream& ofs, const bool& cal_stress) return; } - -void MD_base::write_restart(const std::string& global_out_dir) +void MD_base::write_restart(const std::string &global_out_dir) { if (!my_rank) { @@ -248,10 +261,13 @@ void MD_base::write_restart(const std::string& global_out_dir) return; } - -void MD_base::restart(const std::string& global_readin_dir) +void MD_base::restart(const std::string &global_readin_dir) { - MD_func::current_md_info(mdcell, global_readin_dir, step_rst_, md_tfirst); + MD_func::current_md_info( + mdcell, + global_readin_dir, + step_rst_, + md_tfirst); return; } diff --git a/source/source_md/md_base.h b/source/source_md/md_base.h index 4095c409a0..6bd7dafb3a 100644 --- a/source/source_md/md_base.h +++ b/source/source_md/md_base.h @@ -30,11 +30,12 @@ class MD_base * builds, where MDCell cannot supply it * @param mdcell_in mdcell information */ - MD_base(const MD_para& mdp_in, - const bool cal_stress_in, - const bool init_vel, - const int my_rank_in, - MDCell& mdcell_in); + MD_base( + const MD_para &mdp_in, + const bool cal_stress_in, + const bool init_vel, + const int my_rank_in, + MDCell &mdcell_in); virtual ~MD_base(); /** @@ -43,14 +44,15 @@ class MD_base * @param p_esolver the energy solver used in md * @param global_readin_dir directory of files for reading */ - virtual void setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp); + virtual void setup( + ModuleESolver::ESolver *p_esolver, const std::string &global_readin_dir, DomainDecomposition &decomp); /** * @brief the first half of equation of motion, update velocities and * positions * @param ofs determine the output files */ - virtual void first_half(std::ofstream& ofs); + virtual void first_half(std::ofstream &ofs); /** * @brief the second half of equation of motion, update velocities @@ -62,13 +64,13 @@ class MD_base * @param ofs determine the output files * @param cal_stress whether calculate and output stress */ - virtual void print_md(std::ofstream& ofs, const bool& cal_stress); + virtual void print_md(std::ofstream &ofs, const bool &cal_stress); /** * @brief write the information into files used for MD restarting * @param global_out_dir directory of output files */ - virtual void write_restart(const std::string& global_out_dir); + virtual void write_restart(const std::string &global_out_dir); /** * @brief restart MD when md_restart is true @@ -79,7 +81,7 @@ class MD_base * * @param global_readin_dir directory of files for reading */ - virtual void restart(const std::string& global_readin_dir); + virtual void restart(const std::string &global_readin_dir); protected: /** @@ -95,25 +97,22 @@ class MD_base public: /// @brief the time increment in a.u., converted from mdp.md_dt - double get_md_dt() const - { - return md_dt; - } - - bool stop; ///< MD stop or not - double t_current; ///< current temperature - int step_; ///< the MD step finished in current calculation - int step_rst_; ///< the MD step finished in previous calculations - std::int64_t frozen_freedom_; ///< the fixed freedom of the system - ModuleBase::matrix virial; ///< virial for this lattice - ModuleBase::matrix stress; ///< stress for this lattice - double potential=0.0; ///< potential energy - double kinetic; ///< kinetic energy + double get_md_dt() const { return md_dt; } + + bool stop; ///< MD stop or not + double t_current; ///< current temperature + int step_; ///< the MD step finished in current calculation + int step_rst_; ///< the MD step finished in previous calculations + std::int64_t frozen_freedom_; ///< the fixed freedom of the system + ModuleBase::matrix virial; ///< virial for this lattice + ModuleBase::matrix stress; ///< stress for this lattice + double potential = 0.0; ///< potential energy + double kinetic; ///< kinetic energy protected: - const MD_para& mdp; ///< input parameters used in md - MDCell& mdcell; ///< mdcell information - double energy_=0.0; ///< total energy of the system + const MD_para &mdp; ///< input parameters used in md + MDCell &mdcell; ///< mdcell information + double energy_ = 0.0; ///< total energy of the system bool cal_stress; ///< whether calculate stress int my_rank; ///< MPI rank of the processor diff --git a/source/source_md/md_func.cpp b/source/source_md/md_func.cpp index 3536dd4a0b..e231b24dc7 100644 --- a/source/source_md/md_func.cpp +++ b/source/source_md/md_func.cpp @@ -1,9 +1,9 @@ #include "md_func.h" -#include "source_cell/domain_decomposition.h" #include "source_base/global_variable.h" #include "source_base/parallel_cell.h" #include "source_base/timer.h" +#include "source_cell/domain_decomposition.h" #include "source_io/module_output/output_log.h" #include @@ -14,17 +14,16 @@ #include #include - namespace MD_func { double gaussrand() { - static double v1=0.0; - static double v2=0.0; - static double S=0.0; + static double v1 = 0.0; + static double v2 = 0.0; + static double S = 0.0; static int phase = 0; - double xx=0.0; + double xx = 0.0; if (phase == 0) { @@ -50,11 +49,11 @@ double gaussrand() return xx; } -double kinetic_energy(const int& natom, const ModuleBase::Vector3* vel, const double* allmass) +double kinetic_energy(const int &natom, const ModuleBase::Vector3 *vel, const double *allmass) { double ke = 0; -#pragma omp parallel for reduction(+:ke) schedule(static) if (natom >= 256) +#pragma omp parallel for reduction(+ : ke) schedule(static) if (natom >= 256) for (int ion = 0; ion < natom; ++ion) { ke += 0.5 * allmass[ion] * vel[ion].norm2(); @@ -63,27 +62,47 @@ double kinetic_energy(const int& natom, const ModuleBase::Vector3* vel, return ke; } -void compute_stress(const MDCell& mdcell, - const bool& cal_stress, - const ModuleBase::matrix& virial, - ModuleBase::matrix& stress) +void compute_stress( + const MDCell &mdcell, + const bool &cal_stress, + const ModuleBase::matrix &virial, + ModuleBase::matrix &stress) { - if (!cal_stress) return; + if (!cal_stress) + { + return; + } ModuleBase::matrix t_vector(3, 3); for (std::size_t i = 0; i < mdcell.owned_atoms().size(); ++i) { - const LocalAtom& atom = mdcell.owned_atoms()[i]; - for (int a = 0; a < 3; ++a) for (int b = 0; b < 3; ++b) - t_vector(a, b) += atom.mass * atom.vel[a] * atom.vel[b]; + const LocalAtom &atom = mdcell.owned_atoms()[i]; + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + t_vector(a, b) += atom.mass * atom.vel[a] * atom.vel[b]; + } + } } #ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, t_vector.c, 9, MPI_DOUBLE, MPI_SUM, mdcell.communicator()); + MPI_Allreduce( + MPI_IN_PLACE, + t_vector.c, + 9, + MPI_DOUBLE, + MPI_SUM, + mdcell.communicator()); #endif - for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) - stress(i, j) = virial(i, j) + t_vector(i, j) / mdcell.omega(); + for (int i = 0; i < 3; ++i) + { + for (int j = 0; j < 3; ++j) + { + stress(i, j) = virial(i, j) + t_vector(i, j) / mdcell.omega(); + } + } } -void read_vel(const UnitCell& unit_in, ModuleBase::Vector3* vel) +void read_vel(const UnitCell &unit_in, ModuleBase::Vector3 *vel) { int iat = 0; for (int it = 0; it < unit_in.ntype; ++it) @@ -111,11 +130,12 @@ void read_vel(const UnitCell& unit_in, ModuleBase::Vector3* vel) return; } -void rescale_vel(const int& natom, - const double& temperature, - const double* allmass, - const int& frozen_freedom, - ModuleBase::Vector3* vel) +void rescale_vel( + const int &natom, + const double &temperature, + const double *allmass, + const int &frozen_freedom, + ModuleBase::Vector3 *vel) { double factor = 0.0; if (3 * natom == frozen_freedom || temperature == 0) @@ -134,14 +154,15 @@ void rescale_vel(const int& natom, } } -void rand_vel(const int& natom, - const double& temperature, - const double* allmass, - const int& frozen_freedom, - const ModuleBase::Vector3 frozen, - const ModuleBase::Vector3* ionmbl, - const int& my_rank, - ModuleBase::Vector3* vel) +void rand_vel( + const int &natom, + const double &temperature, + const double *allmass, + const int &frozen_freedom, + const ModuleBase::Vector3 frozen, + const ModuleBase::Vector3 *ionmbl, + const int &my_rank, + ModuleBase::Vector3 *vel) { if (!my_rank) { @@ -181,43 +202,77 @@ void rand_vel(const int& natom, } // rescale the velocity to the target temperature - rescale_vel(natom, temperature, allmass, frozen_freedom, vel); + rescale_vel( + natom, + temperature, + allmass, + frozen_freedom, + vel); } #ifdef __MPI - MPI_Bcast(vel, natom * 3, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast( + vel, + natom * 3, + MPI_DOUBLE, + 0, + MPI_COMM_WORLD); #endif return; } -void init_vel(MDCell& mdcell, - const bool& init_vel, - const bool& restart, - double& temperature, - std::int64_t& frozen_freedom) +void init_vel( + MDCell &mdcell, + const bool &init_vel, + const bool &restart, + double &temperature, + std::int64_t &frozen_freedom) { if (mdcell.has_backing_unitcell()) { - const UnitCell& unit = mdcell.backing_unitcell(); + const UnitCell &unit = mdcell.backing_unitcell(); std::vector masses(static_cast(unit.nat), 0.0); - std::vector > mobile(static_cast(unit.nat)); - std::vector > velocities(static_cast(unit.nat)); + std::vector> mobile(static_cast(unit.nat)); + std::vector> velocities(static_cast(unit.nat)); ModuleBase::Vector3 frozen; - get_mass_mbl(unit, masses.data(), frozen, mobile.data()); + get_mass_mbl( + unit, + masses.data(), + frozen, + mobile.data()); frozen_freedom = frozen.x + frozen.y + frozen.z; - if (frozen.x == 0) ++frozen_freedom; - if (frozen.y == 0) ++frozen_freedom; - if (frozen.z == 0) ++frozen_freedom; + if (frozen.x == 0) + { + ++frozen_freedom; + } + if (frozen.y == 0) + { + ++frozen_freedom; + } + if (frozen.z == 0) + { + ++frozen_freedom; + } if (init_vel) { read_vel(unit, velocities.data()); double kinetic = 0.0; - const double current = current_temp(kinetic, unit.nat, frozen_freedom, masses.data(), velocities.data()); + const double current = current_temp( + kinetic, + unit.nat, + frozen_freedom, + masses.data(), + velocities.data()); if (!restart && temperature >= 0.0 && current > 0.0) { - rescale_vel(unit.nat, temperature, masses.data(), frozen_freedom, velocities.data()); + rescale_vel( + unit.nat, + temperature, + masses.data(), + frozen_freedom, + velocities.data()); } else if (!restart && temperature < 0.0) { @@ -231,14 +286,15 @@ void init_vel(MDCell& mdcell, #else const int rank = 0; #endif - rand_vel(unit.nat, - temperature, - masses.data(), - static_cast(frozen_freedom), - frozen, - mobile.data(), - rank, - velocities.data()); + rand_vel( + unit.nat, + temperature, + masses.data(), + static_cast(frozen_freedom), + frozen, + mobile.data(), + rank, + velocities.data()); } std::vector offsets(static_cast(unit.ntype + 1), 0); @@ -246,27 +302,51 @@ void init_vel(MDCell& mdcell, { offsets[static_cast(it + 1)] = offsets[static_cast(it)] + unit.atoms[it].na; } - for (LocalAtom& atom : mdcell.owned_atoms()) + for (LocalAtom &atom : mdcell.owned_atoms()) { - atom.vel = velocities[static_cast(offsets[static_cast(atom.type)] + atom.type_index)]; + atom.vel = + velocities[static_cast(offsets[static_cast(atom.type)] + atom.type_index)]; } return; } - std::vector& atoms = mdcell.owned_atoms(); + std::vector &atoms = mdcell.owned_atoms(); ModuleBase::Vector3 frozen(0, 0, 0); for (std::size_t i = 0; i < atoms.size(); ++i) - for (int k = 0; k < 3; ++k) if (!atoms[i].mbl[k]) ++frozen[k]; + { + for (int k = 0; k < 3; ++k) + { + if (!atoms[i].mbl[k]) + { + ++frozen[k]; + } + } + } #ifdef __MPI if (mdcell.mpi_size() > 1) { - MPI_Allreduce(MPI_IN_PLACE, &frozen.x, 3, MPI_INT64_T, MPI_SUM, mdcell.communicator()); + MPI_Allreduce( + MPI_IN_PLACE, + &frozen.x, + 3, + MPI_INT64_T, + MPI_SUM, + mdcell.communicator()); } #endif frozen_freedom = frozen.x + frozen.y + frozen.z; - if (!frozen.x) ++frozen_freedom; - if (!frozen.y) ++frozen_freedom; - if (!frozen.z) ++frozen_freedom; + if (!frozen.x) + { + ++frozen_freedom; + } + if (!frozen.y) + { + ++frozen_freedom; + } + if (!frozen.z) + { + ++frozen_freedom; + } if (init_vel) { double kinetic = 0.0; @@ -274,31 +354,52 @@ void init_vel(MDCell& mdcell, if (!restart && current > 0.0 && temperature > 0.0) { const double factor = sqrt(temperature / current); - for (LocalAtom& atom : atoms) atom.vel *= factor; + for (LocalAtom &atom : atoms) + { + atom.vel *= factor; + } } return; } double local_mass = 0.0; ModuleBase::Vector3 momentum(0.0, 0.0, 0.0); - for (LocalAtom& atom : atoms) + for (LocalAtom &atom : atoms) { local_mass += atom.mass; for (int k = 0; k < 3; ++k) { atom.vel[k] = atom.mbl[k] ? gaussrand() * sqrt(temperature / atom.mass) : 0.0; - if (frozen[k] == 0) momentum[k] += atom.mass * atom.vel[k]; + if (frozen[k] == 0) + { + momentum[k] += atom.mass * atom.vel[k]; + } } } #ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, &local_mass, 1, MPI_DOUBLE, MPI_SUM, mdcell.communicator()); - MPI_Allreduce(MPI_IN_PLACE, &momentum.x, 3, MPI_DOUBLE, MPI_SUM, mdcell.communicator()); + MPI_Allreduce( + MPI_IN_PLACE, + &local_mass, + 1, + MPI_DOUBLE, + MPI_SUM, + mdcell.communicator()); + MPI_Allreduce( + MPI_IN_PLACE, + &momentum.x, + 3, + MPI_DOUBLE, + MPI_SUM, + mdcell.communicator()); #endif for (int k = 0; k < 3; ++k) { if (frozen[k] == 0 && local_mass > 0.0) { - for (LocalAtom& atom : atoms) atom.vel[k] -= momentum[k] / local_mass; + for (LocalAtom &atom : atoms) + { + atom.vel[k] -= momentum[k] / local_mass; + } } } double kinetic = 0.0; @@ -306,39 +407,49 @@ void init_vel(MDCell& mdcell, if (current > 0.0 && temperature > 0.0) { const double factor = sqrt(temperature / current); - for (LocalAtom& atom : atoms) atom.vel *= factor; + for (LocalAtom &atom : atoms) + { + atom.vel *= factor; + } } } -void force_virial(ModuleESolver::ESolver* p_esolver, - const int& istep, - MDCell& mdcell, - DomainDecomposition& decomp, - double& potential, - const bool& cal_stress, - ModuleBase::matrix& virial, - const bool& md_out_force) +void force_virial( + ModuleESolver::ESolver *p_esolver, + const int &istep, + MDCell &mdcell, + DomainDecomposition &decomp, + double &potential, + const bool &cal_stress, + ModuleBase::matrix &virial, + const bool &md_out_force) { ModuleBase::TITLE("MD_func", "force_virial"); ModuleBase::timer::start("MD_func", "force_virial"); if (!mdcell.has_backing_unitcell()) { decomp.prepare_neighbors(mdcell); - p_esolver->runner(static_cast(mdcell), istep); + p_esolver->runner(static_cast(mdcell), istep); decomp.accumulate_ghost_forces(mdcell); potential = 0.5 * p_esolver->cal_energy(); - for (LocalAtom& atom : mdcell.owned_atoms()) atom.force *= 0.5; + for (LocalAtom &atom : mdcell.owned_atoms()) + { + atom.force *= 0.5; + } if (md_out_force) { ModuleIO::print_force(GlobalV::ofs_running, mdcell, "TOTAL-FORCE (eV/Angstrom)"); } - if (cal_stress) { p_esolver->cal_stress(static_cast(mdcell), virial); virial *= 0.5; } + if (cal_stress) + { + p_esolver->cal_stress(static_cast(mdcell), virial); + virial *= 0.5; + } } else { - UnitCell& ucell = mdcell.backing_unitcell(); - std::vector>> backing_velocities( - static_cast(ucell.ntype)); + UnitCell &ucell = mdcell.backing_unitcell(); + std::vector>> backing_velocities(static_cast(ucell.ntype)); for (int it = 0; it < ucell.ntype; ++it) { backing_velocities[static_cast(it)] = ucell.atoms[it].vel; @@ -346,21 +457,35 @@ void force_virial(ModuleESolver::ESolver* p_esolver, mdcell.sync_backing_unitcell(); for (int it = 0; it < ucell.ntype; ++it) { - std::copy(backing_velocities[static_cast(it)].begin(), - backing_velocities[static_cast(it)].end(), - ucell.atoms[it].vel.begin()); + std::copy( + backing_velocities[static_cast(it)].begin(), + backing_velocities[static_cast(it)].end(), ucell.atoms[it].vel.begin()); + } + p_esolver->runner(ucell, istep); + potential = 0.5 * p_esolver->cal_energy(); + ModuleBase::matrix full_force(ucell.nat, 3); + p_esolver->cal_force(ucell, full_force); + full_force *= 0.5; + if (cal_stress) + { + p_esolver->cal_stress(ucell, virial); + virial *= 0.5; + } + std::vector offsets(ucell.ntype + 1, 0); + for (int it = 0; it < ucell.ntype; ++it) + { + offsets[it + 1] = offsets[it] + ucell.atoms[it].na; + } + for (LocalAtom &atom : mdcell.owned_atoms()) + { + const int iat = offsets[atom.type] + atom.type_index; + atom.force.set(full_force(iat, 0), full_force(iat, 1), full_force(iat, 2)); } - p_esolver->runner(ucell, istep); potential = 0.5 * p_esolver->cal_energy(); - ModuleBase::matrix full_force(ucell.nat, 3); p_esolver->cal_force(ucell, full_force); full_force *= 0.5; - if (cal_stress) { p_esolver->cal_stress(ucell, virial); virial *= 0.5; } - std::vector offsets(ucell.ntype + 1, 0); for (int it=0; it rows; @@ -424,31 +550,40 @@ void write_dump_rows(std::ostream& output, { ModuleBase::WARNING_QUIT("MD_func::dump_info", "missing MD_dump atoms."); } - for (const auto& row : rows) + for (const auto &row : rows) { output << row.second << "\n"; } } } // namespace -void dump_info(const int& step, - const std::string& global_out_dir, - const MDCell& mdcell, - const MD_para& mdp, - const bool cal_stress, - const ModuleBase::matrix& virial) +void dump_info( + const int &step, + const std::string &global_out_dir, + const MDCell &mdcell, + const MD_para &mdp, + const bool cal_stress, + const ModuleBase::matrix &virial) { const int max_batch_bytes = 1024 * 1024; - dump_info(step, global_out_dir, mdcell, mdp, cal_stress, virial, max_batch_bytes); + dump_info( + step, + global_out_dir, + mdcell, + mdp, + cal_stress, + virial, + max_batch_bytes); } -void dump_info(const int& step, - const std::string& global_out_dir, - const MDCell& mdcell, - const MD_para& mdp, - const bool cal_stress, - const ModuleBase::matrix& virial, - const int max_batch_bytes) +void dump_info( + const int &step, + const std::string &global_out_dir, + const MDCell &mdcell, + const MD_para &mdp, + const bool cal_stress, + const ModuleBase::matrix &virial, + const int max_batch_bytes) { if (max_batch_bytes <= 0) { @@ -470,41 +605,55 @@ void dump_info(const int& step, { header << "VIRIAL (kbar)\n"; for (int i = 0; i < 3; ++i) - header << " " << virial(i, 0) * unit_virial << " " << virial(i, 1) * unit_virial << " " << virial(i, 2) * unit_virial << "\n"; + { + header << " " << virial(i, 0) * unit_virial << " " << virial(i, 1) * unit_virial << " " + << virial(i, 2) * unit_virial << "\n"; + } } header << "INDEX LABEL POSITION (Angstrom)"; - if (mdp.dump_force) header << " FORCE (eV/Angstrom)"; - if (mdp.dump_vel) header << " VELOCITY (Angstrom/fs)"; + if (mdp.dump_force) + { + header << " FORCE (eV/Angstrom)"; + } + if (mdp.dump_vel) + { + header << " VELOCITY (Angstrom/fs)"; + } header << "\n"; std::vector type_offsets(mdcell.type_atom_counts().size() + 1, 0); for (std::size_t it = 0; it < mdcell.type_atom_counts().size(); ++it) { type_offsets[it + 1] = type_offsets[it] + mdcell.type_atom_counts()[it]; } - const auto format_atom = [&](const LocalAtom& atom) -> std::string + const auto format_atom = [&](const LocalAtom &atom) -> std::string { std::ostringstream row; row << std::fixed << std::setprecision(12); - row << " " << type_offsets[static_cast(atom.type)] + atom.type_index - << " " << mdcell.type_labels()[static_cast(atom.type)] - << " " << atom.cart.x * unit_pos << " " << atom.cart.y * unit_pos << " " << atom.cart.z * unit_pos; + row << " " << type_offsets[static_cast(atom.type)] + atom.type_index << " " + << mdcell.type_labels()[static_cast(atom.type)] << " " << atom.cart.x * unit_pos << " " + << atom.cart.y * unit_pos << " " << atom.cart.z * unit_pos; if (mdp.dump_force) - row << " " << atom.force.x * unit_force << " " << atom.force.y * unit_force << " " << atom.force.z * unit_force; + { + row << " " << atom.force.x * unit_force << " " << atom.force.y * unit_force << " " + << atom.force.z * unit_force; + } if (mdp.dump_vel) + { row << " " << atom.vel.x * unit_vel << " " << atom.vel.y * unit_vel << " " << atom.vel.z * unit_vel; + } row << "\n"; return row.str(); }; - std::vector> ordered_atoms; + std::vector> ordered_atoms; const std::size_t local_atom_count = mdcell.owned_atoms().size(); ordered_atoms.reserve(local_atom_count); int local_row_bytes = 0; int local_overflow = 0; const std::size_t maximum_count = std::numeric_limits::max(); - for (const LocalAtom& atom : mdcell.owned_atoms()) + for (const LocalAtom &atom : mdcell.owned_atoms()) { const std::int64_t index = type_offsets[static_cast(atom.type)] + atom.type_index; - const LocalAtom* atom_pointer = &atom; + const LocalAtom *atom_pointer = &atom; ordered_atoms.emplace_back(index, atom_pointer); const std::string row = format_atom(atom); if (row.size() > maximum_count) @@ -517,10 +666,10 @@ void dump_info(const int& step, local_row_bytes = std::max(local_row_bytes, row_bytes); } } - std::sort(ordered_atoms.begin(), ordered_atoms.end(), - [](const std::pair& left, - const std::pair& right) - { return left.first < right.first; }); + std::sort( + ordered_atoms.begin(), ordered_atoms.end(), + [](const std::pair &left, + const std::pair &right) { return left.first < right.first; }); ModuleBase::CommunicationDomain domain; #ifdef __MPI domain.initialize(mdcell.communicator()); @@ -556,7 +705,7 @@ void dump_info(const int& step, std::ostringstream local; while (local_index < ordered_atoms.size() && ordered_atoms[local_index].first < end) { - const LocalAtom& atom = *ordered_atoms[local_index].second; + const LocalAtom &atom = *ordered_atoms[local_index].second; local << format_atom(atom); ++local_index; } @@ -564,7 +713,11 @@ void dump_info(const int& step, const std::string gathered = domain.gather_text(local_text, batch_bytes); if (domain.rank() == 0) { - write_dump_rows(ofs, gathered, begin, end); + write_dump_rows( + ofs, + gathered, + begin, + end); } begin = end; } @@ -579,10 +732,11 @@ void dump_info(const int& step, } } -void get_mass_mbl(const UnitCell& unit_in, - double* allmass, - ModuleBase::Vector3& frozen, - ModuleBase::Vector3* ionmbl) +void get_mass_mbl( + const UnitCell &unit_in, + double *allmass, + ModuleBase::Vector3 &frozen, + ModuleBase::Vector3 *ionmbl) { int ion = 0; frozen.set(0, 0, 0); @@ -592,15 +746,18 @@ void get_mass_mbl(const UnitCell& unit_in, { allmass[ion] = unit_in.atoms[it].mass / ModuleBase::AU_to_MASS; ionmbl[ion] = unit_in.atoms[it].mbl[i]; - if (ionmbl[ion].x == 0) { + if (ionmbl[ion].x == 0) + { ++frozen.x; -} - if (ionmbl[ion].y == 0) { + } + if (ionmbl[ion].y == 0) + { ++frozen.y; -} - if (ionmbl[ion].z == 0) { + } + if (ionmbl[ion].z == 0) + { ++frozen.z; -} + } ion++; } @@ -609,18 +766,23 @@ void get_mass_mbl(const UnitCell& unit_in, return; } -double target_temp(const int& istep, const int& nstep, const double& tfirst, const double& tlast) +double target_temp( + const int &istep, + const int &nstep, + const double &tfirst, + const double &tlast) { - assert(nstep>0); + assert(nstep > 0); double delta = static_cast(istep) / nstep; return tfirst + delta * (tlast - tfirst); } -double current_temp(double& kinetic, - const int& natom, - const int& frozen_freedom, - const double* allmass, - const ModuleBase::Vector3* vel) +double current_temp( + double &kinetic, + const int &natom, + const int &frozen_freedom, + const double *allmass, + const ModuleBase::Vector3 *vel) { if (3 * natom == frozen_freedom) { @@ -631,46 +793,81 @@ double current_temp(double& kinetic, return 2 * kinetic / (3 * natom - frozen_freedom); } -double current_temp(double& kinetic, - const MDCell& mdcell, - const std::int64_t& frozen_freedom) +double current_temp(double &kinetic, const MDCell &mdcell, const std::int64_t &frozen_freedom) { kinetic = 0.0; for (std::size_t i = 0; i < mdcell.owned_atoms().size(); ++i) { - const LocalAtom& atom = mdcell.owned_atoms()[i]; + const LocalAtom &atom = mdcell.owned_atoms()[i]; kinetic += 0.5 * atom.mass * atom.vel.norm2(); } #ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_DOUBLE, MPI_SUM, mdcell.communicator()); + MPI_Allreduce( + MPI_IN_PLACE, + &kinetic, + 1, + MPI_DOUBLE, + MPI_SUM, + mdcell.communicator()); #endif const std::int64_t dof = 3 * mdcell.nat() - frozen_freedom; - if (dof == 0) return 0.0; + if (dof == 0) + { + return 0.0; + } return 2.0 * kinetic / static_cast(dof); } -std::int64_t global_dof(const MDCell& mdcell) +std::int64_t global_dof(const MDCell &mdcell) { std::int64_t local_frozen[3] = {0, 0, 0}; for (int i = 0; i < mdcell.owned_atoms().size(); ++i) { - const ModuleBase::Vector3& mbl = mdcell.owned_atoms()[static_cast(i)].mbl; - if (mbl.x == 0) ++local_frozen[0]; - if (mbl.y == 0) ++local_frozen[1]; - if (mbl.z == 0) ++local_frozen[2]; + const ModuleBase::Vector3 &mbl = mdcell.owned_atoms()[static_cast(i)].mbl; + if (mbl.x == 0) + { + ++local_frozen[0]; + } + if (mbl.y == 0) + { + ++local_frozen[1]; + } + if (mbl.z == 0) + { + ++local_frozen[2]; + } } std::int64_t global_frozen[3] = {local_frozen[0], local_frozen[1], local_frozen[2]}; #ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, global_frozen, 3, MPI_INT64_T, MPI_SUM, mdcell.communicator()); + MPI_Allreduce( + MPI_IN_PLACE, + global_frozen, + 3, + MPI_INT64_T, + MPI_SUM, + mdcell.communicator()); #endif std::int64_t total_frozen = global_frozen[0] + global_frozen[1] + global_frozen[2]; - if (global_frozen[0] == 0) ++total_frozen; - if (global_frozen[1] == 0) ++total_frozen; - if (global_frozen[2] == 0) ++total_frozen; + if (global_frozen[0] == 0) + { + ++total_frozen; + } + if (global_frozen[1] == 0) + { + ++total_frozen; + } + if (global_frozen[2] == 0) + { + ++total_frozen; + } return 3 * mdcell.nat() - total_frozen; } -void current_md_info(const MDCell& mdcell, const std::string& file_dir, int& md_step, double& temperature) +void current_md_info( + const MDCell &mdcell, + const std::string &file_dir, + int &md_step, + double &temperature) { bool ok = true; @@ -695,15 +892,30 @@ void current_md_info(const MDCell& mdcell, const std::string& file_dir, int& md_ } #ifdef __MPI - MPI_Bcast(&ok, 1, MPI_C_BOOL, 0, mdcell.communicator()); + MPI_Bcast( + &ok, + 1, + MPI_C_BOOL, + 0, + mdcell.communicator()); #endif if (!ok) { ModuleBase::WARNING_QUIT("current_md_info", "no Restart_md.txt!"); } #ifdef __MPI - MPI_Bcast(&md_step, 1, MPI_INT, 0, mdcell.communicator()); - MPI_Bcast(&temperature, 1, MPI_DOUBLE, 0, mdcell.communicator()); + MPI_Bcast( + &md_step, + 1, + MPI_INT, + 0, + mdcell.communicator()); + MPI_Bcast( + &temperature, + 1, + MPI_DOUBLE, + 0, + mdcell.communicator()); #endif } diff --git a/source/source_md/md_func.h b/source/source_md/md_func.h index 5dc21f2f85..f07cacbf5b 100644 --- a/source/source_md/md_func.h +++ b/source/source_md/md_func.h @@ -29,11 +29,12 @@ namespace MD_func */ double gaussrand(); -void init_vel(MDCell& mdcell, - const bool& init_vel, - const bool& restart, - double& temperature, - std::int64_t& frozen_freedom); +void init_vel( + MDCell &mdcell, + const bool &init_vel, + const bool &restart, + double &temperature, + std::int64_t &frozen_freedom); /** * @brief read in atomic velocities from STRU @@ -41,7 +42,7 @@ void init_vel(MDCell& mdcell, * @param unit_in unitcell information * @param vel the read-in atomic velocities */ -void read_vel(const UnitCell& unit_in, ModuleBase::Vector3* vel); +void read_vel(const UnitCell &unit_in, ModuleBase::Vector3 *vel); /** * @brief generate atomic velocities that satisfy the Boltzmann distribution @@ -55,14 +56,15 @@ void read_vel(const UnitCell& unit_in, ModuleBase::Vector3* vel); * @param my_rank MPI rank of the processor * @param vel the genarated atomic velocities */ -void rand_vel(const int& natom, - const double& temperature, - const double* allmass, - const int& frozen_freedom, - const ModuleBase::Vector3 frozen, - const ModuleBase::Vector3* ionmbl, - const int& my_rank, - ModuleBase::Vector3* vel); +void rand_vel( + const int &natom, + const double &temperature, + const double *allmass, + const int &frozen_freedom, + const ModuleBase::Vector3 frozen, + const ModuleBase::Vector3 *ionmbl, + const int &my_rank, + ModuleBase::Vector3 *vel); /** * @brief rescale the velocity to the target temperature @@ -73,11 +75,12 @@ void rand_vel(const int& natom, * @param frozen_freedom the fixed freedom * @param vel the genarated atomic velocities */ -void rescale_vel(const int& natom, - const double& temperature, - const double* allmass, - const int& frozen_freedom, - ModuleBase::Vector3* vel); +void rescale_vel( + const int &natom, + const double &temperature, + const double *allmass, + const int &frozen_freedom, + ModuleBase::Vector3 *vel); /** * @brief calculate energy, forces and virial tensor @@ -89,14 +92,15 @@ void rescale_vel(const int& natom, * @param cal_stress whether calculate stress * @param virial lattice virial tensor */ -void force_virial(ModuleESolver::ESolver* p_esolver, - const int& istep, - MDCell& mdcell, - DomainDecomposition& decomp, - double& potential, - const bool& cal_stress, - ModuleBase::matrix& virial, - const bool& md_out_force); +void force_virial( + ModuleESolver::ESolver *p_esolver, + const int &istep, + MDCell &mdcell, + DomainDecomposition &decomp, + double &potential, + const bool &cal_stress, + ModuleBase::matrix &virial, + const bool &md_out_force); /** * @brief calculate the ionic kinetic energy * @@ -105,7 +109,7 @@ void force_virial(ModuleESolver::ESolver* p_esolver, * @param allmass atomic mass * @return the ionic kinetic energy */ -double kinetic_energy(const int& natom, const ModuleBase::Vector3* vel, const double* allmass); +double kinetic_energy(const int &natom, const ModuleBase::Vector3 *vel, const double *allmass); /** * @brief calculate the total stress tensor @@ -115,10 +119,11 @@ double kinetic_energy(const int& natom, const ModuleBase::Vector3* vel, * @param virial lattice virial tensor * @param stress total stress tensor */ -void compute_stress(const MDCell& mdcell, - const bool& cal_stress, - const ModuleBase::matrix& virial, - ModuleBase::matrix& stress); +void compute_stress( + const MDCell &mdcell, + const bool &cal_stress, + const ModuleBase::matrix &virial, + ModuleBase::matrix &stress); /** * @brief output the stress information @@ -127,7 +132,7 @@ void compute_stress(const MDCell& mdcell, * @param virial lattice virial tensor * @param stress total stress tensor */ -void print_stress(std::ofstream& ofs, const ModuleBase::matrix& virial, const ModuleBase::matrix& stress); +void print_stress(std::ofstream &ofs, const ModuleBase::matrix &virial, const ModuleBase::matrix &stress); /** * @brief dump the md information @@ -142,21 +147,23 @@ void print_stress(std::ofstream& ofs, const ModuleBase::matrix& virial, const Mo * @param cal_stress whether stress is calculated * @param virial lattice virial tensor */ -void dump_info(const int& step, - const std::string& global_out_dir, - const MDCell& mdcell, - const MD_para& mdp, - const bool cal_stress, - const ModuleBase::matrix& virial); +void dump_info( + const int &step, + const std::string &global_out_dir, + const MDCell &mdcell, + const MD_para &mdp, + const bool cal_stress, + const ModuleBase::matrix &virial); /// Explicit buffer budget for bounded dump collectives; one row is the minimum batch. -void dump_info(const int& step, - const std::string& global_out_dir, - const MDCell& mdcell, - const MD_para& mdp, - const bool cal_stress, - const ModuleBase::matrix& virial, - const int max_batch_bytes); +void dump_info( + const int &step, + const std::string &global_out_dir, + const MDCell &mdcell, + const MD_para &mdp, + const bool cal_stress, + const ModuleBase::matrix &virial, + const int max_batch_bytes); /** * @brief obtain the atomic mass and whether the freedom is fixed @@ -166,10 +173,11 @@ void dump_info(const int& step, * @param frozen the fixed freedom along three directions * @param ionmbl determine whether the atomic freedom is fixed */ -void get_mass_mbl(const UnitCell& unit_in, - double* allmass, - ModuleBase::Vector3& frozen, - ModuleBase::Vector3* ionmbl); +void get_mass_mbl( + const UnitCell &unit_in, + double *allmass, + ModuleBase::Vector3 &frozen, + ModuleBase::Vector3 *ionmbl); /** * @brief get the target temperature of the current md step @@ -180,7 +188,11 @@ void get_mass_mbl(const UnitCell& unit_in, * @param tlast the final temperature * @return the target temperature */ -double target_temp(const int& istep, const int& nstep, const double& tfirst, const double& tlast); +double target_temp( + const int &istep, + const int &nstep, + const double &tfirst, + const double &tlast); /** * @brief get the current temperature @@ -192,17 +204,20 @@ double target_temp(const int& istep, const int& nstep, const double& tfirst, con * @param vel atomic velocities * @return the current temperature */ -double current_temp(double& kinetic, - const int& natom, - const int& frozen_freedom, - const double* allmass, - const ModuleBase::Vector3* vel); -double current_temp(double& kinetic, - const MDCell& mdcell, - const std::int64_t& frozen_freedom); -std::int64_t global_dof(const MDCell& mdcell); - -void current_md_info(const MDCell& mdcell, const std::string& file_dir, int& md_step, double& temperature); +double current_temp( + double &kinetic, + const int &natom, + const int &frozen_freedom, + const double *allmass, + const ModuleBase::Vector3 *vel); +double current_temp(double &kinetic, const MDCell &mdcell, const std::int64_t &frozen_freedom); +std::int64_t global_dof(const MDCell &mdcell); + +void current_md_info( + const MDCell &mdcell, + const std::string &file_dir, + int &md_step, + double &temperature); } // namespace MD_func diff --git a/source/source_md/module_plumed/plumed_interface.cpp b/source/source_md/module_plumed/plumed_interface.cpp index ca92f8f3d8..1a9ae2fbea 100644 --- a/source/source_md/module_plumed/plumed_interface.cpp +++ b/source/source_md/module_plumed/plumed_interface.cpp @@ -1,9 +1,9 @@ #include "plumed_interface.h" -#include - #include "source_base/tool_quit.h" +#include + #ifdef USE_PLUMED // C API of the PLUMED wrapper library, matching in a // PLUMED installation. The handle is an opaque pointer-sized structure; the @@ -11,7 +11,7 @@ extern "C" { plumed_handle plumed_create(); - void plumed_cmd(plumed_handle plumed, const char* key, const void* val); + void plumed_cmd(plumed_handle plumed, const char *key, const void *val); void plumed_finalize(plumed_handle plumed); } #endif @@ -31,15 +31,13 @@ constexpr double AU_TIME_TO_PS = 2.418884326509e-5; // atomic unit of time constexpr double ELECTRON_MASS_TO_AMU = 5.48579909065e-4; // 1 electron mass = 5.48579909065e-4 amu } // namespace -PlumedInterface::~PlumedInterface() -{ - this->finalize(); -} +PlumedInterface::~PlumedInterface() { this->finalize(); } -bool PlumedInterface::init(const std::string& plumed_file, - const int nat, - const double dt_au, - const double* masses_au) +bool PlumedInterface::init( + const std::string &plumed_file, + const int nat, + const double dt_au, + const double *masses_au) { #ifdef USE_PLUMED this->active_run_ = true; @@ -57,8 +55,8 @@ bool PlumedInterface::init(const std::string& plumed_file, this->has_plumed_ = true; if (this->plumed_.p == nullptr) { - ModuleBase::WARNING_QUIT("ModulePlumed::PlumedInterface::init", - "could not create a PLUMED object, check the PLUMED installation"); + ModuleBase::WARNING_QUIT( + "ModulePlumed::PlumedInterface::init", "could not create a PLUMED object, check the PLUMED installation"); } int real_precision = 8; // double precision reals @@ -67,7 +65,7 @@ bool PlumedInterface::init(const std::string& plumed_file, double energy_units = HARTREE_TO_KJMOL; double time_units = AU_TIME_TO_PS; double timestep = dt_au; - const char* engine = "abacus"; + const char *engine = "abacus"; // commands that configure PLUMED must be issued before "init" plumed_cmd(this->plumed_, "setRealPrecision", &real_precision); @@ -87,14 +85,15 @@ bool PlumedInterface::init(const std::string& plumed_file, (void)dt_au; (void)masses_au; - ModuleBase::WARNING_QUIT("ModulePlumed::PlumedInterface::init", - "plumed was requested in INPUT but ABACUS was not compiled with PLUMED support; " - "reconfigure with -DENABLE_PLUMED=ON"); + ModuleBase::WARNING_QUIT( + "ModulePlumed::PlumedInterface::init", + "plumed was requested in INPUT but ABACUS was not compiled with PLUMED support; " + "reconfigure with -DENABLE_PLUMED=ON"); return false; #endif } -void PlumedInterface::prepare_virial_(const double* cell, const double* virial) +void PlumedInterface::prepare_virial_(const double *cell, const double *virial) { // The Box (cell) data receives forces whenever a periodic collective // variable is biased, so PLUMED always needs a valid force buffer here; @@ -113,21 +112,22 @@ void PlumedInterface::prepare_virial_(const double* cell, const double* virial) { return; } - const double volume = std::fabs(cell[0] * (cell[4] * cell[8] - cell[5] * cell[7]) - - cell[1] * (cell[3] * cell[8] - cell[5] * cell[6]) - + cell[2] * (cell[3] * cell[7] - cell[4] * cell[6])); + const double volume = std::fabs( + cell[0] * (cell[4] * cell[8] - cell[5] * cell[7]) - cell[1] * (cell[3] * cell[8] - cell[5] * cell[6]) + + cell[2] * (cell[3] * cell[7] - cell[4] * cell[6])); for (int i = 0; i < 9; ++i) { this->virial_buf_[i] = -virial[i] * volume; } } -void PlumedInterface::compute_forces_(const int istep, - const ModuleBase::Vector3* pos, - double* force, - const double* cell, - const double potential, - const double* virial) +void PlumedInterface::compute_forces_( + const int istep, + const ModuleBase::Vector3 *pos, + double *force, + const double *cell, + const double potential, + const double *virial) { #ifdef USE_PLUMED int step = istep; @@ -136,8 +136,8 @@ void PlumedInterface::compute_forces_(const int istep, plumed_cmd(this->plumed_, "setStep", &step); // Vector3 is a plain {double x, y, z} structure, so the arrays can // be handed to PLUMED as contiguous [nat][3] buffers - plumed_cmd(this->plumed_, "setPositions", - reinterpret_cast(const_cast*>(pos))); + plumed_cmd( + this->plumed_, "setPositions", reinterpret_cast(const_cast *>(pos))); plumed_cmd(this->plumed_, "setMasses", this->masses_amu_.data()); plumed_cmd(this->plumed_, "setBox", cell); plumed_cmd(this->plumed_, "setEnergy", &energy); @@ -165,19 +165,26 @@ void PlumedInterface::compute_forces_(const int istep, #endif } -void PlumedInterface::move(const int istep, - const ModuleBase::Vector3* pos, - ModuleBase::Vector3* force, - const double* cell, - const double potential, - const double* virial) +void PlumedInterface::move( + const int istep, + const ModuleBase::Vector3 *pos, + ModuleBase::Vector3 *force, + const double *cell, + const double potential, + const double *virial) { #ifdef USE_PLUMED if (!this->active_run_) { return; } - this->compute_forces_(istep, pos, reinterpret_cast(force), cell, potential, virial); + this->compute_forces_( + istep, + pos, + reinterpret_cast(force), + cell, + potential, + virial); #else (void)istep; (void)pos; diff --git a/source/source_md/module_plumed/plumed_interface.h b/source/source_md/module_plumed/plumed_interface.h index 8a1d08e4a3..9ba1ca81af 100644 --- a/source/source_md/module_plumed/plumed_interface.h +++ b/source/source_md/module_plumed/plumed_interface.h @@ -1,17 +1,17 @@ #ifndef PLUMED_INTERFACE_H #define PLUMED_INTERFACE_H +#include "source_base/vector3.h" + #include #include -#include "source_base/vector3.h" - /// opaque handle of the PLUMED wrapper library, see /// in a PLUMED installation; it is a plain pointer-sized structure, so it can /// be stored by value and initialized to a null pointer struct plumed_handle { - void* p; + void *p; }; namespace ModulePlumed @@ -57,8 +57,8 @@ class PlumedInterface PlumedInterface() = default; ~PlumedInterface(); - PlumedInterface(const PlumedInterface&) = delete; - PlumedInterface& operator=(const PlumedInterface&) = delete; + PlumedInterface(const PlumedInterface &) = delete; + PlumedInterface &operator=(const PlumedInterface &) = delete; /** * @brief create a PLUMED object and read the PLUMED input file @@ -72,10 +72,11 @@ class PlumedInterface * masses), nat entries * @return true when the interface is ready to be used in the MD loop */ - bool init(const std::string& plumed_file, - const int nat, - const double dt_au, - const double* masses_au); + bool init( + const std::string &plumed_file, + const int nat, + const double dt_au, + const double *masses_au); /** * @brief one MD step: hand over the configuration, let PLUMED compute the @@ -91,34 +92,33 @@ class PlumedInterface * @param virial [in] 3x3 virial in Hartree (=-stress * volume), or * nullptr when the stress is not available */ - void move(const int istep, - const ModuleBase::Vector3* pos, - ModuleBase::Vector3* force, - const double* cell, - const double potential, - const double* virial); + void move( + const int istep, + const ModuleBase::Vector3 *pos, + ModuleBase::Vector3 *force, + const double *cell, + const double potential, + const double *virial); /// @brief finalize PLUMED (flushes the COLVAR and other output files) void finalize(); /// @brief whether PLUMED is active in this run - bool active() const - { - return this->active_run_; - } + bool active() const { return this->active_run_; } private: /// @brief register the configuration of one step with PLUMED and let it /// compute the collective variables and the biasing forces - void compute_forces_(const int istep, - const ModuleBase::Vector3* pos, - double* force, - const double* cell, - const double potential, - const double* virial); + void compute_forces_( + const int istep, + const ModuleBase::Vector3 *pos, + double *force, + const double *cell, + const double potential, + const double *virial); /// @brief build the 3x3 virial (in energy units) handed over to PLUMED - void prepare_virial_(const double* cell, const double* virial); + void prepare_virial_(const double *cell, const double *virial); plumed_handle plumed_{nullptr}; ///< handle returned by plumed_create() bool has_plumed_ = false; ///< whether a PLUMED object is alive diff --git a/source/source_md/msst.cpp b/source/source_md/msst.cpp index 02f311a739..0f7ac77844 100644 --- a/source/source_md/msst.cpp +++ b/source/source_md/msst.cpp @@ -7,12 +7,18 @@ #endif #include "source_base/timer.h" -MSST::MSST(const MD_para& mdp_in, - const bool cal_stress_in, - const bool init_vel, - const int my_rank_in, - MDCell& mdcell_in) - : MD_base(mdp_in, cal_stress_in, init_vel, my_rank_in, mdcell_in) +MSST::MSST( + const MD_para &mdp_in, + const bool cal_stress_in, + const bool init_vel, + const int my_rank_in, + MDCell &mdcell_in) + : MD_base( + mdp_in, + cal_stress_in, + init_vel, + my_rank_in, + mdcell_in) { msst_qmass = mdp.msst_qmass / pow(ModuleBase::ANGSTROM_AU, 4) / pow(ModuleBase::AU_to_MASS, 2); msst_vel = mdp.msst_vel * ModuleBase::ANGSTROM_AU * ModuleBase::AU_to_FS; @@ -29,17 +35,24 @@ MSST::MSST(const MD_para& mdp_in, lag_pos = 0; vsum = 0; - for (const LocalAtom& atom : mdcell.owned_atoms()) totmass += atom.mass; + for (const LocalAtom &atom : mdcell.owned_atoms()) + { + totmass += atom.mass; + } #ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, &totmass, 1, MPI_DOUBLE, MPI_SUM, mdcell.communicator()); + MPI_Allreduce( + MPI_IN_PLACE, + &totmass, + 1, + MPI_DOUBLE, + MPI_SUM, + mdcell.communicator()); #endif } -MSST::~MSST() -{ -} +MSST::~MSST() {} -void MSST::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp) +void MSST::setup(ModuleESolver::ESolver *p_esolver, const std::string &global_readin_dir, DomainDecomposition &decomp) { ModuleBase::TITLE("MSST", "setup"); ModuleBase::timer::start("MSST", "setup"); @@ -67,10 +80,17 @@ void MSST::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_re std::cout << "initial strain rate = " << fac2 << " msst_tscale = " << mdp.msst_tscale << std::endl; - for (LocalAtom& atom : mdcell.owned_atoms()) atom.vel *= sqrt(1.0 - mdp.msst_tscale); + for (LocalAtom &atom : mdcell.owned_atoms()) + { + atom.vel *= sqrt(1.0 - mdp.msst_tscale); + } } - MD_func::compute_stress(mdcell, cal_stress, virial, stress); + MD_func::compute_stress( + mdcell, + cal_stress, + virial, + stress); t_current = MD_func::current_temp(kinetic, mdcell, frozen_freedom_); } @@ -79,7 +99,7 @@ void MSST::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_re return; } -void MSST::first_half(std::ofstream& ofs) +void MSST::first_half(std::ofstream &ofs) { ModuleBase::TITLE("MSST", "first_half"); ModuleBase::timer::start("MSST", "first_half"); @@ -148,7 +168,11 @@ void MSST::second_half() propagate_vel(); vsum = vel_sum(); - MD_func::compute_stress(mdcell, cal_stress, virial, stress); + MD_func::compute_stress( + mdcell, + cal_stress, + virial, + stress); t_current = MD_func::current_temp(kinetic, mdcell, frozen_freedom_); /// propagate the time derivative of volume 1/2 step @@ -162,15 +186,14 @@ void MSST::second_half() return; } - -void MSST::print_md(std::ofstream& ofs, const bool& cal_stress) +void MSST::print_md(std::ofstream &ofs, const bool &cal_stress) { MD_base::print_md(ofs, cal_stress); return; } -void MSST::write_restart(const std::string& global_out_dir) +void MSST::write_restart(const std::string &global_out_dir) { if (!my_rank) { @@ -195,8 +218,7 @@ void MSST::write_restart(const std::string& global_out_dir) return; } - -void MSST::restart(const std::string& global_readin_dir) +void MSST::restart(const std::string &global_readin_dir) { bool ok = true; @@ -214,7 +236,7 @@ void MSST::restart(const std::string& global_readin_dir) if (ok) { file >> step_rst_ >> md_tfirst >> omega[mdp.msst_direction] >> e0 >> v0 >> p0 >> lag_pos; - if(!file) + if (!file) { ok = false; } @@ -223,7 +245,12 @@ void MSST::restart(const std::string& global_readin_dir) } #ifdef __MPI - MPI_Bcast(&ok, 1, MPI_C_BOOL, 0, mdcell.communicator()); + MPI_Bcast( + &ok, + 1, + MPI_C_BOOL, + 0, + mdcell.communicator()); #endif if (!ok) @@ -232,13 +259,48 @@ void MSST::restart(const std::string& global_readin_dir) } #ifdef __MPI - MPI_Bcast(&step_rst_, 1, MPI_INT, 0, mdcell.communicator()); - MPI_Bcast(&md_tfirst, 1, MPI_DOUBLE, 0, mdcell.communicator()); - MPI_Bcast(&omega[mdp.msst_direction], 1, MPI_DOUBLE, 0, mdcell.communicator()); - MPI_Bcast(&e0, 1, MPI_DOUBLE, 0, mdcell.communicator()); - MPI_Bcast(&v0, 1, MPI_DOUBLE, 0, mdcell.communicator()); - MPI_Bcast(&p0, 1, MPI_DOUBLE, 0, mdcell.communicator()); - MPI_Bcast(&lag_pos, 1, MPI_DOUBLE, 0, mdcell.communicator()); + MPI_Bcast( + &step_rst_, + 1, + MPI_INT, + 0, + mdcell.communicator()); + MPI_Bcast( + &md_tfirst, + 1, + MPI_DOUBLE, + 0, + mdcell.communicator()); + MPI_Bcast( + &omega[mdp.msst_direction], + 1, + MPI_DOUBLE, + 0, + mdcell.communicator()); + MPI_Bcast( + &e0, + 1, + MPI_DOUBLE, + 0, + mdcell.communicator()); + MPI_Bcast( + &v0, + 1, + MPI_DOUBLE, + 0, + mdcell.communicator()); + MPI_Bcast( + &p0, + 1, + MPI_DOUBLE, + 0, + mdcell.communicator()); + MPI_Bcast( + &lag_pos, + 1, + MPI_DOUBLE, + 0, + mdcell.communicator()); #endif return; @@ -248,14 +310,23 @@ double MSST::vel_sum() const { double vsum = 0; - for (const LocalAtom& atom : mdcell.owned_atoms()) vsum += atom.vel.norm2(); + for (const LocalAtom &atom : mdcell.owned_atoms()) + { + vsum += atom.vel.norm2(); + } #ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, &vsum, 1, MPI_DOUBLE, MPI_SUM, mdcell.communicator()); + MPI_Allreduce( + MPI_IN_PLACE, + &vsum, + 1, + MPI_DOUBLE, + MPI_SUM, + mdcell.communicator()); #endif return vsum; } -void MSST::rescale(std::ofstream& ofs, const double& volume) +void MSST::rescale(std::ofstream &ofs, const double &volume) { int sd = mdp.msst_direction; @@ -270,18 +341,20 @@ void MSST::rescale(std::ofstream& ofs, const double& volume) mdcell.refresh_cart_from_frac(); /// rescale velocity - for (LocalAtom& atom : mdcell.owned_atoms()) atom.vel[sd] *= dilation[sd]; + for (LocalAtom &atom : mdcell.owned_atoms()) + { + atom.vel[sd] *= dilation[sd]; + } static_cast(ofs); } - void MSST::propagate_vel() { const int sd = mdp.msst_direction; const double dthalf = 0.5 * md_dt; const double fac = msst_vis * pow(omega[sd], 2) / (vsum * mdcell.omega()); - for (LocalAtom& atom : mdcell.owned_atoms()) + for (LocalAtom &atom : mdcell.owned_atoms()) { ModuleBase::Vector3 const_C = atom.force / atom.mass; ModuleBase::Vector3 const_D; @@ -297,16 +370,15 @@ void MSST::propagate_vel() } else { - atom.vel[k] - += (const_C[k] + const_D[k] * atom.vel[k]) * dthalf - + 0.5 * (const_D[k] * const_D[k] * atom.vel[k] + const_C[k] * const_D[k]) * dthalf * dthalf; + atom.vel[k] += + (const_C[k] + const_D[k] * atom.vel[k]) * dthalf + + 0.5 * (const_D[k] * const_D[k] * atom.vel[k] + const_C[k] * const_D[k]) * dthalf * dthalf; } } } return; } - void MSST::propagate_voldot() { const int sd = mdp.msst_direction; @@ -330,8 +402,8 @@ void MSST::propagate_voldot() } else { - omega[sd] += (const_A - const_B * omega[sd]) * dthalf - + 0.5 * (const_B * const_B * omega[sd] - const_A * const_B) * dthalf * dthalf; + omega[sd] += (const_A - const_B * omega[sd]) * dthalf + + 0.5 * (const_B * const_B * omega[sd] - const_A * const_B) * dthalf * dthalf; } return; diff --git a/source/source_md/msst.h b/source/source_md/msst.h index 00ba2a7914..68f63b613e 100644 --- a/source/source_md/msst.h +++ b/source/source_md/msst.h @@ -2,6 +2,7 @@ #define MSST_H #include "md_base.h" + #include /** @@ -15,50 +16,36 @@ class MSST : public MD_base { public: - MSST(const MD_para& mdp_in, - const bool cal_stress_in, - const bool init_vel, - const int my_rank_in, - MDCell& mdcell_in); + MSST( + const MD_para &mdp_in, + const bool cal_stress_in, + const bool init_vel, + const int my_rank_in, + MDCell &mdcell_in); ~MSST(); /// @brief time derivative of the volume, per lattice direction - const ModuleBase::Vector3& get_omega() const - { - return omega; - } + const ModuleBase::Vector3 &get_omega() const { return omega; } /// @brief energy of the initial configuration - double get_e0() const - { - return e0; - } + double get_e0() const { return e0; } /// @brief volume of the initial configuration - double get_v0() const - { - return v0; - } + double get_v0() const { return v0; } /// @brief pressure of the initial configuration - double get_p0() const - { - return p0; - } + double get_p0() const { return p0; } /// @brief Lagrangian location of the cell - double get_lag_pos() const - { - return lag_pos; - } + double get_lag_pos() const { return lag_pos; } private: - void setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp); - void first_half(std::ofstream& ofs); + void setup(ModuleESolver::ESolver *p_esolver, const std::string &global_readin_dir, DomainDecomposition &decomp); + void first_half(std::ofstream &ofs); void second_half(); - void print_md(std::ofstream& ofs, const bool& cal_stress); - void write_restart(const std::string& global_out_dir); - void restart(const std::string& global_readin_dir); + void print_md(std::ofstream &ofs, const bool &cal_stress); + void write_restart(const std::string &global_out_dir); + void restart(const std::string &global_readin_dir); /** * @brief get the sum of square of velocities @@ -73,7 +60,7 @@ class MSST : public MD_base * @param ofs determine the output files * @param volume the current cell volume */ - void rescale(std::ofstream& ofs, const double& volume); + void rescale(std::ofstream &ofs, const double &volume); /** * @brief propagate atomic velocities @@ -87,18 +74,18 @@ class MSST : public MD_base */ void propagate_voldot(void); - std::vector > old_v; ///< old atomic velocities - ModuleBase::Vector3 dilation; ///< dilation scale - ModuleBase::Vector3 omega; ///< time derivative of volume - double p0; ///< initial pressure - double v0; ///< initial volume - double e0; ///< initial energy - double totmass; ///< total mass of the cell - double lag_pos; ///< Lagrangian location of cell - double vsum; ///< sum over v^2 - double msst_vel; ///< shock msst_vel (\AA/fs) - double msst_qmass; ///< cell mass-like parameter (mass^2/length^4) - double msst_vis; ///< artificial msst_vis (mass/length/time) + std::vector> old_v; ///< old atomic velocities + ModuleBase::Vector3 dilation; ///< dilation scale + ModuleBase::Vector3 omega; ///< time derivative of volume + double p0; ///< initial pressure + double v0; ///< initial volume + double e0; ///< initial energy + double totmass; ///< total mass of the cell + double lag_pos; ///< Lagrangian location of cell + double vsum; ///< sum over v^2 + double msst_vel; ///< shock msst_vel (\AA/fs) + double msst_qmass; ///< cell mass-like parameter (mass^2/length^4) + double msst_vis; ///< artificial msst_vis (mass/length/time) }; #endif diff --git a/source/source_md/nhchain.cpp b/source/source_md/nhchain.cpp index 9a99f1b05e..97151b80d5 100644 --- a/source/source_md/nhchain.cpp +++ b/source/source_md/nhchain.cpp @@ -6,18 +6,24 @@ #endif #include "source_base/timer.h" -Nose_Hoover::Nose_Hoover(const MD_para& mdp_in, - const bool cal_stress_in, - const bool init_vel, - const int my_rank_in, - MDCell& mdcell_in) - : MD_base(mdp_in, cal_stress_in, init_vel, my_rank_in, mdcell_in) +Nose_Hoover::Nose_Hoover( + const MD_para &mdp_in, + const bool cal_stress_in, + const bool init_vel, + const int my_rank_in, + MDCell &mdcell_in) + : MD_base( + mdp_in, + cal_stress_in, + init_vel, + my_rank_in, + mdcell_in) { const double unit_transform = ModuleBase::HARTREE_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; md_tfreq = mdp.md_tfreq * ModuleBase::AU_to_FS; - assert(unit_transform>0.0); + assert(unit_transform > 0.0); md_pfirst = mdp.md_pfirst / unit_transform; md_plast = mdp.md_plast / unit_transform; @@ -64,7 +70,7 @@ Nose_Hoover::Nose_Hoover(const MD_para& mdp_in, */ else if (mdp.md_pmode == "tri") { - const ModuleBase::Matrix3& latvec = mdcell.latvec(); + const ModuleBase::Matrix3 &latvec = mdcell.latvec(); if (latvec.e12 || latvec.e13 || latvec.e23) { ModuleBase::WARNING_QUIT("Nose_Hoover", "the lattice must be lower-triangular when md_pmode == tri!"); @@ -95,7 +101,7 @@ Nose_Hoover::Nose_Hoover(const MD_para& mdp_in, tdof = MD_func::global_dof(mdcell); - assert(mdp.md_tchain>0); + assert(mdp.md_tchain > 0); /// allocate thermostats coupled with particles mass_eta = new double[mdp.md_tchain]; @@ -159,7 +165,8 @@ Nose_Hoover::~Nose_Hoover() } } -void Nose_Hoover::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp) +void Nose_Hoover::setup( + ModuleESolver::ESolver *p_esolver, const std::string &global_readin_dir, DomainDecomposition &decomp) { ModuleBase::TITLE("Nose_Hoover", "setup"); ModuleBase::timer::start("Nose_Hoover", "setup"); @@ -167,7 +174,11 @@ void Nose_Hoover::setup(ModuleESolver::ESolver* p_esolver, const std::string& gl MD_base::setup(p_esolver, global_readin_dir, decomp); /// determine target temperature - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); + t_target = MD_func::target_temp( + step_ + step_rst_, + mdp.md_nstep, + md_tfirst, + md_tlast); /// init thermostats coupled with particles mass_eta[0] = tdof * t_target / md_tfreq / md_tfreq; @@ -214,7 +225,7 @@ void Nose_Hoover::setup(ModuleESolver::ESolver* p_esolver, const std::string& gl return; } -void Nose_Hoover::first_half(std::ofstream& ofs) +void Nose_Hoover::first_half(std::ofstream &ofs) { ModuleBase::TITLE("Nose_Hoover", "first_half"); ModuleBase::timer::start("Nose_Hoover", "first_half"); @@ -226,7 +237,11 @@ void Nose_Hoover::first_half(std::ofstream& ofs) } /// update target T - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); + t_target = MD_func::target_temp( + step_ + step_rst_, + mdp.md_nstep, + md_tfirst, + md_tlast); /// update thermostats coupled with particles particle_thermo(); @@ -235,7 +250,11 @@ void Nose_Hoover::first_half(std::ofstream& ofs) { /// update temperature and stress due to velocity rescaling t_current = MD_func::current_temp(kinetic, mdcell, frozen_freedom_); - MD_func::compute_stress(mdcell, cal_stress, virial, stress); + MD_func::compute_stress( + mdcell, + cal_stress, + virial, + stress); /// couple stress component due to md_pcouple couple_stress(); @@ -273,7 +292,6 @@ void Nose_Hoover::first_half(std::ofstream& ofs) return; } - void Nose_Hoover::second_half() { ModuleBase::TITLE("Nose_Hoover", "second_half"); @@ -294,7 +312,11 @@ void Nose_Hoover::second_half() if (npt_flag) { /// update stress due to velocity rescaling - MD_func::compute_stress(mdcell, cal_stress, virial, stress); + MD_func::compute_stress( + mdcell, + cal_stress, + virial, + stress); /// couple stress component due to md_pcouple couple_stress(); @@ -317,13 +339,9 @@ void Nose_Hoover::second_half() return; } +void Nose_Hoover::print_md(std::ofstream &ofs, const bool &cal_stress) { MD_base::print_md(ofs, cal_stress); } -void Nose_Hoover::print_md(std::ofstream& ofs, const bool& cal_stress) -{ - MD_base::print_md(ofs, cal_stress); -} - -void Nose_Hoover::write_restart(const std::string& global_out_dir) +void Nose_Hoover::write_restart(const std::string &global_out_dir) { if (!my_rank) { @@ -373,7 +391,7 @@ void Nose_Hoover::write_restart(const std::string& global_out_dir) #endif } -void Nose_Hoover::restart(const std::string& global_readin_dir) +void Nose_Hoover::restart(const std::string &global_readin_dir) { bool ok = true; bool ok2 = true; @@ -444,9 +462,24 @@ void Nose_Hoover::restart(const std::string& global_readin_dir) } #ifdef __MPI - MPI_Bcast(&ok, 1, MPI_C_BOOL, 0, mdcell.communicator()); - MPI_Bcast(&ok2, 1, MPI_C_BOOL, 0, mdcell.communicator()); - MPI_Bcast(&ok3, 1, MPI_C_BOOL, 0, mdcell.communicator()); + MPI_Bcast( + &ok, + 1, + MPI_C_BOOL, + 0, + mdcell.communicator()); + MPI_Bcast( + &ok2, + 1, + MPI_C_BOOL, + 0, + mdcell.communicator()); + MPI_Bcast( + &ok3, + 1, + MPI_C_BOOL, + 0, + mdcell.communicator()); #endif if (!ok) @@ -463,15 +496,50 @@ void Nose_Hoover::restart(const std::string& global_readin_dir) } #ifdef __MPI - MPI_Bcast(&step_rst_, 1, MPI_INT, 0, mdcell.communicator()); - MPI_Bcast(&md_tfirst, 1, MPI_DOUBLE, 0, mdcell.communicator()); - MPI_Bcast(eta, mdp.md_tchain, MPI_DOUBLE, 0, mdcell.communicator()); - MPI_Bcast(v_eta, mdp.md_tchain, MPI_DOUBLE, 0, mdcell.communicator()); + MPI_Bcast( + &step_rst_, + 1, + MPI_INT, + 0, + mdcell.communicator()); + MPI_Bcast( + &md_tfirst, + 1, + MPI_DOUBLE, + 0, + mdcell.communicator()); + MPI_Bcast( + eta, + mdp.md_tchain, + MPI_DOUBLE, + 0, + mdcell.communicator()); + MPI_Bcast( + v_eta, + mdp.md_tchain, + MPI_DOUBLE, + 0, + mdcell.communicator()); if (npt_flag) { - MPI_Bcast(v_omega, 6, MPI_DOUBLE, 0, mdcell.communicator()); - MPI_Bcast(peta, mdp.md_pchain, MPI_DOUBLE, 0, mdcell.communicator()); - MPI_Bcast(v_peta, mdp.md_pchain, MPI_DOUBLE, 0, mdcell.communicator()); + MPI_Bcast( + v_omega, + 6, + MPI_DOUBLE, + 0, + mdcell.communicator()); + MPI_Bcast( + peta, + mdp.md_pchain, + MPI_DOUBLE, + 0, + mdcell.communicator()); + MPI_Bcast( + v_peta, + mdp.md_pchain, + MPI_DOUBLE, + 0, + mdcell.communicator()); } #endif } @@ -555,7 +623,10 @@ void Nose_Hoover::particle_thermo() } /// rescale velocity due to thermostats - for (LocalAtom& atom : mdcell.owned_atoms()) atom.vel *= scale; + for (LocalAtom &atom : mdcell.owned_atoms()) + { + atom.vel *= scale; + } } void Nose_Hoover::baro_thermo() @@ -649,12 +720,24 @@ void Nose_Hoover::update_baro() else { ModuleBase::matrix t_vector(3, 3); - for (const LocalAtom& atom : mdcell.owned_atoms()) + for (const LocalAtom &atom : mdcell.owned_atoms()) + { for (int i = 0; i < 3; ++i) + { for (int j = 0; j < 3; ++j) + { t_vector(i, j) += atom.mass * atom.vel[i] * atom.vel[j]; + } + } + } #ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, t_vector.c, t_vector.nr * t_vector.nc, MPI_DOUBLE, MPI_SUM, mdcell.communicator()); + MPI_Allreduce( + MPI_IN_PLACE, + t_vector.c, + t_vector.nr * t_vector.nc, + MPI_DOUBLE, + MPI_SUM, + mdcell.communicator()); #endif for (int i = 0; i < 3; ++i) @@ -700,7 +783,7 @@ void Nose_Hoover::vel_baro() factor[i] = exp(-(v_omega[i] + mtk_term) * md_dt / 4); } - for (LocalAtom& atom : mdcell.owned_atoms()) + for (LocalAtom &atom : mdcell.owned_atoms()) { for (int j = 0; j < 3; ++j) { @@ -724,7 +807,7 @@ void Nose_Hoover::vel_baro() } } -void Nose_Hoover::update_volume(std::ofstream& ofs) +void Nose_Hoover::update_volume(std::ofstream &ofs) { double factor = 0.0; ModuleBase::Matrix3 latvec = mdcell.latvec(); diff --git a/source/source_md/nhchain.h b/source/source_md/nhchain.h index 6e0db3bec6..335e3bcc19 100644 --- a/source/source_md/nhchain.h +++ b/source/source_md/nhchain.h @@ -12,54 +12,40 @@ class Nose_Hoover : public MD_base { public: - Nose_Hoover(const MD_para& mdp_in, - const bool cal_stress_in, - const bool init_vel, - const int my_rank_in, - MDCell& mdcell_in); + Nose_Hoover( + const MD_para &mdp_in, + const bool cal_stress_in, + const bool init_vel, + const int my_rank_in, + MDCell &mdcell_in); ~Nose_Hoover(); /// @brief positions of the thermostats coupled with the particles /// (md_tchain entries) - const double* get_eta() const - { - return eta; - } + const double *get_eta() const { return eta; } /// @brief velocities of the thermostats coupled with the particles /// (md_tchain entries) - const double* get_v_eta() const - { - return v_eta; - } + const double *get_v_eta() const { return v_eta; } /// @brief positions of the thermostats coupled with the barostat /// (md_pchain entries) - const double* get_peta() const - { - return peta; - } + const double *get_peta() const { return peta; } /// @brief velocities of the thermostats coupled with the barostat /// (md_pchain entries) - const double* get_v_peta() const - { - return v_peta; - } + const double *get_v_peta() const { return v_peta; } /// @brief velocities of the six lattice components - const double* get_v_omega() const - { - return v_omega; - } + const double *get_v_omega() const { return v_omega; } private: - void setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp); - void first_half(std::ofstream& ofs); + void setup(ModuleESolver::ESolver *p_esolver, const std::string &global_readin_dir, DomainDecomposition &decomp); + void first_half(std::ofstream &ofs); void second_half(); - void print_md(std::ofstream& ofs, const bool& cal_stress); - void write_restart(const std::string& global_out_dir); - void restart(const std::string& global_readin_dir); + void print_md(std::ofstream &ofs, const bool &cal_stress); + void write_restart(const std::string &global_out_dir); + void restart(const std::string &global_readin_dir); /** * @brief perform half-step update of thermostats coupled with particles @@ -102,40 +88,40 @@ class Nose_Hoover : public MD_base * * @param ofs determine the output files */ - void update_volume(std::ofstream& ofs); + void update_volume(std::ofstream &ofs); const int nc_tchain = 1; ///< the number of nc for nvt const int nc_pchain = 1; ///< the number of nc for npt const static int nys = 7; ///< the number of scale evolution operator double w[nys]; ///< scale evolution operator - std::int64_t tdof; ///< particle degree of freedom - double t_target=0.0;///< target temperature - double* mass_eta = nullptr; ///< mass of thermostats coupled with particles - double* eta = nullptr; ///< position of thermostats coupled with particles - double* v_eta = nullptr; ///< velocity of thermostats coupled with particles - double* g_eta = nullptr; ///< acceleration of thermostats coupled with particles - - int npt_flag; ///< whether NPT ensemble - double mass_omega[6]; ///< mass of lattice component - double v_omega[6]; ///< velocity of lattice component - double pstart[6]; ///< initial stress components - double pstop[6]; ///< final stress components - double pfreq[6]; ///< Oscillation frequency, used to determine qmass of thermostats coupled with barostat - int pflag[6]; ///< control stress components - int pdim; ///< pdim = pflag[0] + pflag[1] + pflag[2], number of barostatted dims - double p_target[6]; ///< target stress components - double p_hydro = 0.0; ///< target hydrostatic target pressure + std::int64_t tdof; ///< particle degree of freedom + double t_target = 0.0; ///< target temperature + double *mass_eta = nullptr; ///< mass of thermostats coupled with particles + double *eta = nullptr; ///< position of thermostats coupled with particles + double *v_eta = nullptr; ///< velocity of thermostats coupled with particles + double *g_eta = nullptr; ///< acceleration of thermostats coupled with particles + + int npt_flag; ///< whether NPT ensemble + double mass_omega[6]; ///< mass of lattice component + double v_omega[6]; ///< velocity of lattice component + double pstart[6]; ///< initial stress components + double pstop[6]; ///< final stress components + double pfreq[6]; ///< Oscillation frequency, used to determine qmass of thermostats coupled with barostat + int pflag[6]; ///< control stress components + int pdim; ///< pdim = pflag[0] + pflag[1] + pflag[2], number of barostatted dims + double p_target[6]; ///< target stress components + double p_hydro = 0.0; ///< target hydrostatic target pressure double p_current[6] = {0.0}; ///< current stress after coupled - double* mass_peta = nullptr; ///< mass of thermostats coupled with barostat - double* peta = nullptr; ///< position of thermostats coupled with barostat - double* v_peta = nullptr; ///< velocity of thermostats coupled with barostat - double* g_peta = nullptr; ///< acceleration of thermostats coupled with barostat - double mtk_term=0; ///< mtk correction - double md_tfreq; ///< Oscillation frequency, used to determine qmass of thermostats coupled with particles - double md_pfirst; ///< Initial pressure - double md_plast; ///< Final pressure - double md_pfreq; ///< Oscillation frequency, used to determine qmass of thermostats coupled with barostat + double *mass_peta = nullptr; ///< mass of thermostats coupled with barostat + double *peta = nullptr; ///< position of thermostats coupled with barostat + double *v_peta = nullptr; ///< velocity of thermostats coupled with barostat + double *g_peta = nullptr; ///< acceleration of thermostats coupled with barostat + double mtk_term = 0; ///< mtk correction + double md_tfreq; ///< Oscillation frequency, used to determine qmass of thermostats coupled with particles + double md_pfirst; ///< Initial pressure + double md_plast; ///< Final pressure + double md_pfreq; ///< Oscillation frequency, used to determine qmass of thermostats coupled with barostat }; #endif diff --git a/source/source_md/run_md.cpp b/source/source_md/run_md.cpp index 114f9cbb3a..426ddf9ebf 100644 --- a/source/source_md/run_md.cpp +++ b/source/source_md/run_md.cpp @@ -1,28 +1,28 @@ #include "run_md.h" +#include "fire.h" +#include "langevin.h" +#include "md_func.h" #include "source_base/constants.h" +#include "source_base/global_file.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" #include "source_base/parallel_cell.h" -#include "source_cell/mdcell_reader.h" -#include "source_cell/mdcell.h" +#include "source_base/timer.h" #include "source_cell/domain_decomposition.h" +#include "source_cell/mdcell.h" +#include "source_cell/mdcell_reader.h" #include "source_io/module_parameter/parameter.h" -#include "fire.h" -#include "langevin.h" -#include "md_func.h" -#include "source_base/global_file.h" -#include "source_base/timer.h" #ifdef __JSON #include "source_io/module_json/output_info.h" #endif -#include "source_io/module_output/print_info.h" #include "msst.h" #include "nhchain.h" -#include "verlet.h" -#include "source_cell/update_cell.h" #include "source_cell/print_cell.h" +#include "source_cell/update_cell.h" +#include "source_io/module_output/print_info.h" #include "source_md/module_plumed/plumed_interface.h" +#include "verlet.h" #include #ifdef __MPI @@ -32,9 +32,9 @@ namespace Run_MD { -void prepare_mdcell(MDCell& mdcell, const Parameter& param_in, DomainDecomposition& decomp) +void prepare_mdcell(MDCell &mdcell, const Parameter ¶m_in, DomainDecomposition &decomp) { - const Input_para& input = param_in.inp; + const Input_para &input = param_in.inp; std::vector effective_replicate = input.cell_replica; if (input.mdp.md_restart) { @@ -42,20 +42,26 @@ void prepare_mdcell(MDCell& mdcell, const Parameter& param_in, DomainDecompositi } const ModuleBase::CommunicationDomain comm_domain = ModuleBase::world_comm_domain(); - mdcell = MDCellReader::read_stru(param_in.globalv.global_in_stru, - effective_replicate, - input.mdp.md_neighbor_skin / ModuleBase::BOHR_TO_A, - comm_domain, - decomp); + mdcell = MDCellReader::read_stru( + param_in.globalv.global_in_stru, + effective_replicate, + input.mdp.md_neighbor_skin / ModuleBase::BOHR_TO_A, + comm_domain, + decomp); GlobalV::ofs_running << std::endl; ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "TOTAL ATOM NUMBER", mdcell.nat()); GlobalV::ofs_running << std::endl; } -void prepare_mdcell(MDCell& mdcell, UnitCell& ucell, DomainDecomposition& decomp) +void prepare_mdcell(MDCell &mdcell, UnitCell &ucell, DomainDecomposition &decomp) { const ModuleBase::CommunicationDomain comm_domain = ModuleBase::world_comm_domain(); - decomp.init(comm_domain, ucell.latvec, ucell.lat0, 0.0, 0.0); + decomp.init( + comm_domain, + ucell.latvec, + ucell.lat0, + 0.0, + 0.0); const std::vector owned_atoms = decomp.split_owned_atoms_from_ucell(ucell); std::vector type_labels; std::vector type_masses; @@ -66,17 +72,18 @@ void prepare_mdcell(MDCell& mdcell, UnitCell& ucell, DomainDecomposition& decomp type_masses.push_back(ucell.atoms[it].mass); type_atom_counts.push_back(ucell.atoms[it].na); } - mdcell.initialize_from_owned_atoms(ucell.latvec, - ucell.GT, - ucell.lat0, - ucell.omega, - ucell.nat, - owned_atoms, - type_labels, - type_masses, - type_atom_counts, - 0.0, - comm_domain); + mdcell.initialize_from_owned_atoms( + ucell.latvec, + ucell.GT, + ucell.lat0, + ucell.omega, + ucell.nat, + owned_atoms, + type_labels, + type_masses, + type_atom_counts, + 0.0, + comm_domain); mdcell.set_backing_unitcell(ucell); mdcell.mutable_stru_meta() = unitcell::make_stru_meta(ucell); } @@ -85,10 +92,10 @@ namespace { /// @brief collect the absolute Cartesian positions (in Bohr) of the atoms /// owned by this rank; the cell stores cart in units of the lattice constant -void collect_plumed_positions(const MDCell& mdcell, std::vector>& positions) +void collect_plumed_positions(const MDCell &mdcell, std::vector> &positions) { const double lat0 = mdcell.lat0(); - const std::vector& atoms = mdcell.owned_atoms(); + const std::vector &atoms = mdcell.owned_atoms(); positions.resize(atoms.size()); for (std::size_t i = 0; i < atoms.size(); ++i) { @@ -97,9 +104,9 @@ void collect_plumed_positions(const MDCell& mdcell, std::vector>& forces) +void collect_plumed_forces(const MDCell &mdcell, std::vector> &forces) { - const std::vector& atoms = mdcell.owned_atoms(); + const std::vector &atoms = mdcell.owned_atoms(); forces.resize(atoms.size()); for (std::size_t i = 0; i < atoms.size(); ++i) { @@ -108,9 +115,9 @@ void collect_plumed_forces(const MDCell& mdcell, std::vector>& forces) +void apply_plumed_forces(MDCell &mdcell, const std::vector> &forces) { - std::vector& atoms = mdcell.owned_atoms(); + std::vector &atoms = mdcell.owned_atoms(); for (std::size_t i = 0; i < atoms.size(); ++i) { atoms[i].force = forces[i]; @@ -118,10 +125,10 @@ void apply_plumed_forces(MDCell& mdcell, const std::vector> plumed_forces; double plumed_cell[9]; - auto run_plumed_step = [&](const int istep) { + auto run_plumed_step = [&](const int istep) + { collect_plumed_positions(mdcell, plumed_positions); collect_plumed_forces(mdcell, plumed_forces); collect_plumed_cell(mdcell, plumed_cell); - plumed.move(istep, - plumed_positions.data(), - plumed_forces.data(), - plumed_cell, - mdrun->potential, - param_in.inp.cal_stress ? mdrun->virial.c : nullptr); + plumed.move( + istep, + plumed_positions.data(), + plumed_forces.data(), + plumed_cell, + mdrun->potential, + param_in.inp.cal_stress ? mdrun->virial.c : nullptr); apply_plumed_forces(mdcell, plumed_forces); }; @@ -212,19 +248,20 @@ void md_line(MDCell& mdcell, #endif if (nproc > 1) { - ModuleBase::WARNING_QUIT("md_line", - "the PLUMED interface currently supports a single MPI rank only; " - "run the MD with one process (OpenMP threads are fine)"); + ModuleBase::WARNING_QUIT( + "md_line", "the PLUMED interface currently supports a single MPI rank only; " + "run the MD with one process (OpenMP threads are fine)"); } std::vector plumed_masses(mdcell.owned_atoms().size()); for (std::size_t i = 0; i < mdcell.owned_atoms().size(); ++i) { plumed_masses[i] = mdcell.owned_atoms()[i].mass; } - plumed.init(param_in.mdp.plumed_file, - static_cast(mdcell.owned_atoms().size()), - mdrun->get_md_dt(), - plumed_masses.data()); + plumed.init( + param_in.mdp.plumed_file, + static_cast(mdcell.owned_atoms().size()), + mdrun->get_md_dt(), + plumed_masses.data()); /// the initial configuration already carries positions and /// forces, so PLUMED can be called right away run_plumed_step(mdrun->step_ + mdrun->step_rst_); @@ -240,14 +277,15 @@ void md_line(MDCell& mdcell, mdrun->first_half(GlobalV::ofs_running); /// update force and virial due to the update of atom positions - MD_func::force_virial(p_esolver, - mdrun->step_, - mdcell, - decomp, - mdrun->potential, - param_in.inp.cal_stress, - mdrun->virial, - param_in.mdp.md_out_force); + MD_func::force_virial( + p_esolver, + mdrun->step_, + mdcell, + decomp, + mdrun->potential, + param_in.inp.cal_stress, + mdrun->virial, + param_in.mdp.md_out_force); if (plumed.active()) { @@ -259,29 +297,27 @@ void md_line(MDCell& mdcell, mdrun->second_half(); - MD_func::compute_stress(mdcell, - param_in.inp.cal_stress, - mdrun->virial, - mdrun->stress); - mdrun->t_current = MD_func::current_temp(mdrun->kinetic, - mdcell, - mdrun->frozen_freedom_); + MD_func::compute_stress( + mdcell, + param_in.inp.cal_stress, + mdrun->virial, + mdrun->stress); + mdrun->t_current = MD_func::current_temp(mdrun->kinetic, mdcell, mdrun->frozen_freedom_); } mdrun->print_md(GlobalV::ofs_running, PARAM.inp.cal_stress); - if (param_in.mdp.md_dumpfreq > 0 - && (mdrun->step_ + mdrun->step_rst_) % param_in.mdp.md_dumpfreq == 0) + if (param_in.mdp.md_dumpfreq > 0 && (mdrun->step_ + mdrun->step_rst_) % param_in.mdp.md_dumpfreq == 0) { - MD_func::dump_info(mdrun->step_ + mdrun->step_rst_, - PARAM.globalv.global_out_dir, - mdcell, - mdp, - cal_stress, - mdrun->virial); + MD_func::dump_info( + mdrun->step_ + mdrun->step_rst_, + PARAM.globalv.global_out_dir, + mdcell, + mdp, + cal_stress, + mdrun->virial); } - if (param_in.mdp.md_restartfreq > 0 - && (mdrun->step_ + mdrun->step_rst_) % param_in.mdp.md_restartfreq == 0) + if (param_in.mdp.md_restartfreq > 0 && (mdrun->step_ + mdrun->step_rst_) % param_in.mdp.md_restartfreq == 0) { if (mdcell.has_backing_unitcell()) { diff --git a/source/source_md/run_md.h b/source/source_md/run_md.h index 8abf531a57..eb46a02e58 100644 --- a/source/source_md/run_md.h +++ b/source/source_md/run_md.h @@ -17,9 +17,9 @@ class ESolver; */ namespace Run_MD { -void prepare_mdcell(MDCell& mdcell, const Parameter& param_in, DomainDecomposition& decomp); +void prepare_mdcell(MDCell &mdcell, const Parameter ¶m_in, DomainDecomposition &decomp); -void prepare_mdcell(MDCell& mdcell, UnitCell& ucell, DomainDecomposition& decomp); +void prepare_mdcell(MDCell &mdcell, UnitCell &ucell, DomainDecomposition &decomp); /** * @brief the md loop line @@ -28,10 +28,11 @@ void prepare_mdcell(MDCell& mdcell, UnitCell& ucell, DomainDecomposition& decomp * @param p_esolver energy solver * @param md_para input parameters used in md */ -void md_line(MDCell& mdcell, - ModuleESolver::ESolver* p_esolver, - const Parameter& param_in, - DomainDecomposition& decomp); +void md_line( + MDCell &mdcell, + ModuleESolver::ESolver *p_esolver, + const Parameter ¶m_in, + DomainDecomposition &decomp); } // namespace Run_MD #endif diff --git a/source/source_md/test/fire_test.cpp b/source/source_md/test/fire_test.cpp index 8655b5aa90..85c253c3c8 100644 --- a/source/source_md/test/fire_test.cpp +++ b/source/source_md/test/fire_test.cpp @@ -1,9 +1,11 @@ +#include "source_md/fire.h" + +#include "setcell.h" #include "source_cell/domain_decomposition.h" +#include "source_esolver/esolver_lj.h" + #include "gmock/gmock.h" #include "gtest/gtest.h" -#include "setcell.h" -#include "source_esolver/esolver_lj.h" -#include "source_md/fire.h" #define doublethreshold 1e-12 /************************************************ @@ -34,12 +36,12 @@ class FIREtest : public testing::Test { protected: - MD_base* mdrun; + MD_base *mdrun; UnitCell ucell; MDCell mdcell; DomainDecomposition decomp; Input_para inp; - ModuleESolver::ESolver* p_esolver; + ModuleESolver::ESolver *p_esolver; /// the md writers take the directory explicitly, so the test owns it const std::string out_dir = "./"; const std::string readin_dir = "./"; @@ -51,9 +53,20 @@ class FIREtest : public testing::Test p_esolver = new ModuleESolver::ESolver_LJ(); mdcell = Setcell::setup_mdcell(ucell); - decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0); + decomp.init( + ModuleBase::world_comm_domain(), + mdcell.latvec(), + mdcell.lat0(), + 0.0, + 0.0); p_esolver->before_all_runners(mdcell, inp); - mdrun = new FIRE(inp.mdp, inp.cal_stress, inp.init_vel, /*my_rank=*/0, inp.force_thr, mdcell); + mdrun = new FIRE( + inp.mdp, + inp.cal_stress, + inp.init_vel, + /*my_rank=*/0, + inp.force_thr, + mdcell); mdrun->setup(p_esolver, readin_dir, decomp); } @@ -82,18 +95,42 @@ TEST_F(FIREtest, FirstHalf) { mdrun->first_half(GlobalV::ofs_running); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00045447059554315662, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00032646833232493271, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.215709523063016e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.0005213674681407162, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00059486888444406608, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00035886062145122004, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00052406920303529794, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, 4.8706739346586155e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00020129054406946794, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00045717233044145918, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00021969381277291936, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00010541298215149392, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00045447059554315662, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00032646833232493271, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.215709523063016e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.0005213674681407162, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00059486888444406608, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00035886062145122004, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00052406920303529794, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, 4.8706739346586155e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00020129054406946794, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00045717233044145918, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00021969381277291936, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00010541298215149392, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -0.00010993118004167345, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 7.8968913216100539e-05, doublethreshold); @@ -114,18 +151,42 @@ TEST_F(FIREtest, SecondHalf) mdrun->first_half(GlobalV::ofs_running); mdrun->second_half(); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00045447059554315662, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00032646833232493271, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.215709523063016e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.0005213674681407162, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00059486888444406608, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00035886062145122004, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00052406920303529794, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, 4.8706739346586155e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00020129054406946794, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00045717233044145918, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00021969381277291936, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00010541298215149392, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00045447059554315662, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00032646833232493271, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.215709523063016e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.0005213674681407162, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00059486888444406608, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00035886062145122004, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00052406920303529794, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, 4.8706739346586155e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00020129054406946794, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00045717233044145918, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00021969381277291936, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00010541298215149392, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -0.00010978976887416819, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 7.9202349471957007e-05, doublethreshold); @@ -169,7 +230,7 @@ TEST_F(FIREtest, Restart) mdrun->restart(readin_dir); remove("Restart_md.txt"); - FIRE* fire = dynamic_cast(mdrun); + FIRE *fire = dynamic_cast(mdrun); EXPECT_EQ(mdrun->step_rst_, 3); EXPECT_EQ(fire->get_alpha(), 0.1); EXPECT_EQ(fire->get_negative_count(), 0); diff --git a/source/source_md/test/langevin_test.cpp b/source/source_md/test/langevin_test.cpp index 784b6a9fd1..601aad1e9f 100644 --- a/source/source_md/test/langevin_test.cpp +++ b/source/source_md/test/langevin_test.cpp @@ -1,9 +1,11 @@ +#include "source_md/langevin.h" + +#include "setcell.h" #include "source_cell/domain_decomposition.h" +#include "source_esolver/esolver_lj.h" + #include "gmock/gmock.h" #include "gtest/gtest.h" -#include "setcell.h" -#include "source_esolver/esolver_lj.h" -#include "source_md/langevin.h" #define doublethreshold 1e-12 /************************************************ @@ -34,12 +36,12 @@ class Langevin_test : public testing::Test { protected: - MD_base* mdrun; + MD_base *mdrun; UnitCell ucell; MDCell mdcell; DomainDecomposition decomp; Input_para inp; - ModuleESolver::ESolver* p_esolver; + ModuleESolver::ESolver *p_esolver; /// the md writers take the directory explicitly, so the test owns it const std::string out_dir = "./"; const std::string readin_dir = "./"; @@ -51,9 +53,19 @@ class Langevin_test : public testing::Test p_esolver = new ModuleESolver::ESolver_LJ(); mdcell = Setcell::setup_mdcell(ucell); - decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0); + decomp.init( + ModuleBase::world_comm_domain(), + mdcell.latvec(), + mdcell.lat0(), + 0.0, + 0.0); p_esolver->before_all_runners(mdcell, inp); - mdrun = new Langevin(inp.mdp, inp.cal_stress, inp.init_vel, /*my_rank=*/0, mdcell); + mdrun = new Langevin( + inp.mdp, + inp.cal_stress, + inp.init_vel, + /*my_rank=*/0, + mdcell); mdrun->setup(p_esolver, readin_dir, decomp); } @@ -82,18 +94,42 @@ TEST_F(Langevin_test, first_half) { mdrun->first_half(GlobalV::ofs_running); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, 0.00012104549072633688, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 2.6272991724490339e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, 0.0002984728051383459, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00066077703137157329, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, 0.00017245549939737259, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, -0.00015046260270490386, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00046755850510571406, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, 0.00030490494761200812, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00036886672854369307, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00029624371924650601, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00013444493932002199, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, 9.0496812405138627e-05, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, 0.00012104549072633688, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 2.6272991724490339e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, 0.0002984728051383459, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00066077703137157329, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, 0.00017245549939737259, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, -0.00015046260270490386, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00046755850510571406, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, 0.00030490494761200812, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00036886672854369307, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00029624371924650601, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00013444493932002199, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, 9.0496812405138627e-05, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, 2.9279504031204254e-05, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 6.3551327892765255e-06, doublethreshold); @@ -115,18 +151,42 @@ TEST_F(Langevin_test, second_half) mdrun->second_half(); ; - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, 0.00012104549072633688, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 2.6272991724490339e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, 0.0002984728051383459, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00066077703137157329, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, 0.00017245549939737259, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, -0.00015046260270490386, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00046755850510571406, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, 0.00030490494761200812, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00036886672854369307, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00029624371924650601, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00013444493932002199, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, 9.0496812405138627e-05, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, 0.00012104549072633688, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 2.6272991724490339e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, 0.0002984728051383459, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00066077703137157329, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, 0.00017245549939737259, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, -0.00015046260270490386, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00046755850510571406, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, 0.00030490494761200812, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00036886672854369307, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00029624371924650601, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00013444493932002199, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, 9.0496812405138627e-05, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -2.3049731761587064e-05, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 7.1603385162874621e-06, doublethreshold); diff --git a/source/source_md/test/lj_pot_test.cpp b/source/source_md/test/lj_pot_test.cpp index 0ae49fb2b4..ce712e2708 100644 --- a/source/source_md/test/lj_pot_test.cpp +++ b/source/source_md/test/lj_pot_test.cpp @@ -1,9 +1,10 @@ -#include "source_cell/domain_decomposition.h" -#include "gtest/gtest.h" #include "setcell.h" +#include "source_cell/domain_decomposition.h" #include "source_esolver/esolver_lj.h" #include "source_io/module_parameter/parameter.h" #include "source_md/md_func.h" + +#include "gtest/gtest.h" #define doublethreshold 1e-12 /************************************************ @@ -35,26 +36,37 @@ class LJ_pot_test : public testing::Test Setcell::parameters(inp); } - void TearDown() - { - } + void TearDown() {} }; TEST_F(LJ_pot_test, potential) { - ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); + ModuleESolver::ESolver *p_esolver = new ModuleESolver::ESolver_LJ(); MDCell mdcell; mdcell = Setcell::setup_mdcell(ucell); DomainDecomposition decomp; - decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0); + decomp.init( + ModuleBase::world_comm_domain(), + mdcell.latvec(), + mdcell.lat0(), + 0.0, + 0.0); p_esolver->before_all_runners(mdcell, inp); - MD_func::force_virial(p_esolver, 0, mdcell, decomp, potential, true, stress, false); + MD_func::force_virial( + p_esolver, + 0, + mdcell, + decomp, + potential, + true, + stress, + false); EXPECT_NEAR(potential, -0.011957818623534381, doublethreshold); } TEST_F(LJ_pot_test, unitcell_compatibility) { - ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); + ModuleESolver::ESolver *p_esolver = new ModuleESolver::ESolver_LJ(); p_esolver->before_all_runners(ucell, inp); p_esolver->runner(ucell, 0); p_esolver->cal_force(ucell, stress); @@ -67,14 +79,27 @@ TEST_F(LJ_pot_test, unitcell_compatibility) TEST_F(LJ_pot_test, force) { - ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); + ModuleESolver::ESolver *p_esolver = new ModuleESolver::ESolver_LJ(); MDCell mdcell; mdcell = Setcell::setup_mdcell(ucell); DomainDecomposition decomp; - decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0); + decomp.init( + ModuleBase::world_comm_domain(), + mdcell.latvec(), + mdcell.lat0(), + 0.0, + 0.0); p_esolver->before_all_runners(mdcell, inp); - MD_func::force_virial(p_esolver, 0, mdcell, decomp, potential, true, stress, false); - const std::vector& atoms = mdcell.owned_atoms(); + MD_func::force_virial( + p_esolver, + 0, + mdcell, + decomp, + potential, + true, + stress, + false); + const std::vector &atoms = mdcell.owned_atoms(); EXPECT_NEAR(atoms[0].force.x, 0.00049817733089377704, doublethreshold); EXPECT_NEAR(atoms[0].force.y, 0.00082237246837022328, doublethreshold); EXPECT_NEAR(atoms[0].force.z, 0.0, doublethreshold); @@ -95,7 +120,12 @@ TEST_F(LJ_pot_test, mdcell_cal_force) MDCell mdcell; mdcell = Setcell::setup_mdcell(ucell); DomainDecomposition decomp; - decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0); + decomp.init( + ModuleBase::world_comm_domain(), + mdcell.latvec(), + mdcell.lat0(), + 0.0, + 0.0); p_esolver.before_all_runners(mdcell, inp); decomp.prepare_neighbors(mdcell); p_esolver.runner(mdcell, 0); @@ -104,7 +134,7 @@ TEST_F(LJ_pot_test, mdcell_cal_force) p_esolver.cal_force(mdcell, force); for (int iat = 0; iat < mdcell.owned_atoms().size(); ++iat) { - const LocalAtom& atom = mdcell.owned_atoms()[static_cast(iat)]; + const LocalAtom &atom = mdcell.owned_atoms()[static_cast(iat)]; EXPECT_DOUBLE_EQ(force(iat, 0), atom.force.x); EXPECT_DOUBLE_EQ(force(iat, 1), atom.force.y); EXPECT_DOUBLE_EQ(force(iat, 2), atom.force.z); @@ -113,13 +143,26 @@ TEST_F(LJ_pot_test, mdcell_cal_force) TEST_F(LJ_pot_test, stress) { - ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); + ModuleESolver::ESolver *p_esolver = new ModuleESolver::ESolver_LJ(); MDCell mdcell; mdcell = Setcell::setup_mdcell(ucell); DomainDecomposition decomp; - decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0); + decomp.init( + ModuleBase::world_comm_domain(), + mdcell.latvec(), + mdcell.lat0(), + 0.0, + 0.0); p_esolver->before_all_runners(mdcell, inp); - MD_func::force_virial(p_esolver, 0, mdcell, decomp, potential, true, stress, false); + MD_func::force_virial( + p_esolver, + 0, + mdcell, + decomp, + potential, + true, + stress, + false); EXPECT_NEAR(stress(0, 0), 8.0360222227631859e-07, doublethreshold); EXPECT_NEAR(stress(0, 1), 1.7207745586539077e-07, doublethreshold); EXPECT_NEAR(stress(0, 2), 0.0, doublethreshold); @@ -137,7 +180,12 @@ TEST_F(LJ_pot_test, mdcell_stress_includes_external_pressure) MDCell mdcell; mdcell = Setcell::setup_mdcell(ucell); DomainDecomposition decomp; - decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0); + decomp.init( + ModuleBase::world_comm_domain(), + mdcell.latvec(), + mdcell.lat0(), + 0.0, + 0.0); Input_para input = inp; p_esolver.before_all_runners(mdcell, input); decomp.prepare_neighbors(mdcell); @@ -164,7 +212,7 @@ TEST_F(LJ_pot_test, mdcell_stress_includes_external_pressure) TEST_F(LJ_pot_test, RcutSearchRadius) { - ModuleESolver::ESolver_LJ* p_esolver = new ModuleESolver::ESolver_LJ(); + ModuleESolver::ESolver_LJ *p_esolver = new ModuleESolver::ESolver_LJ(); ucell.ntype = 2; std::vector rcut = {3.0}; p_esolver->rcut_search_radius_for_testing(ucell.ntype, rcut); @@ -189,7 +237,7 @@ TEST_F(LJ_pot_test, RcutSearchRadius) TEST_F(LJ_pot_test, SetC6C12) { - ModuleESolver::ESolver_LJ* p_esolver = new ModuleESolver::ESolver_LJ(); + ModuleESolver::ESolver_LJ *p_esolver = new ModuleESolver::ESolver_LJ(); ucell.ntype = 2; // no rule @@ -197,7 +245,11 @@ TEST_F(LJ_pot_test, SetC6C12) std::vector lj_epsilon = {0.1, 0.2, 0.3}; std::vector lj_sigma = {0.2, 0.4, 0.6}; - p_esolver->set_c6_c12_for_testing(ucell.ntype, rule, lj_epsilon, lj_sigma); + p_esolver->set_c6_c12_for_testing( + ucell.ntype, + rule, + lj_epsilon, + lj_sigma); for (int i = 0; i < ucell.ntype; i++) { @@ -205,7 +257,8 @@ TEST_F(LJ_pot_test, SetC6C12) { int k = i * (i + 1) / 2 + j; double temp = pow(lj_sigma[k] * ModuleBase::ANGSTROM_AU, 6); - EXPECT_NEAR(p_esolver->get_lj_c6()(i, j), 4.0 * lj_epsilon[k] * temp / ModuleBase::Ry_to_eV, doublethreshold); + EXPECT_NEAR( + p_esolver->get_lj_c6()(i, j), 4.0 * lj_epsilon[k] * temp / ModuleBase::Ry_to_eV, doublethreshold); EXPECT_NEAR(p_esolver->get_lj_c12()(i, j), p_esolver->get_lj_c6()(i, j) * temp, doublethreshold); EXPECT_DOUBLE_EQ(p_esolver->get_lj_c6()(i, j), p_esolver->get_lj_c6()(j, i)); EXPECT_DOUBLE_EQ(p_esolver->get_lj_c12()(i, j), p_esolver->get_lj_c12()(j, i)); @@ -217,7 +270,11 @@ TEST_F(LJ_pot_test, SetC6C12) lj_epsilon = {0.1, 0.2}; lj_sigma = {0.2, 0.4}; - p_esolver->set_c6_c12_for_testing(ucell.ntype, rule, lj_epsilon, lj_sigma); + p_esolver->set_c6_c12_for_testing( + ucell.ntype, + rule, + lj_epsilon, + lj_sigma); for (int i = 0; i < ucell.ntype; i++) { @@ -227,12 +284,12 @@ TEST_F(LJ_pot_test, SetC6C12) for (int j = 0; j < i; j++) { - EXPECT_NEAR(p_esolver->get_lj_c6()(i, j), - std::sqrt(p_esolver->get_lj_c6()(i, i) * p_esolver->get_lj_c6()(j, j)), - doublethreshold); - EXPECT_NEAR(p_esolver->get_lj_c12()(i, j), - std::sqrt(p_esolver->get_lj_c12()(i, i) * p_esolver->get_lj_c12()(j, j)), - doublethreshold); + EXPECT_NEAR( + p_esolver->get_lj_c6()(i, j), std::sqrt(p_esolver->get_lj_c6()(i, i) * p_esolver->get_lj_c6()(j, j)), + doublethreshold); + EXPECT_NEAR( + p_esolver->get_lj_c12()(i, j), std::sqrt(p_esolver->get_lj_c12()(i, i) * p_esolver->get_lj_c12()(j, j)), + doublethreshold); EXPECT_DOUBLE_EQ(p_esolver->get_lj_c6()(i, j), p_esolver->get_lj_c6()(j, i)); EXPECT_DOUBLE_EQ(p_esolver->get_lj_c12()(i, j), p_esolver->get_lj_c12()(j, i)); } @@ -243,16 +300,20 @@ TEST_F(LJ_pot_test, SetC6C12) lj_epsilon = {0.1, 0.2}; lj_sigma = {0.2, 0.4}; - p_esolver->set_c6_c12_for_testing(ucell.ntype, rule, lj_epsilon, lj_sigma); + p_esolver->set_c6_c12_for_testing( + ucell.ntype, + rule, + lj_epsilon, + lj_sigma); for (int i = 0; i < ucell.ntype; i++) { for (int j = 0; j <= i; j++) { double temp = pow((lj_sigma[i] + lj_sigma[j]) / 2 * ModuleBase::ANGSTROM_AU, 6); - EXPECT_NEAR(p_esolver->get_lj_c6()(i, j), - 4.0 * std::sqrt(lj_epsilon[i] * lj_epsilon[j]) * temp / ModuleBase::Ry_to_eV, - doublethreshold); + EXPECT_NEAR( + p_esolver->get_lj_c6()(i, j), + 4.0 * std::sqrt(lj_epsilon[i] * lj_epsilon[j]) * temp / ModuleBase::Ry_to_eV, doublethreshold); EXPECT_NEAR(p_esolver->get_lj_c12()(i, j), p_esolver->get_lj_c6()(i, j) * temp, doublethreshold); EXPECT_DOUBLE_EQ(p_esolver->get_lj_c6()(i, j), p_esolver->get_lj_c6()(j, i)); EXPECT_DOUBLE_EQ(p_esolver->get_lj_c12()(i, j), p_esolver->get_lj_c12()(j, i)); @@ -262,7 +323,7 @@ TEST_F(LJ_pot_test, SetC6C12) TEST_F(LJ_pot_test, CalEnShift) { - ModuleESolver::ESolver_LJ* p_esolver = new ModuleESolver::ESolver_LJ(); + ModuleESolver::ESolver_LJ *p_esolver = new ModuleESolver::ESolver_LJ(); ucell.ntype = 2; std::vector rcut = {3.0}; @@ -271,7 +332,11 @@ TEST_F(LJ_pot_test, CalEnShift) int rule = 1; std::vector lj_epsilon = {0.1, 0.2, 0.3}; std::vector lj_sigma = {0.2, 0.4, 0.6}; - p_esolver->set_c6_c12_for_testing(ucell.ntype, rule, lj_epsilon, lj_sigma); + p_esolver->set_c6_c12_for_testing( + ucell.ntype, + rule, + lj_epsilon, + lj_sigma); // false p_esolver->cal_en_shift_for_testing(ucell.ntype, false); diff --git a/source/source_md/test/md_func_test.cpp b/source/source_md/test/md_func_test.cpp index 7eeccd7869..adf1423f3c 100644 --- a/source/source_md/test/md_func_test.cpp +++ b/source/source_md/test/md_func_test.cpp @@ -1,17 +1,24 @@ +#include "source_md/md_func.h" + +#include "md_test_fixture.h" #include "source_cell/domain_decomposition.h" + #include "gmock/gmock.h" #include "gtest/gtest.h" -#include "md_test_fixture.h" -#include "source_md/md_func.h" #define doublethreshold 1e-12 namespace { -void initialize_mdcell_from_ucell(MDCell& mdcell, UnitCell& ucell, DomainDecomposition& decomp) +void initialize_mdcell_from_ucell(MDCell &mdcell, UnitCell &ucell, DomainDecomposition &decomp) { const ModuleBase::CommunicationDomain comm_domain = ModuleBase::world_comm_domain(); - decomp.init(comm_domain, ucell.latvec, ucell.lat0, 0.0, 0.0); + decomp.init( + comm_domain, + ucell.latvec, + ucell.lat0, + 0.0, + 0.0); const std::vector owned_atoms = decomp.split_owned_atoms_from_ucell(ucell); std::vector type_labels; std::vector type_masses; @@ -22,17 +29,18 @@ void initialize_mdcell_from_ucell(MDCell& mdcell, UnitCell& ucell, DomainDecompo type_masses.push_back(ucell.atoms[it].mass); type_atom_counts.push_back(ucell.atoms[it].na); } - mdcell.initialize_from_owned_atoms(ucell.latvec, - ucell.GT, - ucell.lat0, - ucell.omega, - ucell.nat, - owned_atoms, - type_labels, - type_masses, - type_atom_counts, - 0.0, - comm_domain); + mdcell.initialize_from_owned_atoms( + ucell.latvec, + ucell.GT, + ucell.lat0, + ucell.omega, + ucell.nat, + owned_atoms, + type_labels, + type_masses, + type_atom_counts, + 0.0, + comm_domain); mdcell.set_backing_unitcell(ucell); } } // namespace @@ -114,7 +122,12 @@ TEST_F(MD_func_test, RescaleVel) vel[i].z = 0.3; } temperature = 300.0 / ModuleBase::Hartree_to_K; - MD_func::rescale_vel(natom, temperature, allmass, frozen_freedom, vel); + MD_func::rescale_vel( + natom, + temperature, + allmass, + frozen_freedom, + vel); EXPECT_DOUBLE_EQ(vel[0].x, 4.579010775069125e-05); EXPECT_DOUBLE_EQ(vel[0].y, 9.15802155013825e-05); @@ -136,11 +149,15 @@ TEST_F(MD_func_test, compute_stress) MDCell mdcell; DomainDecomposition decomp; initialize_mdcell_from_ucell(mdcell, ucell, decomp); - for (LocalAtom& atom : mdcell.owned_atoms()) + for (LocalAtom &atom : mdcell.owned_atoms()) { atom.vel = test_velocity; } - MD_func::compute_stress(mdcell, true, virial, stress); + MD_func::compute_stress( + mdcell, + true, + virial, + stress); EXPECT_DOUBLE_EQ(stress(0, 0), 2.9128300667662788); EXPECT_DOUBLE_EQ(stress(0, 1), 5.8256601335325575); EXPECT_DOUBLE_EQ(stress(0, 2), 8.7384902002988341); @@ -157,11 +174,17 @@ TEST_F(MD_func_test, dump_info) MDCell mdcell; DomainDecomposition decomp; initialize_mdcell_from_ucell(mdcell, ucell, decomp); - for (LocalAtom& atom : mdcell.owned_atoms()) + for (LocalAtom &atom : mdcell.owned_atoms()) { atom.vel = ModuleBase::Vector3(0.0, 0.0, 0.0); } - MD_func::dump_info(0, out_dir, mdcell, inp.mdp, inp.cal_stress, virial); + MD_func::dump_info( + 0, + out_dir, + mdcell, + inp.mdp, + inp.cal_stress, + virial); std::ifstream ifs("MD_dump"); std::string output_str; getline(ifs, output_str); @@ -189,25 +212,35 @@ TEST_F(MD_func_test, dump_info) output_str, testing::HasSubstr("INDEX LABEL POSITION (Angstrom) FORCE (eV/Angstrom) VELOCITY (Angstrom/fs)")); getline(ifs, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" 0 Ar 0.000000000000 0.000000000000 0.000000000000 0.000000000000 " - "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); + EXPECT_THAT( + output_str, testing::HasSubstr( + " 0 Ar 0.000000000000 0.000000000000 0.000000000000 0.000000000000 " + "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); getline(ifs, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" 1 Ar 2.751720222021 2.751720222021 0.000000000000 0.000000000000 " - "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); + EXPECT_THAT( + output_str, testing::HasSubstr( + " 1 Ar 2.751720222021 2.751720222021 0.000000000000 0.000000000000 " + "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); getline(ifs, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" 2 Ar 2.698802525444 0.000000000000 2.645884828867 0.000000000000 " - "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); + EXPECT_THAT( + output_str, testing::HasSubstr( + " 2 Ar 2.698802525444 0.000000000000 2.645884828867 0.000000000000 " + "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); getline(ifs, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" 3 Ar 0.000000000000 2.804637918599 2.645884828867 0.000000000000 " - "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); + EXPECT_THAT( + output_str, testing::HasSubstr( + " 3 Ar 0.000000000000 2.804637918599 2.645884828867 0.000000000000 " + "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); ifs.close(); // append - MD_func::dump_info(1, out_dir, mdcell, inp.mdp, inp.cal_stress, virial); + MD_func::dump_info( + 1, + out_dir, + mdcell, + inp.mdp, + inp.cal_stress, + virial); std::ifstream ifs2("MD_dump"); getline(ifs2, output_str); EXPECT_THAT(output_str, testing::HasSubstr("MDSTEP: 0")); @@ -234,21 +267,25 @@ TEST_F(MD_func_test, dump_info) output_str, testing::HasSubstr("INDEX LABEL POSITION (Angstrom) FORCE (eV/Angstrom) VELOCITY (Angstrom/fs)")); getline(ifs2, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" 0 Ar 0.000000000000 0.000000000000 0.000000000000 0.000000000000 " - "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); + EXPECT_THAT( + output_str, testing::HasSubstr( + " 0 Ar 0.000000000000 0.000000000000 0.000000000000 0.000000000000 " + "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); getline(ifs2, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" 1 Ar 2.751720222021 2.751720222021 0.000000000000 0.000000000000 " - "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); + EXPECT_THAT( + output_str, testing::HasSubstr( + " 1 Ar 2.751720222021 2.751720222021 0.000000000000 0.000000000000 " + "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); getline(ifs2, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" 2 Ar 2.698802525444 0.000000000000 2.645884828867 0.000000000000 " - "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); + EXPECT_THAT( + output_str, testing::HasSubstr( + " 2 Ar 2.698802525444 0.000000000000 2.645884828867 0.000000000000 " + "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); getline(ifs2, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" 3 Ar 0.000000000000 2.804637918599 2.645884828867 0.000000000000 " - "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); + EXPECT_THAT( + output_str, testing::HasSubstr( + " 3 Ar 0.000000000000 2.804637918599 2.645884828867 0.000000000000 " + "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); getline(ifs2, output_str); getline(ifs2, output_str); getline(ifs2, output_str); @@ -276,21 +313,25 @@ TEST_F(MD_func_test, dump_info) output_str, testing::HasSubstr("INDEX LABEL POSITION (Angstrom) FORCE (eV/Angstrom) VELOCITY (Angstrom/fs)")); getline(ifs2, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" 0 Ar 0.000000000000 0.000000000000 0.000000000000 0.000000000000 " - "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); + EXPECT_THAT( + output_str, testing::HasSubstr( + " 0 Ar 0.000000000000 0.000000000000 0.000000000000 0.000000000000 " + "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); getline(ifs2, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" 1 Ar 2.751720222021 2.751720222021 0.000000000000 0.000000000000 " - "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); + EXPECT_THAT( + output_str, testing::HasSubstr( + " 1 Ar 2.751720222021 2.751720222021 0.000000000000 0.000000000000 " + "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); getline(ifs2, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" 2 Ar 2.698802525444 0.000000000000 2.645884828867 0.000000000000 " - "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); + EXPECT_THAT( + output_str, testing::HasSubstr( + " 2 Ar 2.698802525444 0.000000000000 2.645884828867 0.000000000000 " + "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); getline(ifs2, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" 3 Ar 0.000000000000 2.804637918599 2.645884828867 0.000000000000 " - "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); + EXPECT_THAT( + output_str, testing::HasSubstr( + " 3 Ar 0.000000000000 2.804637918599 2.645884828867 0.000000000000 " + "0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000")); ifs2.close(); remove("MD_dump"); @@ -309,27 +350,27 @@ TEST_F(MD_func_test, print_stress) EXPECT_THAT(output_str, testing::HasSubstr("IONIC (KINETIC) PART OF STRESS: 0 kbar")); getline(ifs, output_str); EXPECT_THAT(output_str, testing::HasSubstr("MD PRESSURE (ELECTRONS+IONS) : 0 kbar")); -/* - getline(ifs, output_str); - getline(ifs, output_str); - getline(ifs, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" ><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><")); - getline(ifs, output_str); - getline(ifs, output_str); - EXPECT_THAT(output_str, testing::HasSubstr(" MD STRESS (kbar)")); - getline(ifs, output_str); - getline(ifs, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr(" ><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><")); - getline(ifs, output_str); - getline(ifs, output_str); - EXPECT_THAT(output_str, testing::HasSubstr(" 0 0 0")); - getline(ifs, output_str); - EXPECT_THAT(output_str, testing::HasSubstr(" 0 0 0")); - getline(ifs, output_str); - EXPECT_THAT(output_str, testing::HasSubstr(" 0 0 0")); -*/ + /* + getline(ifs, output_str); + getline(ifs, output_str); + getline(ifs, output_str); + EXPECT_THAT(output_str, + testing::HasSubstr(" ><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><")); + getline(ifs, output_str); + getline(ifs, output_str); + EXPECT_THAT(output_str, testing::HasSubstr(" MD STRESS (kbar)")); + getline(ifs, output_str); + getline(ifs, output_str); + EXPECT_THAT(output_str, + testing::HasSubstr(" ><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><")); + getline(ifs, output_str); + getline(ifs, output_str); + EXPECT_THAT(output_str, testing::HasSubstr(" 0 0 0")); + getline(ifs, output_str); + EXPECT_THAT(output_str, testing::HasSubstr(" 0 0 0")); + getline(ifs, output_str); + EXPECT_THAT(output_str, testing::HasSubstr(" 0 0 0")); + */ ifs.close(); remove("running.log"); @@ -346,7 +387,11 @@ TEST_F(MD_func_test, current_md_info_mdcell_accepts_step_only_restart) initialize_mdcell_from_ucell(mdcell, ucell, decomp); int istep = -1; double temperature = 0.0; - MD_func::current_md_info(mdcell, "./", istep, temperature); + MD_func::current_md_info( + mdcell, + "./", + istep, + temperature); EXPECT_EQ(istep, 123); EXPECT_DOUBLE_EQ(temperature, 0.0); @@ -360,7 +405,7 @@ TEST_F(MD_func_test, global_dof_mdcell) initialize_mdcell_from_ucell(mdcell, ucell, decomp); EXPECT_EQ(MD_func::global_dof(mdcell), 9); - for (LocalAtom& atom : mdcell.owned_atoms()) + for (LocalAtom &atom : mdcell.owned_atoms()) { atom.mbl.x = 0; } @@ -374,5 +419,11 @@ TEST_F(MD_func_test, current_step_warning) initialize_mdcell_from_ucell(mdcell, ucell, decomp); int istep = 0; double temperature = 0.0; - EXPECT_EXIT(MD_func::current_md_info(mdcell, "./", istep, temperature), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT( + MD_func::current_md_info( + mdcell, + "./", + istep, + temperature), + ::testing::ExitedWithCode(1), ""); } diff --git a/source/source_md/test/md_test_fixture.h b/source/source_md/test/md_test_fixture.h index 853ccd1a3f..ed7df1726a 100644 --- a/source/source_md/test/md_test_fixture.h +++ b/source/source_md/test/md_test_fixture.h @@ -1,13 +1,13 @@ #ifndef MD_TEST_FIXTURE_H #define MD_TEST_FIXTURE_H -#include "gtest/gtest.h" +#include "setcell.h" +#include "source_cell/domain_decomposition.h" #include "source_esolver/esolver_lj.h" #include "source_io/module_parameter/input_parameter.h" #include "source_md/md_base.h" -#include "source_cell/domain_decomposition.h" -#include "setcell.h" +#include "gtest/gtest.h" #include #include #include @@ -29,8 +29,7 @@ class MdTestBase : public testing::Test } }; -template -class MdIntegratorFixture : public MdTestBase +template class MdIntegratorFixture : public MdTestBase { protected: MDCell mdcell; @@ -43,9 +42,19 @@ class MdIntegratorFixture : public MdTestBase { MdTestBase::SetUp(); mdcell = Setcell::setup_mdcell(ucell); - decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0); + decomp.init( + ModuleBase::world_comm_domain(), + mdcell.latvec(), + mdcell.lat0(), + 0.0, + 0.0); p_esolver->before_all_runners(mdcell, inp); - mdrun.reset(new Integrator(this->inp.mdp, this->inp.cal_stress, this->inp.init_vel, /*my_rank=*/0, mdcell)); + mdrun.reset(new Integrator( + this->inp.mdp, + this->inp.cal_stress, + this->inp.init_vel, + /*my_rank=*/0, + mdcell)); mdrun->setup(p_esolver.get(), readin_dir, decomp); } }; @@ -59,11 +68,11 @@ class MdFuncTestFixture : public testing::Test std::vector> vel_store; std::vector> ionmbl_store; std::vector> force_store; - double* allmass = nullptr; - ModuleBase::Vector3* pos = nullptr; - ModuleBase::Vector3* vel = nullptr; - ModuleBase::Vector3* ionmbl = nullptr; - ModuleBase::Vector3* force = nullptr; + double *allmass = nullptr; + ModuleBase::Vector3 *pos = nullptr; + ModuleBase::Vector3 *vel = nullptr; + ModuleBase::Vector3 *ionmbl = nullptr; + ModuleBase::Vector3 *force = nullptr; ModuleBase::matrix virial; ModuleBase::matrix stress; double potential = 0.0; @@ -99,7 +108,7 @@ class LjPotTestFixture : public testing::Test { protected: std::vector> force_store; - ModuleBase::Vector3* force = nullptr; + ModuleBase::Vector3 *force = nullptr; ModuleBase::matrix stress; double potential = 0.0; int natom = 0; diff --git a/source/source_md/test/msst_test.cpp b/source/source_md/test/msst_test.cpp index 6a7e7cd3d2..65983593f3 100644 --- a/source/source_md/test/msst_test.cpp +++ b/source/source_md/test/msst_test.cpp @@ -1,9 +1,11 @@ +#include "source_md/msst.h" + +#include "setcell.h" #include "source_cell/domain_decomposition.h" +#include "source_esolver/esolver_lj.h" + #include "gmock/gmock.h" #include "gtest/gtest.h" -#include "setcell.h" -#include "source_esolver/esolver_lj.h" -#include "source_md/msst.h" #define doublethreshold 1e-12 /************************************************ @@ -34,12 +36,12 @@ class MSST_test : public testing::Test { protected: - MD_base* mdrun; + MD_base *mdrun; UnitCell ucell; MDCell mdcell; DomainDecomposition decomp; Input_para inp; - ModuleESolver::ESolver* p_esolver; + ModuleESolver::ESolver *p_esolver; /// the md writers take the directory explicitly, so the test owns it const std::string out_dir = "./"; const std::string readin_dir = "./"; @@ -51,9 +53,19 @@ class MSST_test : public testing::Test p_esolver = new ModuleESolver::ESolver_LJ(); mdcell = Setcell::setup_mdcell(ucell); - decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0); + decomp.init( + ModuleBase::world_comm_domain(), + mdcell.latvec(), + mdcell.lat0(), + 0.0, + 0.0); p_esolver->before_all_runners(mdcell, inp); - mdrun = new MSST(inp.mdp, inp.cal_stress, inp.init_vel, /*my_rank=*/0, mdcell); + mdrun = new MSST( + inp.mdp, + inp.cal_stress, + inp.init_vel, + /*my_rank=*/0, + mdcell); mdrun->setup(p_esolver, readin_dir, decomp); } @@ -107,18 +119,42 @@ TEST_F(MSST_test, first_half) EXPECT_NEAR(mdcell.latvec().e33, 9.9959581179144905, doublethreshold); EXPECT_NEAR(mdcell.omega(), 999.59581179144902, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054271823071484467, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029442816868202821, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7685149290774873e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00062875654254500096, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060353746208327032, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.0003968957326219519, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055074716824834991, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1576283073263842e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022262503373940808, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046470885642230719, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032068557647491725, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011658554959218052, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054271823071484467, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029442816868202821, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7685149290774873e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00062875654254500096, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060353746208327032, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.0003968957326219519, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055074716824834991, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1576283073263842e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022262503373940808, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046470885642230719, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032068557647491725, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011658554959218052, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -0.00013127726219846624, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 7.121876825065284e-05, doublethreshold); @@ -153,18 +189,42 @@ TEST_F(MSST_test, second_half) EXPECT_NEAR(mdcell.latvec().e33, 9.9959581179144905, doublethreshold); EXPECT_NEAR(mdcell.omega(), 999.59581179144902, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054271823071484467, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029442816868202821, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7685149290774873e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00062875654254500096, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060353746208327032, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.0003968957326219519, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055074716824834991, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1576283073263842e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022262503373940808, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046470885642230719, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032068557647491725, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011658554959218052, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054271823071484467, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029442816868202821, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7685149290774873e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00062875654254500096, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060353746208327032, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.0003968957326219519, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055074716824834991, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1576283073263842e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022262503373940808, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046470885642230719, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032068557647491725, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011658554959218052, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -0.00013113585103096098, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 7.1452204506509308e-05, doublethreshold); @@ -210,7 +270,7 @@ TEST_F(MSST_test, restart) mdrun->restart(readin_dir); remove("Restart_md.txt"); - MSST* msst = dynamic_cast(mdrun); + MSST *msst = dynamic_cast(mdrun); EXPECT_EQ(mdrun->step_rst_, 3); EXPECT_EQ(msst->get_omega()[inp.mdp.msst_direction], -0.00977662); EXPECT_EQ(msst->get_e0(), -0.00768262); @@ -230,13 +290,15 @@ TEST_F(MSST_test, print_md) getline(ifs, output_str); EXPECT_THAT(output_str, testing::HasSubstr(" ELECTRONIC PART OF STRESS: 0.24609992 kbar")); getline(ifs, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr( - " IONIC (KINETIC) PART OF STRESS: 0.8301538 kbar")); // result different from other MD methods + EXPECT_THAT( + output_str, + testing::HasSubstr( + " IONIC (KINETIC) PART OF STRESS: 0.8301538 kbar")); // result different from other MD methods getline(ifs, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr( - " MD PRESSURE (ELECTRONS+IONS) : 1.0762537 kbar")); // result different from other MD methods + EXPECT_THAT( + output_str, + testing::HasSubstr( + " MD PRESSURE (ELECTRONS+IONS) : 1.0762537 kbar")); // result different from other MD methods getline(ifs, output_str); getline(ifs, output_str); EXPECT_THAT( diff --git a/source/source_md/test/nhchain_test.cpp b/source/source_md/test/nhchain_test.cpp index a73b9811c4..a26597cd66 100644 --- a/source/source_md/test/nhchain_test.cpp +++ b/source/source_md/test/nhchain_test.cpp @@ -1,9 +1,11 @@ +#include "source_md/nhchain.h" + +#include "setcell.h" #include "source_cell/domain_decomposition.h" +#include "source_esolver/esolver_lj.h" + #include "gmock/gmock.h" #include "gtest/gtest.h" -#include "setcell.h" -#include "source_esolver/esolver_lj.h" -#include "source_md/nhchain.h" #define doublethreshold 1e-12 /************************************************ * unit test of functions in nhchain.h @@ -32,12 +34,12 @@ class NHC_test : public testing::Test { protected: - MD_base* mdrun; + MD_base *mdrun; UnitCell ucell; MDCell mdcell; DomainDecomposition decomp; Input_para inp; - ModuleESolver::ESolver* p_esolver; + ModuleESolver::ESolver *p_esolver; /// the md writers take the directory explicitly, so the test owns it const std::string out_dir = "./"; const std::string readin_dir = "./"; @@ -53,9 +55,19 @@ class NHC_test : public testing::Test inp.mdp.md_pfirst = 1; inp.mdp.md_plast = 1; mdcell = Setcell::setup_mdcell(ucell); - decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0); + decomp.init( + ModuleBase::world_comm_domain(), + mdcell.latvec(), + mdcell.lat0(), + 0.0, + 0.0); p_esolver->before_all_runners(mdcell, inp); - mdrun = new Nose_Hoover(inp.mdp, inp.cal_stress, inp.init_vel, /*my_rank=*/0, mdcell); + mdrun = new Nose_Hoover( + inp.mdp, + inp.cal_stress, + inp.init_vel, + /*my_rank=*/0, + mdcell); mdrun->setup(p_esolver, readin_dir, decomp); } @@ -84,18 +96,42 @@ TEST_F(NHC_test, first_half) { mdrun->first_half(GlobalV::ofs_running); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00035596392702161582, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00026566987683715606, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -6.4082739615824722e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00037007414441809518, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00052501803299631633, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00044091358349508534, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00036876922955593201, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -2.6151466573228018e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00024731533582713971, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00035465901216238645, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00028549962273273618, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00012951550805257814, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00035596392702161582, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00026566987683715606, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -6.4082739615824722e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00037007414441809518, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00052501803299631633, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00044091358349508534, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00036876922955593201, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -2.6151466573228018e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00024731533582713971, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00035465901216238645, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00028549962273273618, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00012951550805257814, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -0.00010335325828338315, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 6.6973537793984337e-05, doublethreshold); @@ -116,18 +152,42 @@ TEST_F(NHC_test, second_half) mdrun->first_half(GlobalV::ofs_running); mdrun->second_half(); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00035596392702161582, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00026566987683715606, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -6.4082739615824722e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00037007414441809518, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00052501803299631633, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00044091358349508534, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00036876922955593201, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -2.6151466573228018e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00024731533582713971, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00035465901216238645, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00028549962273273618, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00012951550805257814, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00035596392702161582, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00026566987683715606, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -6.4082739615824722e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00037007414441809518, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00052501803299631633, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00044091358349508534, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00036876922955593201, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -2.6151466573228018e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00024731533582713971, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00035465901216238645, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00028549962273273618, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00012951550805257814, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -8.4972683205367143e-05, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 6.6834262571392232e-05, doublethreshold); @@ -165,8 +225,8 @@ TEST_F(NHC_test, write_restart) getline(ifs, output_str); EXPECT_THAT(output_str, testing::HasSubstr("-0.00658882 -0.0304055 -0.0188618 -0.0175663 ")); getline(ifs, output_str); - EXPECT_THAT(output_str, - testing::HasSubstr("0.583152 -0.106519 -0.895936 -0.634424 0.627532 -0.473422 ")); + EXPECT_THAT( + output_str, testing::HasSubstr("0.583152 -0.106519 -0.895936 -0.634424 0.627532 -0.473422 ")); getline(ifs, output_str); EXPECT_THAT(output_str, testing::HasSubstr("4")); getline(ifs, output_str); @@ -181,7 +241,7 @@ TEST_F(NHC_test, restart) mdrun->restart(readin_dir); remove("Restart_md.txt"); - Nose_Hoover* nhc = dynamic_cast(mdrun); + Nose_Hoover *nhc = dynamic_cast(mdrun); EXPECT_EQ(mdrun->step_rst_, 3); EXPECT_EQ(inp.mdp.md_tchain, 4); EXPECT_EQ(inp.mdp.md_pchain, 4); diff --git a/source/source_md/test/run_md_test.cpp b/source/source_md/test/run_md_test.cpp index 9c459d3e88..6eb8fcb676 100644 --- a/source/source_md/test/run_md_test.cpp +++ b/source/source_md/test/run_md_test.cpp @@ -1,9 +1,10 @@ -#include "source_cell/domain_decomposition.h" -#include "gtest/gtest.h" +#include "source_md/run_md.h" +#include "source_cell/domain_decomposition.h" #include "source_cell/mdcell.h" #include "source_cell/unitcell.h" -#include "source_md/run_md.h" + +#include "gtest/gtest.h" TEST(RunMDTest, prepare_mdcell_from_unitcell) { @@ -51,7 +52,7 @@ TEST(RunMDTest, prepare_mdcell_from_unitcell) ASSERT_EQ(mdcell.owned_atoms().size(), 1); EXPECT_TRUE(mdcell.has_backing_unitcell()); EXPECT_EQ(&mdcell.backing_unitcell(), &ucell); - const LocalAtom& atom = mdcell.owned_atoms()[0]; + const LocalAtom &atom = mdcell.owned_atoms()[0]; EXPECT_NEAR(atom.frac.x, 0.8, 1.0e-12); EXPECT_NEAR(atom.frac.y, 0.25, 1.0e-12); EXPECT_NEAR(atom.frac.z, 0.5, 1.0e-12); diff --git a/source/source_md/test/setcell.h b/source/source_md/test/setcell.h index 0cd0809c1a..b268ec5f95 100644 --- a/source/source_md/test/setcell.h +++ b/source/source_md/test/setcell.h @@ -4,12 +4,12 @@ #ifndef SETCELL_H #define SETCELL_H +#include "source_base/constants.h" +#include "source_base/parallel_cell.h" +#include "source_cell/mdcell.h" #include "source_cell/module_neighbor/sltk_atom_arrange.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" -#include "source_cell/mdcell.h" #include "source_cell/unitcell.h" -#include "source_base/constants.h" -#include "source_base/parallel_cell.h" #include "source_io/module_parameter/input_parameter.h" #include @@ -22,20 +22,18 @@ Magnetism::Magnetism() this->tot_mag = 0.0; this->abs_mag = 0.0; } -Magnetism::~Magnetism() -{ -} +Magnetism::~Magnetism() {} class Setcell { public: - static void setupcell(UnitCell& ucell) + static void setupcell(UnitCell &ucell) { ucell.ntype = 1; ucell.atoms = new Atom[ucell.ntype]; ucell.set_atom_flag = true; - + ucell.atoms[0].mass = 39.948; ucell.atoms[0].label = "Ar"; @@ -104,7 +102,7 @@ class Setcell /// from the global parameter singleton, and the output/readin /// directories are passed explicitly to setup(), write_restart() /// and restart(). - static void parameters(Input_para& input) + static void parameters(Input_para &input) { input.search_radius = 8.5 * ModuleBase::ANGSTROM_AU; input.esolver_type = "lj"; @@ -140,7 +138,7 @@ class Setcell input.mdp.md_tolerance = 0; }; - static MDCell setup_mdcell(UnitCell& ucell) + static MDCell setup_mdcell(UnitCell &ucell) { std::vector owned_atoms; std::vector type_labels; @@ -153,41 +151,41 @@ class Setcell type_atom_counts.push_back(ucell.atoms[it].na); for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { - owned_atoms.push_back(LocalAtom(ucell.atoms[it].tau[ia], - ucell.atoms[it].taud[ia], - ucell.atoms[it].vel[ia], - ModuleBase::Vector3(0.0, 0.0, 0.0), - ucell.atoms[it].mbl[ia], - ucell.atoms[it].mass / ModuleBase::AU_to_MASS, - it, - ia, - 0)); + owned_atoms.push_back(LocalAtom( + ucell.atoms[it].tau[ia], + ucell.atoms[it].taud[ia], + ucell.atoms[it].vel[ia], + ModuleBase::Vector3(0.0, 0.0, 0.0), + ucell.atoms[it].mbl[ia], + ucell.atoms[it].mass / ModuleBase::AU_to_MASS, + it, + ia, + 0)); } } MDCell mdcell; - mdcell.initialize_from_owned_atoms(ucell.latvec, - ucell.GT, - ucell.lat0, - ucell.omega, - ucell.nat, - owned_atoms, - type_labels, - type_masses, - type_atom_counts, - 0.0, - ModuleBase::world_comm_domain()); + mdcell.initialize_from_owned_atoms( + ucell.latvec, + ucell.GT, + ucell.lat0, + ucell.omega, + ucell.nat, + owned_atoms, + type_labels, + type_masses, + type_atom_counts, + 0.0, + ModuleBase::world_comm_domain()); return mdcell; } - static ModuleBase::Vector3 fractional_displacement(const LocalAtom& atom) + static ModuleBase::Vector3 fractional_displacement(const LocalAtom &atom) { - const ModuleBase::Vector3 initial_frac[] = { - ModuleBase::Vector3(0.0, 0.0, 0.0), + const ModuleBase::Vector3 initial_frac[] = {ModuleBase::Vector3(0.0, 0.0, 0.0), ModuleBase::Vector3(0.52, 0.52, 0.0), ModuleBase::Vector3(0.51, 0.0, 0.5), - ModuleBase::Vector3(0.0, 0.53, 0.5) - }; + ModuleBase::Vector3(0.0, 0.53, 0.5)}; ModuleBase::Vector3 displacement = atom.frac - initial_frac[atom.type_index]; displacement.x -= std::floor(displacement.x + 0.5); displacement.y -= std::floor(displacement.y + 0.5); diff --git a/source/source_md/test/test_md_base.cpp b/source/source_md/test/test_md_base.cpp index 2d0f69cfdb..b09e5fb33c 100644 --- a/source/source_md/test/test_md_base.cpp +++ b/source/source_md/test/test_md_base.cpp @@ -1,10 +1,9 @@ -#include "source_md/md_base.h" #include "source_base/constants.h" #include "source_base/parallel_cell.h" - -#include +#include "source_md/md_base.h" #include +#include #include namespace @@ -30,7 +29,18 @@ MDCell make_position_test_cell() const std::vector masses(1, 1.0); const std::vector counts(1, 2); MDCell cell; - cell.initialize_from_owned_atoms(lattice, inverse, 1.0, 1.0, 2, atoms, labels, masses, counts, 0.0, domain); + cell.initialize_from_owned_atoms( + lattice, + inverse, + 1.0, + 1.0, + 2, + atoms, + labels, + masses, + counts, + 0.0, + domain); return cell; } @@ -41,8 +51,13 @@ TEST(MdPositionTest, TinyNegativeDisplacementRemainsInsideHalfOpenCell) parameters.md_dt = ModuleBase::AU_to_FS; parameters.md_tfirst = 0.0; parameters.md_tlast = 0.0; - MD_base integrator(parameters, false, true, 0, cell); - LocalAtom& moving = cell.owned_atoms()[0]; + MD_base integrator( + parameters, + false, + true, + 0, + cell); + LocalAtom &moving = cell.owned_atoms()[0]; moving.frac.set(0.0, 0.0, 0.0); moving.cart.set(0.0, 0.0, 0.0); moving.vel.set(-1.0e-21, -2.0e-21, -3.0e-21); @@ -64,17 +79,22 @@ TEST(MdPositionTest, FiniteBoundaryCrossingsAndFrozenComponentsKeepTheirPosition parameters.md_dt = ModuleBase::AU_to_FS; parameters.md_tfirst = 0.0; parameters.md_tlast = 0.0; - MD_base integrator(parameters, false, true, 0, cell); + MD_base integrator( + parameters, + false, + true, + 0, + cell); cell.owned_atoms()[0].vel.set(-0.5, 1.75, -1.25); cell.owned_atoms()[1].mbl.set(0, 0, 0); cell.owned_atoms()[1].vel.set(-0.5, 1.75, -1.25); std::ofstream log; integrator.first_half(log); - const LocalAtom& moving = cell.owned_atoms()[0]; + const LocalAtom &moving = cell.owned_atoms()[0]; EXPECT_DOUBLE_EQ(moving.frac.x, 0.75); EXPECT_DOUBLE_EQ(moving.frac.y, 0.25); EXPECT_DOUBLE_EQ(moving.frac.z, 0.5); - const LocalAtom& frozen = cell.owned_atoms()[1]; + const LocalAtom &frozen = cell.owned_atoms()[1]; EXPECT_DOUBLE_EQ(frozen.frac.x, 0.25); EXPECT_DOUBLE_EQ(frozen.frac.y, 0.5); EXPECT_DOUBLE_EQ(frozen.frac.z, 0.75); diff --git a/source/source_md/test/test_md_func.cpp b/source/source_md/test/test_md_func.cpp index 9e0aef2f6e..2a2399e240 100644 --- a/source/source_md/test/test_md_func.cpp +++ b/source/source_md/test/test_md_func.cpp @@ -1,8 +1,8 @@ -#include "source_md/md_func.h" -#include "source_cell/mdcell.h" -#include "source_cell/domain_decomposition.h" #include "source_base/parallel_cell.h" +#include "source_cell/domain_decomposition.h" +#include "source_cell/mdcell.h" #include "source_esolver/esolver.h" +#include "source_md/md_func.h" #include @@ -10,30 +10,30 @@ namespace { class GhostForceSolver : public ModuleESolver::ESolver { -public: - void before_all_runners(BaseCell&, const Input_para&) override {} - void after_all_runners(BaseCell&) override {} - void runner(BaseCell& base, const int) override + public: + void before_all_runners(BaseCell &, const Input_para &) override {} + void after_all_runners(BaseCell &) override {} + void runner(BaseCell &base, const int) override { - MDCell& cell = static_cast(base); + MDCell &cell = static_cast(base); EXPECT_TRUE(cell.has_neighbor_search()); EXPECT_GT(cell.ghost_atoms().size(), 0); - for (LocalAtom& atom : cell.owned_atoms()) + for (LocalAtom &atom : cell.owned_atoms()) { atom.force.set(4.0, 0.0, 0.0); } - for (LocalAtom& atom : cell.ghost_atoms()) + for (LocalAtom &atom : cell.ghost_atoms()) { EXPECT_DOUBLE_EQ(atom.force.norm2(), 0.0); atom.force.set(2.0, 0.0, 0.0); } } double cal_energy() override { return 6.0; } - void cal_force(BaseCell&, ModuleBase::matrix&) override + void cal_force(BaseCell &, ModuleBase::matrix &) override { ADD_FAILURE() << "Direct MDCell forces must not use cal_force."; } - void cal_stress(BaseCell&, ModuleBase::matrix& stress) override + void cal_stress(BaseCell &, ModuleBase::matrix &stress) override { stress.zero_out(); stress(0, 0) = 8.0; @@ -52,24 +52,44 @@ TEST(MdForceVirialTest, ReturnsGhostForcesBeforeConvertingUnitsExactlyOnce) LocalAtom atom; atom.frac.set(0.05, 0.5, 0.5); atom.cart = atom.frac * lattice; - cell.initialize_from_owned_atoms(lattice, lattice.Inverse(), 1.0, 64.0, 1, - std::vector(1, atom), - std::vector(1, "X"), - std::vector(1, 1.0), - std::vector(1, 1), 0.4, domain); - decomp.init(domain, lattice, 1.0, 0.0, 0.4); + cell.initialize_from_owned_atoms( + lattice, + lattice.Inverse(), + 1.0, + 64.0, + 1, + std::vector(1, atom), + std::vector(1, "X"), + std::vector(1, 1.0), + std::vector(1, 1), + 0.4, + domain); + decomp.init( + domain, + lattice, + 1.0, + 0.0, + 0.4); cell.set_neighbor_cutoff(0.5); GhostForceSolver solver; double potential = 0.0; ModuleBase::matrix virial(3, 3); for (int step = 0; step < 2; ++step) { - MD_func::force_virial(&solver, step, cell, decomp, potential, true, virial, false); + MD_func::force_virial( + &solver, + step, + cell, + decomp, + potential, + true, + virial, + false); EXPECT_DOUBLE_EQ(potential, 3.0); EXPECT_DOUBLE_EQ(virial(0, 0), 4.0); ASSERT_EQ(cell.owned_atoms().size(), 1); EXPECT_DOUBLE_EQ(cell.owned_atoms()[0].force.x, 2.0 + cell.ghost_atoms().size()); - for (const LocalAtom& ghost : cell.ghost_atoms()) + for (const LocalAtom &ghost : cell.ghost_atoms()) { EXPECT_DOUBLE_EQ(ghost.force.norm2(), 0.0); } diff --git a/source/source_md/test/verlet_test.cpp b/source/source_md/test/verlet_test.cpp index 8d960f04d3..37d7bb5e8d 100644 --- a/source/source_md/test/verlet_test.cpp +++ b/source/source_md/test/verlet_test.cpp @@ -1,10 +1,11 @@ -#include "source_cell/domain_decomposition.h" -#include "gmock/gmock.h" -#include "gtest/gtest.h" +#include "source_md/verlet.h" + #include "setcell.h" +#include "source_cell/domain_decomposition.h" #include "source_esolver/esolver_lj.h" -#include "source_md/verlet.h" +#include "gmock/gmock.h" +#include "gtest/gtest.h" #include #define doublethreshold 1e-12 @@ -36,12 +37,12 @@ class Verlet_test : public testing::Test { protected: - MD_base* mdrun; + MD_base *mdrun; UnitCell ucell; MDCell mdcell; DomainDecomposition decomp; Input_para inp; - ModuleESolver::ESolver* p_esolver; + ModuleESolver::ESolver *p_esolver; /// the md writers take the directory explicitly, so the test owns it const std::string out_dir = "./"; const std::string readin_dir = "./"; @@ -53,9 +54,19 @@ class Verlet_test : public testing::Test p_esolver = new ModuleESolver::ESolver_LJ(); mdcell = Setcell::setup_mdcell(ucell); - decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0); + decomp.init( + ModuleBase::world_comm_domain(), + mdcell.latvec(), + mdcell.lat0(), + 0.0, + 0.0); p_esolver->before_all_runners(mdcell, inp); - mdrun = new Verlet(inp.mdp, inp.cal_stress, inp.init_vel, /*my_rank=*/0, mdcell); + mdrun = new Verlet( + inp.mdp, + inp.cal_stress, + inp.init_vel, + /*my_rank=*/0, + mdcell); mdrun->setup(p_esolver, readin_dir, decomp); } @@ -84,18 +95,42 @@ TEST_F(Verlet_test, first_half) { mdrun->first_half(GlobalV::ofs_running); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00063192793031220879, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060657401578200095, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00039873402383468892, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055351963726126224, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1639385612741475e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022365616007718661, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046704699702541431, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032230681977380224, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011712553572388214, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00063192793031220879, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060657401578200095, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00039873402383468892, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055351963726126224, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1639385612741475e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022365616007718661, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046704699702541431, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032230681977380224, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011712553572388214, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -0.00013193932519649473, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 7.1576379239356465e-05, doublethreshold); @@ -118,18 +153,42 @@ TEST_F(Verlet_test, NVE) mdrun->second_half(); ; - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00063192793031220879, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060657401578200095, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00039873402383468892, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055351963726126224, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1639385612741475e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022365616007718661, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046704699702541431, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032230681977380224, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011712553572388214, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00063192793031220879, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060657401578200095, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00039873402383468892, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055351963726126224, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1639385612741475e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022365616007718661, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046704699702541431, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032230681977380224, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011712553572388214, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -0.00013179791402898947, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 7.1809815495212933e-05, doublethreshold); @@ -153,18 +212,42 @@ TEST_F(Verlet_test, Anderson) mdrun->second_half(); ; - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00063192793031220879, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060657401578200095, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00039873402383468892, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055351963726126224, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1639385612741475e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022365616007718661, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046704699702541431, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032230681977380224, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011712553572388214, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00063192793031220879, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060657401578200095, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00039873402383468892, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055351963726126224, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1639385612741475e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022365616007718661, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046704699702541431, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032230681977380224, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011712553572388214, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -0.00013179791402898947, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 7.1809815495212933e-05, doublethreshold); @@ -188,18 +271,42 @@ TEST_F(Verlet_test, Berendsen) mdrun->second_half(); ; - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00063192793031220879, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060657401578200095, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00039873402383468892, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055351963726126224, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1639385612741475e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022365616007718661, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046704699702541431, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032230681977380224, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011712553572388214, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00063192793031220879, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060657401578200095, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00039873402383468892, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055351963726126224, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1639385612741475e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022365616007718661, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046704699702541431, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032230681977380224, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011712553572388214, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -0.00013179175250738632, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 7.1806458403162173e-05, doublethreshold); @@ -223,18 +330,42 @@ TEST_F(Verlet_test, rescaling) mdrun->second_half(); ; - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00063192793031220879, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060657401578200095, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00039873402383468892, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055351963726126224, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1639385612741475e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022365616007718661, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046704699702541431, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032230681977380224, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011712553572388214, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00063192793031220879, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060657401578200095, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00039873402383468892, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055351963726126224, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1639385612741475e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022365616007718661, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046704699702541431, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032230681977380224, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011712553572388214, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -0.00013178559069770653, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 7.1803101154153484e-05, doublethreshold); @@ -258,18 +389,42 @@ TEST_F(Verlet_test, rescale_v) mdrun->second_half(); ; - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00063192793031220879, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060657401578200095, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00039873402383468892, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055351963726126224, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1639385612741475e-05, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022365616007718661, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046704699702541431, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032230681977380224, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011712553572388214, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).x, 0.00063192793031220879, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).y, -0.00060657401578200095, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(1)]).z, 0.00039873402383468892, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).x, -0.00055351963726126224, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).y, -1.1639385612741475e-05, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(2)]).z, -0.00022365616007718661, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).x, 0.00046704699702541431, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).y, 0.00032230681977380224, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(3)]).z, -0.00011712553572388214, + doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -0.00013178559069770653, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 7.1803101154153484e-05, doublethreshold); @@ -296,9 +451,15 @@ TEST_F(Verlet_test, CSVR) mdrun->second_half(); // Check that positions are updated correctly - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, doublethreshold); - EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, -0.00054545529007222658, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).y, 0.00029590658162135359, + doublethreshold); + EXPECT_NEAR( + Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).z, -5.7952328034033513e-05, + doublethreshold); // Check that temperature is in reasonable range double temp = mdrun->t_current * ModuleBase::Hartree_to_K; diff --git a/source/source_md/unittests/test_md_func.cpp b/source/source_md/unittests/test_md_func.cpp index 986dec221a..c0da47e4d2 100644 --- a/source/source_md/unittests/test_md_func.cpp +++ b/source/source_md/unittests/test_md_func.cpp @@ -1,14 +1,13 @@ -#include "source_md/md_func.h" #include "source_base/parallel_cell.h" #include "source_base/parallel_common.h" #include "source_base/parallel_global.h" - -#include +#include "source_md/md_func.h" #include #include #include #include +#include #include #include #include @@ -17,13 +16,13 @@ namespace { -std::string make_dump_directory(const ModuleBase::CommunicationDomain& domain) +std::string make_dump_directory(const ModuleBase::CommunicationDomain &domain) { std::string directory; if (domain.rank() == 0) { std::string pattern = testing::TempDir() + "abacus-md-dump-XXXXXX"; - char* created = mkdtemp(&pattern[0]); + char *created = mkdtemp(&pattern[0]); if (created != nullptr) { directory = created; @@ -68,14 +67,38 @@ TEST_P(MdDumpBatchTest, OrdersAtomsAcrossRanksAndAppendsCompleteFrames) const std::vector masses(2, 1.0); const std::vector counts = {6, 7}; MDCell cell; - cell.initialize_from_owned_atoms(lattice, inverse, 1.0, 1.0, total_atoms, atoms, labels, masses, counts, 0.0, domain); + cell.initialize_from_owned_atoms( + lattice, + inverse, + 1.0, + 1.0, + total_atoms, + atoms, + labels, + masses, + counts, + 0.0, + domain); const std::string directory = make_dump_directory(domain); ASSERT_FALSE(directory.empty()); MD_para parameters; ModuleBase::matrix virial(3, 3); - MD_func::dump_info(0, directory, cell, parameters, false, virial); + MD_func::dump_info( + 0, + directory, + cell, + parameters, + false, + virial); const int batch_bytes = GetParam(); - MD_func::dump_info(1, directory, cell, parameters, false, virial, batch_bytes); + MD_func::dump_info( + 1, + directory, + cell, + parameters, + false, + virial, + batch_bytes); if (domain.rank() == 0) { const std::string filename = directory + "MD_dump"; @@ -165,7 +188,13 @@ TEST_P(MdDumpBatchTest, OrdersAtomsAcrossRanksAndAppendsCompleteFrames) } // One-row, several-row, and normal budgets must produce identical complete frames. -INSTANTIATE_TEST_SUITE_P(BufferBudgets, MdDumpBatchTest, testing::Values(1, 512, 1024, 1048576)); +INSTANTIATE_TEST_SUITE_P( + BufferBudgets, MdDumpBatchTest, + testing::Values( + 1, + 512, + 1024, + 1048576)); TEST(MdDumpMpiTest, BoundedGatherRejectsOversizedBatchesCollectively) { @@ -197,12 +226,29 @@ TEST(MdDumpMpiTest, EmptyRanksCanWriteAnEmptyFrame) const std::vector masses(1, 1.0); const std::vector counts(1, 0); MDCell cell; - cell.initialize_from_owned_atoms(lattice, inverse, 1.0, 1.0, 0, atoms, labels, masses, counts, 0.0, domain); + cell.initialize_from_owned_atoms( + lattice, + inverse, + 1.0, + 1.0, + 0, + atoms, + labels, + masses, + counts, + 0.0, + domain); const std::string directory = make_dump_directory(domain); ASSERT_FALSE(directory.empty()); MD_para parameters; ModuleBase::matrix virial(3, 3); - MD_func::dump_info(0, directory, cell, parameters, false, virial); + MD_func::dump_info( + 0, + directory, + cell, + parameters, + false, + virial); if (domain.rank() == 0) { const std::string filename = directory + "MD_dump"; @@ -218,12 +264,17 @@ TEST(MdDumpMpiTest, EmptyRanksCanWriteAnEmptyFrame) } } // namespace -int main(int argc, char** argv) +int main(int argc, char **argv) { int process_count = 1; int thread_count = 1; int rank = 0; - Parallel_Global::read_pal_param(argc, argv, process_count, thread_count, rank); + Parallel_Global::read_pal_param( + argc, + argv, + process_count, + thread_count, + rank); // These application communicators are not created by this focused test. POOL_WORLD = MPI_COMM_NULL; KP_WORLD = MPI_COMM_NULL; diff --git a/source/source_md/verlet.cpp b/source/source_md/verlet.cpp index 05c9c0f38a..b4db033128 100644 --- a/source/source_md/verlet.cpp +++ b/source/source_md/verlet.cpp @@ -7,21 +7,24 @@ #include #endif -Verlet::Verlet(const MD_para& mdp_in, - const bool cal_stress_in, - const bool init_vel, - const int my_rank_in, - MDCell& mdcell_in) - : MD_base(mdp_in, cal_stress_in, init_vel, my_rank_in, mdcell_in) +Verlet::Verlet( + const MD_para &mdp_in, + const bool cal_stress_in, + const bool init_vel, + const int my_rank_in, + MDCell &mdcell_in) + : MD_base( + mdp_in, + cal_stress_in, + init_vel, + my_rank_in, + mdcell_in) { } -Verlet::~Verlet() -{ -} +Verlet::~Verlet() {} - -void Verlet::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp) +void Verlet::setup(ModuleESolver::ESolver *p_esolver, const std::string &global_readin_dir, DomainDecomposition &decomp) { ModuleBase::TITLE("Verlet", "setup"); ModuleBase::timer::start("Verlet", "setup"); @@ -31,8 +34,7 @@ void Verlet::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_ ModuleBase::timer::end("Verlet", "setup"); } - -void Verlet::first_half(std::ofstream& ofs) +void Verlet::first_half(std::ofstream &ofs) { ModuleBase::TITLE("Verlet", "first_half"); ModuleBase::timer::start("Verlet", "first_half"); @@ -43,7 +45,6 @@ void Verlet::first_half(std::ofstream& ofs) ModuleBase::timer::end("Verlet", "first_half"); } - void Verlet::second_half() { ModuleBase::TITLE("Verlet", "second_half"); @@ -55,7 +56,6 @@ void Verlet::second_half() ModuleBase::timer::end("Verlet", "second_half"); } - void Verlet::apply_thermostat(void) { double t_target = 0.0; @@ -66,7 +66,11 @@ void Verlet::apply_thermostat(void) } else if (mdp.md_thermostat == "rescaling") { - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); + t_target = MD_func::target_temp( + step_ + step_rst_, + mdp.md_nstep, + md_tfirst, + md_tlast); if (std::abs(t_target - t_current) * ModuleBase::Hartree_to_K > mdp.md_tolerance) { thermalize(0, t_current, t_target); @@ -76,13 +80,17 @@ void Verlet::apply_thermostat(void) { if ((step_ + step_rst_) % mdp.md_nraise == 0) { - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); + t_target = MD_func::target_temp( + step_ + step_rst_, + mdp.md_nstep, + md_tfirst, + md_tlast); thermalize(0, t_current, t_target); } } else if (mdp.md_thermostat == "anderson") { - for (LocalAtom& atom : mdcell.owned_atoms()) + for (LocalAtom &atom : mdcell.owned_atoms()) { if (static_cast(std::rand()) / RAND_MAX <= 1.0 / mdp.md_nraise) { @@ -99,12 +107,20 @@ void Verlet::apply_thermostat(void) } else if (mdp.md_thermostat == "berendsen") { - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); + t_target = MD_func::target_temp( + step_ + step_rst_, + mdp.md_nstep, + md_tfirst, + md_tlast); thermalize(mdp.md_nraise, t_current, t_target); } else if (mdp.md_thermostat == "csvr") { - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); + t_target = MD_func::target_temp( + step_ + step_rst_, + mdp.md_nstep, + md_tfirst, + md_tlast); apply_csvr(t_current, t_target); } else @@ -113,8 +129,7 @@ void Verlet::apply_thermostat(void) } } - -void Verlet::thermalize(const int& nraise, const double& current_temp, const double& target_temp) +void Verlet::thermalize(const int &nraise, const double ¤t_temp, const double &target_temp) { double fac = 0.0; if (nraise > 0 && current_temp > 0 && target_temp > 0) @@ -126,11 +141,13 @@ void Verlet::thermalize(const int& nraise, const double& current_temp, const dou fac = sqrt(target_temp / current_temp); } - for (LocalAtom& atom : mdcell.owned_atoms()) atom.vel *= fac; + for (LocalAtom &atom : mdcell.owned_atoms()) + { + atom.vel *= fac; + } } - -void Verlet::apply_csvr(const double& current_temp, const double& target_temp) +void Verlet::apply_csvr(const double ¤t_temp, const double &target_temp) { // CSVR thermostat: Canonical Sampling through Velocity Rescaling // Reference: G. Bussi, D. Donadio, M. Parrinello, J. Chem. Phys. 126, 014101 (2007) @@ -144,8 +161,8 @@ void Verlet::apply_csvr(const double& current_temp, const double& target_temp) std::int64_t ndeg = MD_func::global_dof(mdcell); // Calculate kinetic energies - double kin_energy = current_temp * static_cast(ndeg) * 0.5; // in Hartree - double kin_target = target_temp * ndeg * 0.5; // in Hartree + double kin_energy = current_temp * static_cast(ndeg) * 0.5; // in Hartree + double kin_target = target_temp * ndeg * 0.5; // in Hartree // Calculate tau parameter (characteristic time scale / dt) double taut = mdp.md_csvr_tau / mdp.md_dt; @@ -171,12 +188,16 @@ void Verlet::apply_csvr(const double& current_temp, const double& target_temp) sumnoises += random_value * random_value; } const double factor2 = (1.0 - factor) * kin_target / kin_energy / ndeg; - const double resample = std::max(0.0, - factor + factor2 * (rr * rr + sumnoises) - + 2.0 * rr * sqrt(factor * factor2)); + const double resample = + std::max(0.0, factor + factor2 * (rr * rr + sumnoises) + 2.0 * rr * sqrt(factor * factor2)); scale = sqrt(resample); } - MPI_Bcast(&scale, 1, MPI_DOUBLE, 0, mdcell.communicator()); + MPI_Bcast( + &scale, + 1, + MPI_DOUBLE, + 0, + mdcell.communicator()); } else #endif @@ -189,32 +210,31 @@ void Verlet::apply_csvr(const double& current_temp, const double& target_temp) sumnoises += random_value * random_value; } const double factor2 = (1.0 - factor) * kin_target / kin_energy / ndeg; - const double resample = std::max(0.0, - factor + factor2 * (rr * rr + sumnoises) - + 2.0 * rr * sqrt(factor * factor2)); + const double resample = + std::max(0.0, factor + factor2 * (rr * rr + sumnoises) + 2.0 * rr * sqrt(factor * factor2)); scale = sqrt(resample); } // Apply velocity scaling - for (LocalAtom& atom : mdcell.owned_atoms()) atom.vel *= scale; + for (LocalAtom &atom : mdcell.owned_atoms()) + { + atom.vel *= scale; + } } - -void Verlet::print_md(std::ofstream& ofs, const bool& cal_stress) +void Verlet::print_md(std::ofstream &ofs, const bool &cal_stress) { MD_base::print_md(ofs, cal_stress); return; } - -void Verlet::write_restart(const std::string& global_out_dir) +void Verlet::write_restart(const std::string &global_out_dir) { MD_base::write_restart(global_out_dir); return; } - -void Verlet::restart(const std::string& global_readin_dir) +void Verlet::restart(const std::string &global_readin_dir) { MD_base::restart(global_readin_dir); return; diff --git a/source/source_md/verlet.h b/source/source_md/verlet.h index 670688251e..6bd1563ebe 100644 --- a/source/source_md/verlet.h +++ b/source/source_md/verlet.h @@ -10,20 +10,21 @@ class Verlet : public MD_base { public: - Verlet(const MD_para& mdp_in, - const bool cal_stress_in, - const bool init_vel, - const int my_rank_in, - MDCell& mdcell_in); + Verlet( + const MD_para &mdp_in, + const bool cal_stress_in, + const bool init_vel, + const int my_rank_in, + MDCell &mdcell_in); ~Verlet(); private: - void setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp); - void first_half(std::ofstream& ofs); + void setup(ModuleESolver::ESolver *p_esolver, const std::string &global_readin_dir, DomainDecomposition &decomp); + void first_half(std::ofstream &ofs); void second_half(); - void restart(const std::string& global_readin_dir); - void print_md(std::ofstream& ofs, const bool& cal_stress); - void write_restart(const std::string& global_out_dir); + void restart(const std::string &global_readin_dir); + void print_md(std::ofstream &ofs, const bool &cal_stress); + void write_restart(const std::string &global_out_dir); /** * @brief apply specifical thermostats according to the input para @@ -38,7 +39,7 @@ class Verlet : public MD_base * @param current_temp the current temperature * @param target_temp the target temperature */ - void thermalize(const int& nraise, const double& current_temp, const double& target_temp); + void thermalize(const int &nraise, const double ¤t_temp, const double &target_temp); /** * @brief apply CSVR thermostat @@ -46,7 +47,7 @@ class Verlet : public MD_base * @param current_temp the current temperature * @param target_temp the target temperature */ - void apply_csvr(const double& current_temp, const double& target_temp); + void apply_csvr(const double ¤t_temp, const double &target_temp); }; #endif