From 90a2f71ebe43a1a895de3aa01d5ffbf9d410d663 Mon Sep 17 00:00:00 2001 From: SY Wang Date: Tue, 20 Oct 2026 20:49:44 +0800 Subject: [PATCH 1/3] Precommit: Prettify source_md --- source/source_md/fire.cpp | 24 +- source/source_md/fire.h | 1 - source/source_md/langevin.cpp | 24 +- source/source_md/langevin.h | 4 +- source/source_md/md_base.cpp | 33 +- source/source_md/md_base.h | 28 +- source/source_md/md_func.cpp | 221 ++++++++----- source/source_md/md_func.h | 4 +- .../module_plumed/plumed_interface.cpp | 12 +- .../module_plumed/plumed_interface.h | 9 +- source/source_md/msst.cpp | 20 +- source/source_md/msst.h | 31 +- source/source_md/nhchain.cpp | 17 +- source/source_md/nhchain.h | 52 +-- source/source_md/run_md.cpp | 40 +-- source/source_md/run_md.h | 5 +- source/source_md/test/fire_test.cpp | 104 ++++-- source/source_md/test/langevin_test.cpp | 104 ++++-- source/source_md/test/lj_pot_test.cpp | 9 +- source/source_md/test/md_func_test.cpp | 54 +-- source/source_md/test/md_test_fixture.h | 6 +- source/source_md/test/msst_test.cpp | 104 ++++-- source/source_md/test/nhchain_test.cpp | 104 ++++-- source/source_md/test/run_md_test.cpp | 7 +- source/source_md/test/setcell.h | 34 +- source/source_md/test/test_md_func.cpp | 37 ++- source/source_md/test/verlet_test.cpp | 309 +++++++++++++----- source/source_md/verlet.cpp | 31 +- 28 files changed, 901 insertions(+), 527 deletions(-) diff --git a/source/source_md/fire.cpp b/source/source_md/fire.cpp index cb4611c6e9c..bd037eaf757 100644 --- a/source/source_md/fire.cpp +++ b/source/source_md/fire.cpp @@ -66,7 +66,6 @@ void FIRE::first_half(std::ofstream& ofs) return; } - void FIRE::second_half(void) { ModuleBase::TITLE("FIRE", "second_half"); @@ -81,7 +80,6 @@ void FIRE::second_half(void) return; } - void FIRE::print_md(std::ofstream& ofs, const bool& cal_stress) { MD_base::print_md(ofs, cal_stress); @@ -94,7 +92,6 @@ void FIRE::print_md(std::ofstream& ofs, const bool& cal_stress) return; } - void FIRE::write_restart(const std::string& global_out_dir) { if (!my_rank) @@ -118,7 +115,6 @@ void FIRE::write_restart(const std::string& global_out_dir) return; } - void FIRE::restart(const std::string& global_readin_dir) { bool ok = true; @@ -137,10 +133,10 @@ 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(); } } @@ -172,7 +168,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,7 +193,7 @@ void FIRE::check_force(void) } } - #ifdef __MPI +#ifdef __MPI if (mdcell.mpi_size() > 1) { double global_max = 0.0; @@ -224,7 +220,6 @@ void FIRE::check_force(void) return; } - void FIRE::check_fire(void) { double P = 0.0; @@ -242,7 +237,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,7 +256,7 @@ void FIRE::check_fire(void) } } - #ifdef __MPI +#ifdef __MPI if (mdcell.mpi_size() > 1) { double local_values[3] = {P, sumforce, normvel}; @@ -289,7 +284,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 +294,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 +313,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 5b84c3be9ce..4649ef85c1c 100644 --- a/source/source_md/fire.h +++ b/source/source_md/fire.h @@ -42,7 +42,6 @@ class FIRE : public MD_base } private: - void setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp); void first_half(std::ofstream& ofs); diff --git a/source/source_md/langevin.cpp b/source/source_md/langevin.cpp index 204e946bbc6..f09df148594 100644 --- a/source/source_md/langevin.cpp +++ b/source/source_md/langevin.cpp @@ -1,8 +1,8 @@ #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, @@ -12,15 +12,16 @@ Langevin::Langevin(const MD_para& mdp_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,7 +34,6 @@ void Langevin::setup(ModuleESolver::ESolver* p_esolver, const std::string& globa return; } - void Langevin::first_half(std::ofstream& ofs) { ModuleBase::TITLE("Langevin", "first_half"); @@ -44,7 +44,8 @@ void Langevin::first_half(std::ofstream& ofs) 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 +54,6 @@ void Langevin::first_half(std::ofstream& ofs) return; } - void Langevin::second_half() { ModuleBase::TITLE("Langevin", "second_half"); @@ -65,7 +65,8 @@ void Langevin::second_half() 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,28 +74,24 @@ void Langevin::second_half() return; } - 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) { MD_base::write_restart(global_out_dir); return; } - 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); @@ -111,8 +108,7 @@ void Langevin::post_force() 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 58ad4b63a66..48e1185c270 100644 --- a/source/source_md/langevin.h +++ b/source/source_md/langevin.h @@ -41,8 +41,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 cb00647a865..fa7ac21f4a6 100644 --- a/source/source_md/md_base.cpp +++ b/source/source_md/md_base.cpp @@ -1,10 +1,12 @@ #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 @@ -13,7 +15,7 @@ MD_base::MD_base(const MD_para& mdp_in, const bool init_vel, const int my_rank_in, MDCell& mdcell_in) -: mdp(mdp_in), mdcell(mdcell_in) + : mdp(mdp_in), mdcell(mdcell_in) { #ifdef __MPI my_rank = mdcell.mpi_rank(); @@ -28,8 +30,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; @@ -43,11 +45,13 @@ MD_base::MD_base(const MD_para& mdp_in, t_current = MD_func::current_temp(kinetic, mdcell, frozen_freedom_); } +MD_base::~MD_base() +{ +} -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,7 +63,7 @@ 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); + 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); @@ -71,7 +75,6 @@ void MD_base::setup(ModuleESolver::ESolver* p_esolver, const std::string& global return; } - void MD_base::first_half(std::ofstream& ofs) { update_vel(); @@ -80,7 +83,6 @@ void MD_base::first_half(std::ofstream& ofs) return; } - void MD_base::second_half() { update_vel(); @@ -88,7 +90,6 @@ void MD_base::second_half() return; } - void MD_base::update_pos() { std::vector& atoms = mdcell.owned_atoms(); @@ -118,7 +119,6 @@ void MD_base::update_pos() return; } - void MD_base::update_vel() { std::vector& atoms = mdcell.owned_atoms(); @@ -136,17 +136,16 @@ void MD_base::update_vel() return; } - 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; @@ -190,7 +189,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 << " ------------------------------------------------------------------------------------------------" @@ -222,7 +221,6 @@ void MD_base::print_md(std::ofstream& ofs, const bool& cal_stress) return; } - void MD_base::write_restart(const std::string& global_out_dir) { if (!my_rank) @@ -242,7 +240,6 @@ void MD_base::write_restart(const std::string& global_out_dir) return; } - void MD_base::restart(const std::string& global_readin_dir) { MD_func::current_md_info(mdcell, global_readin_dir, step_rst_, md_tfirst); diff --git a/source/source_md/md_base.h b/source/source_md/md_base.h index 4095c409a0c..900392b8d96 100644 --- a/source/source_md/md_base.h +++ b/source/source_md/md_base.h @@ -43,7 +43,9 @@ 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 @@ -100,20 +102,20 @@ class MD_base 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 + 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 a549e87732c..23e25e0e779 100644 --- a/source/source_md/md_func.cpp +++ b/source/source_md/md_func.cpp @@ -1,8 +1,8 @@ #include "md_func.h" -#include "source_cell/domain_decomposition.h" #include "source_base/global_variable.h" #include "source_base/timer.h" +#include "source_cell/domain_decomposition.h" #include "source_io/module_output/output_log.h" #include @@ -10,7 +10,6 @@ #include #include - namespace MD_func { #ifdef __MPI @@ -21,23 +20,24 @@ bool write_dump_at(const int file, const std::string& data, const MPI_Offset off std::size_t written = 0; while (written < data.size()) { - const ssize_t count = pwrite(file, data.data() + written, data.size() - written, - static_cast(offset + written)); - if (count <= 0) return false; + const ssize_t count + = pwrite(file, data.data() + written, data.size() - written, static_cast(offset + written)); + if (count <= 0) + return false; written += static_cast(count); } return true; } -} +} // namespace #endif 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) { @@ -67,7 +67,7 @@ double kinetic_energy(const int& natom, const ModuleBase::Vector3* vel, { 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(); @@ -81,19 +81,22 @@ void compute_stress(const MDCell& mdcell, 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]; + 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()); #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) @@ -214,14 +217,17 @@ void init_vel(MDCell& mdcell, { 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()); 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) { @@ -259,9 +265,10 @@ 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; } @@ -269,7 +276,9 @@ void init_vel(MDCell& mdcell, 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) { @@ -277,9 +286,12 @@ void init_vel(MDCell& mdcell, } #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; @@ -287,20 +299,22 @@ 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 @@ -311,7 +325,8 @@ void init_vel(MDCell& mdcell, { 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; @@ -319,7 +334,8 @@ 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; } } @@ -340,18 +356,22 @@ void force_virial(ModuleESolver::ESolver* p_esolver, 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)); + 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; @@ -363,16 +383,28 @@ void force_virial(ModuleESolver::ESolver* p_esolver, 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; itrunner(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)); + } } ModuleBase::timer::end("MD_func", "force_virial"); } - void print_stress(std::ofstream& ofs, const ModuleBase::matrix& virial, const ModuleBase::matrix& stress) { double stress_scalar = 0.0; @@ -386,23 +418,23 @@ void print_stress(std::ofstream& ofs, const ModuleBase::matrix& virial, const Mo const double unit_transform = ModuleBase::HARTREE_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; - ofs << " ELECTRONIC PART OF STRESS: " << virial_scalar * unit_transform << " kbar" << std::endl; - ofs << " IONIC (KINETIC) PART OF STRESS: " << (stress_scalar - virial_scalar) * unit_transform << " kbar" << std::endl; + ofs << " IONIC (KINETIC) PART OF STRESS: " << (stress_scalar - virial_scalar) * unit_transform << " kbar" + << std::endl; ofs << " MD PRESSURE (ELECTRONS+IONS) : " << stress_scalar * unit_transform << " kbar" << std::endl; // one should use 'print_stress' function in ../source/source_io/output_log.cpp -/* - ofs.unsetf(std::ios::fixed); - ofs << std::setprecision(8) << std::endl; - ModuleBase::GlobalFunc::NEW_PART("MD STRESS (kbar)"); - for (int i = 0; i < 3; i++) - { - ofs << std::setw(15) << stress(i, 0) * unit_transform << std::setw(15) << stress(i, 1) * unit_transform - << std::setw(15) << stress(i, 2) * unit_transform << std::endl; - } - ofs << std::setiosflags(std::ios::left); -*/ + /* + ofs.unsetf(std::ios::fixed); + ofs << std::setprecision(8) << std::endl; + ModuleBase::GlobalFunc::NEW_PART("MD STRESS (kbar)"); + for (int i = 0; i < 3; i++) + { + ofs << std::setw(15) << stress(i, 0) * unit_transform << std::setw(15) << stress(i, 1) * unit_transform + << std::setw(15) << stress(i, 2) * unit_transform << std::endl; + } + ofs << std::setiosflags(std::ios::left); + */ return; } @@ -430,11 +462,14 @@ 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) @@ -446,11 +481,12 @@ void dump_info(const int& step, for (int i = 0; i < mdcell.owned_atoms().size(); ++i) { const LocalAtom& atom = mdcell.owned_atoms()[static_cast(i)]; - local << " " << 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; + local << " " << 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) - local << " " << atom.force.x * unit_force << " " << atom.force.y * unit_force << " " << atom.force.z * unit_force; + local << " " << atom.force.x * unit_force << " " << atom.force.y * unit_force << " " + << atom.force.z * unit_force; if (mdp.dump_vel) local << " " << atom.vel.x * unit_vel << " " << atom.vel.y * unit_vel << " " << atom.vel.z * unit_vel; local << "\n"; @@ -458,23 +494,30 @@ void dump_info(const int& step, local << "\n\n"; #ifdef __MPI const MPI_Comm comm = mdcell.communicator(); - int rank = 0; MPI_Comm_rank(comm, &rank); + int rank = 0; + MPI_Comm_rank(comm, &rank); MPI_Offset base = 0; if (rank == 0) { - const int file_descriptor = open(file.str().c_str(), O_CREAT | O_WRONLY | (step == 0 ? O_TRUNC : O_APPEND), 0666); - if (file_descriptor < 0) ModuleBase::WARNING_QUIT("MD_func::dump_info", "cannot open MD_dump."); + const int file_descriptor + = open(file.str().c_str(), O_CREAT | O_WRONLY | (step == 0 ? O_TRUNC : O_APPEND), 0666); + if (file_descriptor < 0) + ModuleBase::WARNING_QUIT("MD_func::dump_info", "cannot open MD_dump."); base = lseek(file_descriptor, 0, SEEK_END); - if (!write_dump_at(file_descriptor, header.str(), base) || close(file_descriptor) != 0) ModuleBase::WARNING_QUIT("MD_func::dump_info", "cannot write MD_dump header."); + if (!write_dump_at(file_descriptor, header.str(), base) || close(file_descriptor) != 0) + ModuleBase::WARNING_QUIT("MD_func::dump_info", "cannot write MD_dump header."); } MPI_Bcast(&base, 1, MPI_OFFSET, 0, comm); const MPI_Offset atom_offset = base + static_cast(header.str().size()); const MPI_Offset local_size = static_cast(local.str().size()); MPI_Offset rank_offset = 0; MPI_Exscan(&local_size, &rank_offset, 1, MPI_OFFSET, MPI_SUM, comm); - if (rank == 0) rank_offset = 0; + if (rank == 0) + rank_offset = 0; const int file_descriptor = open(file.str().c_str(), O_WRONLY); - if (file_descriptor < 0 || !write_dump_at(file_descriptor, local.str(), atom_offset + rank_offset) || close(file_descriptor) != 0) ModuleBase::WARNING_QUIT("MD_func::dump_info", "cannot write MD_dump atoms."); + if (file_descriptor < 0 || !write_dump_at(file_descriptor, local.str(), atom_offset + rank_offset) + || close(file_descriptor) != 0) + ModuleBase::WARNING_QUIT("MD_func::dump_info", "cannot write MD_dump atoms."); #else std::ofstream ofs(file.str().c_str(), step == 0 ? std::ios::trunc : std::ios::app); ofs << header.str() << local.str(); @@ -494,15 +537,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++; } @@ -513,7 +559,7 @@ void get_mass_mbl(const UnitCell& unit_in, 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); } @@ -533,9 +579,7 @@ 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) @@ -547,7 +591,8 @@ double current_temp(double& kinetic, 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); } @@ -557,18 +602,24 @@ std::int64_t global_dof(const MDCell& mdcell) 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]; + 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()); #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; } diff --git a/source/source_md/md_func.h b/source/source_md/md_func.h index c8095d069e8..f31d1b93fdc 100644 --- a/source/source_md/md_func.h +++ b/source/source_md/md_func.h @@ -188,9 +188,7 @@ double current_temp(double& kinetic, 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); +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); diff --git a/source/source_md/module_plumed/plumed_interface.cpp b/source/source_md/module_plumed/plumed_interface.cpp index ca92f8f3d81..4588e1ce9a2 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 @@ -36,10 +36,7 @@ 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; @@ -136,8 +133,7 @@ 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); diff --git a/source/source_md/module_plumed/plumed_interface.h b/source/source_md/module_plumed/plumed_interface.h index 8a1d08e4a3c..5749781fe2c 100644 --- a/source/source_md/module_plumed/plumed_interface.h +++ b/source/source_md/module_plumed/plumed_interface.h @@ -1,11 +1,11 @@ #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 @@ -72,10 +72,7 @@ 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 diff --git a/source/source_md/msst.cpp b/source/source_md/msst.cpp index 02f311a7395..03fb73aa1be 100644 --- a/source/source_md/msst.cpp +++ b/source/source_md/msst.cpp @@ -29,7 +29,8 @@ 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()); #endif @@ -67,7 +68,8 @@ 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); @@ -162,7 +164,6 @@ void MSST::second_half() return; } - void MSST::print_md(std::ofstream& ofs, const bool& cal_stress) { MD_base::print_md(ofs, cal_stress); @@ -195,7 +196,6 @@ void MSST::write_restart(const std::string& global_out_dir) return; } - void MSST::restart(const std::string& global_readin_dir) { bool ok = true; @@ -214,7 +214,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; } @@ -248,7 +248,8 @@ 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()); #endif @@ -270,18 +271,18 @@ 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; @@ -306,7 +307,6 @@ void MSST::propagate_vel() return; } - void MSST::propagate_voldot() { const int sd = mdp.msst_direction; diff --git a/source/source_md/msst.h b/source/source_md/msst.h index 00ba2a7914a..f230522141c 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,11 +16,7 @@ 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 @@ -87,18 +84,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 9a99f1b05e8..4c4bd3a8502 100644 --- a/source/source_md/nhchain.cpp +++ b/source/source_md/nhchain.cpp @@ -17,7 +17,7 @@ Nose_Hoover::Nose_Hoover(const MD_para& mdp_in, 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; @@ -95,7 +95,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 +159,9 @@ 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"); @@ -273,7 +275,6 @@ void Nose_Hoover::first_half(std::ofstream& ofs) return; } - void Nose_Hoover::second_half() { ModuleBase::TITLE("Nose_Hoover", "second_half"); @@ -317,7 +318,6 @@ 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); @@ -555,7 +555,8 @@ 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,7 +650,7 @@ 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]; @@ -700,7 +701,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) { diff --git a/source/source_md/nhchain.h b/source/source_md/nhchain.h index 6e0db3bec60..ee8330e7f2b 100644 --- a/source/source_md/nhchain.h +++ b/source/source_md/nhchain.h @@ -109,33 +109,33 @@ class Nose_Hoover : public MD_base 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 114f9cbb3a5..c176aa8212e 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 @@ -129,10 +129,7 @@ void collect_plumed_cell(const MDCell& mdcell, double* cell) } } // namespace -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& param_in, DomainDecomposition& decomp) { ModuleBase::TITLE("Run_MD", "md_line"); ModuleBase::timer::start("Run_MD", "md_line"); @@ -259,18 +256,12 @@ 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, @@ -280,8 +271,7 @@ void md_line(MDCell& mdcell, 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 8abf531a571..2beaf7c93a4 100644 --- a/source/source_md/run_md.h +++ b/source/source_md/run_md.h @@ -28,10 +28,7 @@ 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& param_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 8655b5aa90b..05ae519c433 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 /************************************************ @@ -82,18 +84,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 +140,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); diff --git a/source/source_md/test/langevin_test.cpp b/source/source_md/test/langevin_test.cpp index 784b6a9fd18..cbb62d1a355 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 /************************************************ @@ -82,18 +84,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 +141,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 0ae49fb2b41..05d619658cf 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 /************************************************ @@ -205,7 +206,9 @@ 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)); diff --git a/source/source_md/test/md_func_test.cpp b/source/source_md/test/md_func_test.cpp index 7eeccd78691..0ae45e9b06b 100644 --- a/source/source_md/test/md_func_test.cpp +++ b/source/source_md/test/md_func_test.cpp @@ -1,8 +1,10 @@ +#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 @@ -136,7 +138,7 @@ 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; } @@ -157,7 +159,7 @@ 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); } @@ -309,27 +311,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"); @@ -360,7 +362,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; } diff --git a/source/source_md/test/md_test_fixture.h b/source/source_md/test/md_test_fixture.h index 853ccd1a3fc..dac429f978a 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 diff --git a/source/source_md/test/msst_test.cpp b/source/source_md/test/msst_test.cpp index 6a7e7cd3d29..b037304de20 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 /************************************************ @@ -107,18 +109,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 +179,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); diff --git a/source/source_md/test/nhchain_test.cpp b/source/source_md/test/nhchain_test.cpp index a73b9811c41..e0756fe9675 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 @@ -84,18 +86,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 +142,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); diff --git a/source/source_md/test/run_md_test.cpp b/source/source_md/test/run_md_test.cpp index 9c459d3e886..52b4e4f14f0 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) { diff --git a/source/source_md/test/setcell.h b/source/source_md/test/setcell.h index 0cd0809c1ab..b6027cdb293 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 @@ -35,7 +35,7 @@ class Setcell ucell.atoms = new Atom[ucell.ntype]; ucell.set_atom_flag = true; - + ucell.atoms[0].mass = 39.948; ucell.atoms[0].label = "Ar"; @@ -154,14 +154,14 @@ class Setcell 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)); + 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)); } } @@ -182,12 +182,10 @@ class Setcell static ModuleBase::Vector3 fractional_displacement(const LocalAtom& atom) { - 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) - }; + 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 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_func.cpp b/source/source_md/test/test_md_func.cpp index 9e0aef2f6ee..c4d36889641 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,25 +10,32 @@ namespace { class GhostForceSolver : public ModuleESolver::ESolver { -public: - void before_all_runners(BaseCell&, const Input_para&) override {} - void after_all_runners(BaseCell&) 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); 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; } + double cal_energy() override + { + return 6.0; + } void cal_force(BaseCell&, ModuleBase::matrix&) override { ADD_FAILURE() << "Direct MDCell forces must not use cal_force."; @@ -52,11 +59,17 @@ 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, + 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); + 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; @@ -69,7 +82,7 @@ TEST(MdForceVirialTest, ReturnsGhostForcesBeforeConvertingUnitsExactlyOnce) 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 8d960f04d38..bb285efb4c6 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 @@ -84,18 +85,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 +143,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 +202,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 +261,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 +320,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 +379,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 +441,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/verlet.cpp b/source/source_md/verlet.cpp index 05c9c0f38ae..bf4bc11a768 100644 --- a/source/source_md/verlet.cpp +++ b/source/source_md/verlet.cpp @@ -20,7 +20,6 @@ Verlet::~Verlet() { } - void Verlet::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, DomainDecomposition& decomp) { ModuleBase::TITLE("Verlet", "setup"); @@ -31,7 +30,6 @@ void Verlet::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_ ModuleBase::timer::end("Verlet", "setup"); } - void Verlet::first_half(std::ofstream& ofs) { ModuleBase::TITLE("Verlet", "first_half"); @@ -43,7 +41,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 +52,6 @@ void Verlet::second_half() ModuleBase::timer::end("Verlet", "second_half"); } - void Verlet::apply_thermostat(void) { double t_target = 0.0; @@ -82,7 +78,7 @@ void Verlet::apply_thermostat(void) } 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) { @@ -113,7 +109,6 @@ void Verlet::apply_thermostat(void) } } - void Verlet::thermalize(const int& nraise, const double& current_temp, const double& target_temp) { double fac = 0.0; @@ -126,10 +121,10 @@ 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) { // CSVR thermostat: Canonical Sampling through Velocity Rescaling @@ -144,8 +139,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,9 +166,8 @@ 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()); @@ -189,31 +183,28 @@ 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) { MD_base::print_md(ofs, cal_stress); return; } - 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) { MD_base::restart(global_readin_dir); From cca2cb3eacd38e7e99437af951acc3209b4d7c1f Mon Sep 17 00:00:00 2001 From: SY Wang Date: Tue, 27 Oct 2026 23:07:27 +0800 Subject: [PATCH 2/3] Try LLVM style, minor adjustment (modify .clang-format) --- .clang-format | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/.clang-format b/.clang-format index 8970df99c29..ced5885879f 100644 --- a/.clang-format +++ b/.clang-format @@ -1,6 +1,6 @@ --- Language: Cpp -BasedOnStyle: Microsoft +BasedOnStyle: LLVM AlwaysBreakTemplateDeclarations: Yes @@ -9,7 +9,7 @@ AllowAllParametersOfDeclarationOnNextLine: false BinPackArguments: false BinPackParameters: false -BreakBeforeBinaryOperators: All +BreakBeforeBinaryOperators: NonAssignment BreakBeforeTernaryOperators: true Cpp11BracedListStyle: true @@ -25,6 +25,11 @@ IncludeIsMainSourceRegex: '' IndentWidth: 4 IndentWrappedFunctionNames: true +ColumnLimit: 120 + +# Avoid awkward line breaks before member access operators (e.g., obj.method()). +PenaltyBreakBeforeMemberAccess: 1000 + # About Point DerivePointerAlignment: false PointerAlignment: Left From cf8459fa22bd5cb9e76ffbb771486a2e6804269c Mon Sep 17 00:00:00 2001 From: SY Wang Date: Tue, 27 Oct 2026 23:08:37 +0800 Subject: [PATCH 3/3] Try LLVM style, minor adjustment (formatting) --- source/source_md/fire.cpp | 116 ++---- source/source_md/fire.h | 18 +- source/source_md/langevin.cpp | 45 +- source/source_md/langevin.h | 3 +- source/source_md/md_base.cpp | 83 ++-- source/source_md/md_base.h | 13 +- source/source_md/md_func.cpp | 264 ++++-------- source/source_md/md_func.h | 3 +- source/source_md/md_statistics.h | 9 +- .../module_plumed/plumed_interface.cpp | 62 +-- .../module_plumed/plumed_interface.h | 14 +- source/source_md/msst.cpp | 99 ++--- source/source_md/msst.h | 28 +- source/source_md/nhchain.cpp | 389 ++++++------------ source/source_md/nhchain.h | 28 +- source/source_md/run_md.cpp | 106 ++--- source/source_md/run_md.h | 6 +- source/source_md/test/fire_test.cpp | 27 +- source/source_md/test/langevin_test.cpp | 27 +- source/source_md/test/lj_pot_test.cpp | 70 ++-- source/source_md/test/md_func_test.cpp | 52 +-- source/source_md/test/md_test_fixture.h | 24 +- source/source_md/test/msst_test.cpp | 27 +- source/source_md/test/nhchain_test.cpp | 27 +- source/source_md/test/run_md_test.cpp | 3 +- source/source_md/test/setcell.h | 31 +- source/source_md/test/test_md_func.cpp | 43 +- source/source_md/test/verlet_test.cpp | 42 +- source/source_md/verlet.cpp | 118 ++---- source/source_md/verlet.h | 3 +- 30 files changed, 578 insertions(+), 1202 deletions(-) diff --git a/source/source_md/fire.cpp b/source/source_md/fire.cpp index bd037eaf757..f95f51a7ff4 100644 --- a/source/source_md/fire.cpp +++ b/source/source_md/fire.cpp @@ -12,8 +12,7 @@ FIRE::FIRE(const MD_para& mdp_in, 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) -{ + : MD_base(mdp_in, cal_stress_in, init_vel, my_rank_in, mdcell_in) { force_thr = force_thr_in; dt_max = -1.0; alpha_start = 0.10; @@ -32,12 +31,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,8 +46,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,8 +61,7 @@ void FIRE::first_half(std::ofstream& ofs) return; } -void FIRE::second_half(void) -{ +void FIRE::second_half(void) { ModuleBase::TITLE("FIRE", "second_half"); ModuleBase::timer::start("FIRE", "second_half"); @@ -80,8 +74,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); const double max_force = max * ModuleBase::Hartree_to_eV * ModuleBase::ANGSTROM_AU; @@ -92,10 +85,8 @@ void FIRE::print_md(std::ofstream& ofs, const bool& cal_stress) return; } -void FIRE::write_restart(const std::string& global_out_dir) -{ - if (!my_rank) - { +void FIRE::write_restart(const std::string& global_out_dir) { + if (!my_rank) { std::stringstream ssc; ssc << global_out_dir << "Restart_md.txt"; std::ofstream file(ssc.str().c_str()); @@ -115,26 +106,21 @@ 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; - if (!my_rank) - { + if (!my_rank) { std::stringstream ssc; ssc << global_readin_dir << "Restart_md.txt"; std::ifstream file(ssc.str().c_str()); - if (!file) - { + if (!file) { ok = false; } - if (ok) - { + if (ok) { file >> step_rst_ >> md_tfirst >> alpha >> negative_count >> dt_max >> md_dt; - if (!file) - { + if (!file) { ok = false; } file.close(); @@ -145,8 +131,7 @@ void FIRE::restart(const std::string& global_readin_dir) MPI_Bcast(&ok, 1, MPI_C_BOOL, 0, mdcell.communicator()); #endif - if (!ok) - { + if (!ok) { ModuleBase::WARNING_QUIT("mdrun", "no Restart_md.txt !"); } @@ -162,16 +147,13 @@ void FIRE::restart(const std::string& global_readin_dir) return; } -void FIRE::check_force(void) -{ +void FIRE::check_force(void) { max = 0.0; std::int64_t movable_dof = 0; - for (const LocalAtom& atom: mdcell.owned_atoms()) - { - for (int j = 0; j < 3; ++j) - { + for (const LocalAtom& atom: mdcell.owned_atoms()) { + for (int j = 0; j < 3; ++j) { // Only movable degrees of freedom should be used // in the FIRE convergence criterion. // @@ -179,23 +161,20 @@ void FIRE::check_force(void) // m 1 1 1 -> x/y/z are included. // m 1 0 1 -> y is excluded. // m 0 0 0 -> this atom contributes no DOF to convergence. - if (!atom.mbl[j]) - { + if (!atom.mbl[j]) { continue; } ++movable_dof; - if (max < std::abs(atom.force[j])) - { + if (max < std::abs(atom.force[j])) { max = std::abs(atom.force[j]); } } } #ifdef __MPI - if (mdcell.mpi_size() > 1) - { + 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()); @@ -206,29 +185,25 @@ void FIRE::check_force(void) #endif // If there are no movable degrees of freedom, there is nothing to optimize. - if (movable_dof == 0) - { + if (movable_dof == 0) { stop = true; return; } - if (2.0 * max < force_thr) - { + if (2.0 * max < force_thr) { stop = true; } return; } -void FIRE::check_fire(void) -{ +void FIRE::check_fire(void) { double P = 0.0; double sumforce = 0.0; double normvel = 0.0; /// initial dt_max - if (dt_max < 0) - { + if (dt_max < 0) { dt_max = 2.5 * md_dt; } @@ -237,12 +212,9 @@ 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 (int j = 0; j < 3; ++j) - { - if (!atom.mbl[j]) - { + for (LocalAtom& atom: mdcell.owned_atoms()) { + for (int j = 0; j < 3; ++j) { + if (!atom.mbl[j]) { // Keep frozen components clean. atom.vel[j] = 0.0; continue; @@ -257,8 +229,7 @@ void FIRE::check_fire(void) } #ifdef __MPI - if (mdcell.mpi_size() > 1) - { + 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; @@ -272,8 +243,7 @@ void FIRE::check_fire(void) #endif // No movable degrees of freedom: nothing to update. - if (movable_dof == 0) - { + if (movable_dof == 0) { return; } @@ -282,14 +252,10 @@ void FIRE::check_fire(void) // If force or velocity norm is zero, the velocity projection is undefined. // 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 (int j = 0; j < 3; ++j) - { - if (!atom.mbl[j]) - { + if (sumforce > 0.0 && normvel > 0.0) { + for (LocalAtom& atom: mdcell.owned_atoms()) { + for (int j = 0; j < 3; ++j) { + if (!atom.mbl[j]) { atom.vel[j] = 0.0; continue; } @@ -299,24 +265,18 @@ void FIRE::check_fire(void) } } - if (P > 0) - { + if (P > 0) { negative_count++; - if (negative_count >= n_min) - { + if (negative_count >= n_min) { md_dt = std::min(md_dt * finc, dt_max); alpha *= f_alpha; } - } - else - { + } else { md_dt *= fdec; negative_count = 0; - for (LocalAtom& atom: mdcell.owned_atoms()) - { - for (int j = 0; j < 3; ++j) - { + for (LocalAtom& atom: mdcell.owned_atoms()) { + for (int j = 0; j < 3; ++j) { atom.vel[j] = 0; } } diff --git a/source/source_md/fire.h b/source/source_md/fire.h index 4649ef85c1c..ac78316a2f3 100644 --- a/source/source_md/fire.h +++ b/source/source_md/fire.h @@ -10,8 +10,7 @@ * It is based on conventional molecular dynamics with additional velocity modifications and adaptive time steps. * The MD trajectory will descend to an energy-minimum. */ -class FIRE : public MD_base -{ +class FIRE : public MD_base { public: /// @param force_thr_in force convergence threshold used by the FIRE descent FIRE(const MD_para& mdp_in, @@ -24,22 +23,13 @@ class FIRE : public MD_base ~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); diff --git a/source/source_md/langevin.cpp b/source/source_md/langevin.cpp index f09df148594..55d47a6237e 100644 --- a/source/source_md/langevin.cpp +++ b/source/source_md/langevin.cpp @@ -9,8 +9,7 @@ Langevin::Langevin(const MD_para& mdp_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) -{ + : 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); @@ -21,8 +20,7 @@ Langevin::Langevin(const MD_para& mdp_in, void Langevin::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, - DomainDecomposition& decomp) -{ + DomainDecomposition& decomp) { ModuleBase::TITLE("Langevin", "setup"); ModuleBase::timer::start("Langevin", "setup"); @@ -34,16 +32,13 @@ void Langevin::setup(ModuleESolver::ESolver* p_esolver, 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) - { + for (int i = 0; i < mdcell.owned_atoms().size(); ++i) { LocalAtom& atom = mdcell.owned_atoms()[static_cast(i)]; - for (int k = 0; k < 3; ++k) - { + for (int k = 0; k < 3; ++k) { if (atom.mbl[k]) atom.vel[k] += 0.5 * total_force[i][k] * md_dt / atom.mass; } @@ -54,17 +49,14 @@ void Langevin::first_half(std::ofstream& ofs) return; } -void Langevin::second_half() -{ +void Langevin::second_half() { ModuleBase::TITLE("Langevin", "second_half"); ModuleBase::timer::start("Langevin", "second_half"); post_force(); - for (int i = 0; i < mdcell.owned_atoms().size(); ++i) - { + for (int i = 0; i < mdcell.owned_atoms().size(); ++i) { LocalAtom& atom = mdcell.owned_atoms()[static_cast(i)]; - for (int k = 0; k < 3; ++k) - { + for (int k = 0; k < 3; ++k) { if (atom.mbl[k]) atom.vel[k] += 0.5 * total_force[i][k] * md_dt / atom.mass; } @@ -74,40 +66,33 @@ 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() -{ +void Langevin::post_force() { 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) - { + for (int i = 0; i < mdcell.owned_atoms().size(); ++i) { ModuleBase::Vector3 random_value; - for (int k = 0; k < 3; ++k) - { + for (int k = 0; k < 3; ++k) { random_value[k] = static_cast(std::rand()) / RAND_MAX - 0.5; } 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) - { + for (int j = 0; j < 3; ++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 48e1185c270..ab0793e9bea 100644 --- a/source/source_md/langevin.h +++ b/source/source_md/langevin.h @@ -13,8 +13,7 @@ * and by imposing a frictional drag force on the motion of the nanoparticle in the solvent. * The damping factor and the random force combine to give the correct NVT ensemble. */ -class Langevin : public MD_base -{ +class Langevin : public MD_base { public: Langevin(const MD_para& mdp_in, const bool cal_stress_in, diff --git a/source/source_md/md_base.cpp b/source/source_md/md_base.cpp index fa7ac21f4a6..b65b5d2ee81 100644 --- a/source/source_md/md_base.cpp +++ b/source/source_md/md_base.cpp @@ -15,8 +15,7 @@ MD_base::MD_base(const MD_para& mdp_in, const bool init_vel, const int my_rank_in, MDCell& mdcell_in) - : mdp(mdp_in), mdcell(mdcell_in) -{ + : mdp(mdp_in), mdcell(mdcell_in) { #ifdef __MPI my_rank = mdcell.mpi_rank(); #else @@ -45,16 +44,12 @@ MD_base::MD_base(const MD_para& mdp_in, t_current = MD_func::current_temp(kinetic, mdcell, frozen_freedom_); } -MD_base::~MD_base() -{ -} +MD_base::~MD_base() {} void MD_base::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, - DomainDecomposition& decomp) -{ - if (mdp.md_restart) - { + DomainDecomposition& decomp) { + if (mdp.md_restart) { restart(global_readin_dir); } @@ -67,44 +62,35 @@ void MD_base::setup(ModuleESolver::ESolver* p_esolver, 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()) - { + if (mdcell.has_backing_unitcell()) { mdcell.backing_unitcell().ionic_position_updated = true; } return; } -void MD_base::first_half(std::ofstream& ofs) -{ +void MD_base::first_half(std::ofstream& ofs) { update_vel(); update_pos(); return; } -void MD_base::second_half() -{ +void MD_base::second_half() { update_vel(); return; } -void MD_base::update_pos() -{ +void MD_base::update_pos() { std::vector& atoms = mdcell.owned_atoms(); - for (std::size_t i = 0; i < atoms.size(); ++i) - { + for (std::size_t i = 0; i < atoms.size(); ++i) { LocalAtom& atom = atoms[i]; ModuleBase::Vector3 pos; - for (int k = 0; k < 3; ++k) - { - if (atom.mbl[k]) - { + for (int k = 0; k < 3; ++k) { + if (atom.mbl[k]) { pos[k] = atom.vel[k] * md_dt / mdcell.lat0(); - } - else - { + } else { pos[k] = 0; } } @@ -119,16 +105,12 @@ void MD_base::update_pos() return; } -void MD_base::update_vel() -{ +void MD_base::update_vel() { std::vector& atoms = mdcell.owned_atoms(); - for (std::size_t i = 0; i < atoms.size(); ++i) - { + for (std::size_t i = 0; i < atoms.size(); ++i) { LocalAtom& atom = atoms[i]; - for (int k = 0; k < 3; ++k) - { - if (atom.mbl[k]) - { + for (int k = 0; k < 3; ++k) { + if (atom.mbl[k]) { atom.vel[k] += 0.5 * atom.force[k] * md_dt / atom.mass; } } @@ -136,12 +118,10 @@ 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; } @@ -149,8 +129,7 @@ void MD_base::print_md(std::ofstream& ofs, const bool& cal_stress) const double unit_transform = ModuleBase::HARTREE_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; double press = 0.0; - for (int i = 0; i < 3; i++) - { + for (int i = 0; i < 3; i++) { press += stress(i, i) / 3; } @@ -161,8 +140,7 @@ void MD_base::print_md(std::ofstream& ofs, const bool& cal_stress) std::cout << " " << std::left << std::setw(20) << "Energy (Ry)" << std::left << std::setw(20) << "Potential (Ry)" << std::left << std::setw(20) << "Kinetic (Ry)" << std::left << std::setw(20) << "Temperature (K)"; - if (cal_stress) - { + if (cal_stress) { std::cout << std::left << std::setw(20) << "Pressure (kbar)"; } @@ -173,8 +151,7 @@ void MD_base::print_md(std::ofstream& ofs, const bool& cal_stress) << t_current * ModuleBase::Hartree_to_K; std::cout << std::setprecision(8); - if (cal_stress) - { + if (cal_stress) { std::cout << std::left << std::setw(20) << press * unit_transform; } @@ -186,8 +163,7 @@ void MD_base::print_md(std::ofstream& ofs, const bool& cal_stress) ofs.unsetf(std::ios::fixed); ofs << std::setprecision(8); - if (cal_stress) - { + if (cal_stress) { MD_func::print_stress(ofs, virial, stress); ofs << std::endl; } @@ -197,8 +173,7 @@ void MD_base::print_md(std::ofstream& ofs, const bool& cal_stress) ofs << " " << std::left << std::setw(20) << "Energy (Ry)" << std::left << std::setw(20) << "Potential (Ry)" << std::left << std::setw(20) << "Kinetic (Ry)" << std::left << std::setw(20) << "Temperature (K)"; - if (cal_stress) - { + if (cal_stress) { ofs << std::left << std::setw(20) << "Pressure (kbar)"; } @@ -209,8 +184,7 @@ void MD_base::print_md(std::ofstream& ofs, const bool& cal_stress) << t_current * ModuleBase::Hartree_to_K; ofs << std::setprecision(8); - if (cal_stress) - { + if (cal_stress) { ofs << std::left << std::setw(20) << press * unit_transform; } @@ -221,10 +195,8 @@ void MD_base::print_md(std::ofstream& ofs, const bool& cal_stress) return; } -void MD_base::write_restart(const std::string& global_out_dir) -{ - if (!my_rank) - { +void MD_base::write_restart(const std::string& global_out_dir) { + if (!my_rank) { std::stringstream ssc; ssc << global_out_dir << "Restart_md.txt"; std::ofstream file(ssc.str().c_str()); @@ -240,8 +212,7 @@ 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); return; diff --git a/source/source_md/md_base.h b/source/source_md/md_base.h index 900392b8d96..5b4b7bc19b1 100644 --- a/source/source_md/md_base.h +++ b/source/source_md/md_base.h @@ -17,8 +17,7 @@ class DomainDecomposition; * energy or particles with its environment, so that the energy of the system * does not change with time. */ -class MD_base -{ +class MD_base { public: /** * @brief construct the integrator from the values it actually uses @@ -43,9 +42,8 @@ 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 @@ -97,10 +95,7 @@ class MD_base public: /// @brief the time increment in a.u., converted from mdp.md_dt - double get_md_dt() const - { - return md_dt; - } + double get_md_dt() const { return md_dt; } bool stop; ///< MD stop or not double t_current; ///< current temperature diff --git a/source/source_md/md_func.cpp b/source/source_md/md_func.cpp index 23e25e0e779..69f9222e2b9 100644 --- a/source/source_md/md_func.cpp +++ b/source/source_md/md_func.cpp @@ -10,18 +10,14 @@ #include #include -namespace MD_func -{ +namespace MD_func { #ifdef __MPI -namespace -{ -bool write_dump_at(const int file, const std::string& data, const MPI_Offset offset) -{ +namespace { +bool write_dump_at(const int file, const std::string& data, const MPI_Offset offset) { std::size_t written = 0; - while (written < data.size()) - { - const ssize_t count - = pwrite(file, data.data() + written, data.size() - written, static_cast(offset + written)); + while (written < data.size()) { + const ssize_t count = + pwrite(file, data.data() + written, data.size() - written, static_cast(offset + written)); if (count <= 0) return false; written += static_cast(count); @@ -31,18 +27,15 @@ bool write_dump_at(const int file, const std::string& data, const MPI_Offset off } // namespace #endif -double gaussrand() -{ +double gaussrand() { static double v1 = 0.0; static double v2 = 0.0; static double S = 0.0; static int phase = 0; double xx = 0.0; - if (phase == 0) - { - do - { + if (phase == 0) { + do { double U1 = static_cast(std::rand()) / RAND_MAX; double U2 = static_cast(std::rand()) / RAND_MAX; @@ -52,9 +45,7 @@ double gaussrand() } while (S >= 1 || S == 0); xx = v1 * sqrt(-2.0 * log(S) / S); - } - else - { + } else { xx = v2 * sqrt(-2.0 * log(S) / S); } @@ -63,13 +54,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) - for (int ion = 0; ion < natom; ++ion) - { + for (int ion = 0; ion < natom; ++ion) { ke += 0.5 * allmass[ion] * vel[ion].norm2(); } @@ -79,13 +68,11 @@ double kinetic_energy(const int& natom, const ModuleBase::Vector3* vel, void compute_stress(const MDCell& mdcell, const bool& cal_stress, const ModuleBase::matrix& virial, - ModuleBase::matrix& stress) -{ + ModuleBase::matrix& stress) { if (!cal_stress) return; ModuleBase::matrix t_vector(3, 3); - for (std::size_t i = 0; i < mdcell.owned_atoms().size(); ++i) - { + 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) @@ -99,24 +86,18 @@ void compute_stress(const MDCell& mdcell, 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) - { - for (int ia = 0; ia < unit_in.atoms[it].na; ++ia) - { + for (int it = 0; it < unit_in.ntype; ++it) { + for (int ia = 0; ia < unit_in.atoms[it].na; ++ia) { vel[iat] = unit_in.atoms[it].vel[ia]; - if (unit_in.atoms[it].mbl[ia].x == 0) - { + if (unit_in.atoms[it].mbl[ia].x == 0) { vel[iat].x = 0; } - if (unit_in.atoms[it].mbl[ia].y == 0) - { + if (unit_in.atoms[it].mbl[ia].y == 0) { vel[iat].y = 0; } - if (unit_in.atoms[it].mbl[ia].z == 0) - { + if (unit_in.atoms[it].mbl[ia].z == 0) { vel[iat].z = 0; } ++iat; @@ -131,21 +112,16 @@ void rescale_vel(const int& natom, const double& temperature, const double* allmass, const int& frozen_freedom, - ModuleBase::Vector3* vel) -{ + ModuleBase::Vector3* vel) { double factor = 0.0; - if (3 * natom == frozen_freedom || temperature == 0) - { + if (3 * natom == frozen_freedom || temperature == 0) { factor = 0; - } - else - { + } else { factor = 0.5 * (3 * natom - frozen_freedom) * temperature / kinetic_energy(natom, vel, allmass); } #pragma omp parallel for schedule(static) if (natom >= 256) - for (int i = 0; i < natom; i++) - { + for (int i = 0; i < natom; i++) { vel[i] = vel[i] * sqrt(factor); } } @@ -157,40 +133,29 @@ void rand_vel(const int& natom, const ModuleBase::Vector3 frozen, const ModuleBase::Vector3* ionmbl, const int& my_rank, - ModuleBase::Vector3* vel) -{ - if (!my_rank) - { + ModuleBase::Vector3* vel) { + if (!my_rank) { double tot_mass = 0; ModuleBase::Vector3 tot_momentum; - for (int i = 0; i < natom; i++) - { + for (int i = 0; i < natom; i++) { tot_mass += allmass[i]; double sigma = sqrt(temperature / allmass[i]); - for (int k = 0; k < 3; ++k) - { - if (ionmbl[i][k] == 0) - { + for (int k = 0; k < 3; ++k) { + if (ionmbl[i][k] == 0) { vel[i][k] = 0; - } - else - { + } else { vel[i][k] = gaussrand() * sigma; } - if (frozen[k] == 0) - { + if (frozen[k] == 0) { tot_momentum[k] += allmass[i] * vel[i][k]; } } } - for (int k = 0; k < 3; ++k) - { - if (frozen[k] == 0) - { - for (int i = 0; i < natom; i++) - { + for (int k = 0; k < 3; ++k) { + if (frozen[k] == 0) { + for (int i = 0; i < natom; i++) { vel[i][k] -= tot_momentum[k] / tot_mass; } } @@ -211,10 +176,8 @@ void init_vel(MDCell& mdcell, const bool& init_vel, const bool& restart, double& temperature, - std::int64_t& frozen_freedom) -{ - if (mdcell.has_backing_unitcell()) - { + std::int64_t& frozen_freedom) { + if (mdcell.has_backing_unitcell()) { const UnitCell& unit = mdcell.backing_unitcell(); std::vector masses(static_cast(unit.nat), 0.0); std::vector> mobile(static_cast(unit.nat)); @@ -229,22 +192,16 @@ void init_vel(MDCell& mdcell, if (frozen.z == 0) ++frozen_freedom; - if (init_vel) - { + 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()); - if (!restart && temperature >= 0.0 && current > 0.0) - { + if (!restart && temperature >= 0.0 && current > 0.0) { rescale_vel(unit.nat, temperature, masses.data(), frozen_freedom, velocities.data()); - } - else if (!restart && temperature < 0.0) - { + } else if (!restart && temperature < 0.0) { temperature = current; } - } - else - { + } else { #ifdef __MPI const int rank = mdcell.mpi_rank(); #else @@ -261,14 +218,12 @@ void init_vel(MDCell& mdcell, } std::vector offsets(static_cast(unit.ntype + 1), 0); - for (int it = 0; it < unit.ntype; ++it) - { + for (int it = 0; it < unit.ntype; ++it) { offsets[static_cast(it + 1)] = offsets[static_cast(it)] + unit.atoms[it].na; } - for (LocalAtom& atom: mdcell.owned_atoms()) - { - atom.vel - = velocities[static_cast(offsets[static_cast(atom.type)] + atom.type_index)]; + for (LocalAtom& atom: mdcell.owned_atoms()) { + atom.vel = + velocities[static_cast(offsets[static_cast(atom.type)] + atom.type_index)]; } return; } @@ -280,8 +235,7 @@ void init_vel(MDCell& mdcell, if (!atoms[i].mbl[k]) ++frozen[k]; #ifdef __MPI - if (mdcell.mpi_size() > 1) - { + if (mdcell.mpi_size() > 1) { MPI_Allreduce(MPI_IN_PLACE, &frozen.x, 3, MPI_INT64_T, MPI_SUM, mdcell.communicator()); } #endif @@ -292,12 +246,10 @@ void init_vel(MDCell& mdcell, ++frozen_freedom; if (!frozen.z) ++frozen_freedom; - if (init_vel) - { + if (init_vel) { double kinetic = 0.0; const double current = current_temp(kinetic, mdcell, frozen_freedom); - if (!restart && current > 0.0 && temperature > 0.0) - { + if (!restart && current > 0.0 && temperature > 0.0) { const double factor = sqrt(temperature / current); for (LocalAtom& atom: atoms) atom.vel *= factor; @@ -307,11 +259,9 @@ void init_vel(MDCell& mdcell, 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) - { + 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]; @@ -321,18 +271,15 @@ void init_vel(MDCell& mdcell, 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 (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; } } double kinetic = 0.0; const double current = current_temp(kinetic, mdcell, frozen_freedom); - if (current > 0.0 && temperature > 0.0) - { + if (current > 0.0 && temperature > 0.0) { const double factor = sqrt(temperature / current); for (LocalAtom& atom: atoms) atom.vel *= factor; @@ -346,39 +293,31 @@ void force_virial(ModuleESolver::ESolver* p_esolver, double& potential, const bool& cal_stress, ModuleBase::matrix& virial, - const bool& md_out_force) -{ + const bool& md_out_force) { ModuleBase::TITLE("MD_func", "force_virial"); ModuleBase::timer::start("MD_func", "force_virial"); - if (!mdcell.has_backing_unitcell()) - { + if (!mdcell.has_backing_unitcell()) { decomp.prepare_neighbors(mdcell); 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; - if (md_out_force) - { + if (md_out_force) { ModuleIO::print_force(GlobalV::ofs_running, mdcell, "TOTAL-FORCE (eV/Angstrom)"); } - if (cal_stress) - { + if (cal_stress) { p_esolver->cal_stress(static_cast(mdcell), virial); virial *= 0.5; } - } - else - { + } else { UnitCell& ucell = mdcell.backing_unitcell(); std::vector>> backing_velocities(static_cast(ucell.ntype)); - for (int it = 0; it < ucell.ntype; ++it) - { + for (int it = 0; it < ucell.ntype; ++it) { backing_velocities[static_cast(it)] = ucell.atoms[it].vel; } mdcell.sync_backing_unitcell(); - for (int it = 0; it < ucell.ntype; ++it) - { + 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()); @@ -388,16 +327,14 @@ void force_virial(ModuleESolver::ESolver* p_esolver, ModuleBase::matrix full_force(ucell.nat, 3); p_esolver->cal_force(ucell, full_force); full_force *= 0.5; - if (cal_stress) - { + 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()) - { + 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)); } @@ -405,13 +342,11 @@ void force_virial(ModuleESolver::ESolver* p_esolver, ModuleBase::timer::end("MD_func", "force_virial"); } -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) { double stress_scalar = 0.0; double virial_scalar = 0.0; - for (int i = 0; i < 3; i++) - { + for (int i = 0; i < 3; i++) { stress_scalar += stress(i, i) / 3.0; virial_scalar += virial(i, i) / 3.0; } @@ -444,8 +379,7 @@ void dump_info(const int& step, const MDCell& mdcell, const MD_para& mdp, const bool cal_stress, - const ModuleBase::matrix& virial) -{ + const ModuleBase::matrix& virial) { std::stringstream file; file << global_out_dir << "MD_dump"; const double unit_pos = mdcell.lat0() / ModuleBase::ANGSTROM_AU; @@ -458,8 +392,7 @@ void dump_info(const int& step, header << " " << mdcell.latvec().e11 << " " << mdcell.latvec().e12 << " " << mdcell.latvec().e13 << "\n"; header << " " << mdcell.latvec().e21 << " " << mdcell.latvec().e22 << " " << mdcell.latvec().e23 << "\n"; header << " " << mdcell.latvec().e31 << " " << mdcell.latvec().e32 << " " << mdcell.latvec().e33 << "\n"; - if (cal_stress && mdp.dump_virial) - { + if (cal_stress && mdp.dump_virial) { header << "VIRIAL (kbar)\n"; for (int i = 0; i < 3; ++i) header << " " << virial(i, 0) * unit_virial << " " << virial(i, 1) * unit_virial << " " @@ -472,14 +405,12 @@ void dump_info(const int& step, 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) - { + 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]; } std::ostringstream local; local << std::fixed << std::setprecision(12); - for (int i = 0; i < mdcell.owned_atoms().size(); ++i) - { + for (int i = 0; i < mdcell.owned_atoms().size(); ++i) { const LocalAtom& atom = mdcell.owned_atoms()[static_cast(i)]; local << " " << type_offsets[static_cast(atom.type)] + atom.type_index << " " << mdcell.type_labels()[static_cast(atom.type)] << " " << atom.cart.x * unit_pos << " " @@ -497,10 +428,9 @@ void dump_info(const int& step, int rank = 0; MPI_Comm_rank(comm, &rank); MPI_Offset base = 0; - if (rank == 0) - { - const int file_descriptor - = open(file.str().c_str(), O_CREAT | O_WRONLY | (step == 0 ? O_TRUNC : O_APPEND), 0666); + if (rank == 0) { + const int file_descriptor = + open(file.str().c_str(), O_CREAT | O_WRONLY | (step == 0 ? O_TRUNC : O_APPEND), 0666); if (file_descriptor < 0) ModuleBase::WARNING_QUIT("MD_func::dump_info", "cannot open MD_dump."); base = lseek(file_descriptor, 0, SEEK_END); @@ -527,26 +457,20 @@ void dump_info(const int& step, void get_mass_mbl(const UnitCell& unit_in, double* allmass, ModuleBase::Vector3& frozen, - ModuleBase::Vector3* ionmbl) -{ + ModuleBase::Vector3* ionmbl) { int ion = 0; frozen.set(0, 0, 0); - for (int it = 0; it < unit_in.ntype; it++) - { - for (int i = 0; i < unit_in.atoms[it].na; i++) - { + for (int it = 0; it < unit_in.ntype; it++) { + for (int i = 0; i < unit_in.atoms[it].na; i++) { 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; } @@ -557,8 +481,7 @@ 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); double delta = static_cast(istep) / nstep; return tfirst + delta * (tlast - tfirst); @@ -568,10 +491,8 @@ double current_temp(double& kinetic, const int& natom, const int& frozen_freedom, const double* allmass, - const ModuleBase::Vector3* vel) -{ - if (3 * natom == frozen_freedom) - { + const ModuleBase::Vector3* vel) { + if (3 * natom == frozen_freedom) { kinetic = 0.0; return 0.0; } @@ -579,11 +500,9 @@ 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) - { + for (std::size_t i = 0; i < mdcell.owned_atoms().size(); ++i) { const LocalAtom& atom = mdcell.owned_atoms()[i]; kinetic += 0.5 * atom.mass * atom.vel.norm2(); } @@ -596,11 +515,9 @@ double current_temp(double& kinetic, const MDCell& mdcell, const std::int64_t& f 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) - { + 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]; @@ -623,8 +540,7 @@ std::int64_t global_dof(const MDCell& mdcell) 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; #ifdef __MPI @@ -632,17 +548,14 @@ void current_md_info(const MDCell& mdcell, const std::string& file_dir, int& md_ #else const int rank = 0; #endif - if (rank == 0) - { + if (rank == 0) { std::stringstream ssc; ssc << file_dir << "Restart_md.txt"; std::ifstream file(ssc.str().c_str()); - if (!file) - { + if (!file) { ok = false; } - if (ok) - { + if (ok) { file >> md_step >> temperature; } } @@ -650,8 +563,7 @@ 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()); #endif - if (!ok) - { + if (!ok) { ModuleBase::WARNING_QUIT("current_md_info", "no Restart_md.txt!"); } #ifdef __MPI diff --git a/source/source_md/md_func.h b/source/source_md/md_func.h index f31d1b93fdc..6b7c22cf48a 100644 --- a/source/source_md/md_func.h +++ b/source/source_md/md_func.h @@ -19,8 +19,7 @@ class DomainDecomposition; * * Namespace MD_func contains several functions used in md. */ -namespace MD_func -{ +namespace MD_func { /** * @brief generate a Gaussian random number diff --git a/source/source_md/md_statistics.h b/source/source_md/md_statistics.h index e7bef175be5..8aa16ebeccb 100644 --- a/source/source_md/md_statistics.h +++ b/source/source_md/md_statistics.h @@ -3,17 +3,14 @@ #include "source_base/matrix.h" -namespace MD_func -{ +namespace MD_func { -struct MDKineticState -{ +struct MDKineticState { double kinetic = 0.0; double temperature = 0.0; }; -struct MDStressState -{ +struct MDStressState { ModuleBase::matrix stress; ModuleBase::matrix temperature_tensor; }; diff --git a/source/source_md/module_plumed/plumed_interface.cpp b/source/source_md/module_plumed/plumed_interface.cpp index 4588e1ce9a2..19caf1bd720 100644 --- a/source/source_md/module_plumed/plumed_interface.cpp +++ b/source/source_md/module_plumed/plumed_interface.cpp @@ -8,19 +8,16 @@ // C API of the PLUMED wrapper library, matching in a // PLUMED installation. The handle is an opaque pointer-sized structure; the // wrapper library takes care of loading the PLUMED kernel. -extern "C" -{ - plumed_handle plumed_create(); - void plumed_cmd(plumed_handle plumed, const char* key, const void* val); - void plumed_finalize(plumed_handle plumed); +extern "C" { +plumed_handle plumed_create(); +void plumed_cmd(plumed_handle plumed, const char* key, const void* val); +void plumed_finalize(plumed_handle plumed); } #endif -namespace ModulePlumed -{ +namespace ModulePlumed { -namespace -{ +namespace { // Conversion factors from ABACUS atomic units to the PLUMED default units // (length: nm, energy: kJ/mol, time: ps, mass: amu). They are declared to // PLUMED once in init() with setMD*Units; afterwards all the arrays handed @@ -31,13 +28,9 @@ 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; this->nat_ = nat; @@ -45,15 +38,13 @@ bool PlumedInterface::init(const std::string& plumed_file, const int nat, const // ABACUS stores atomic masses in atomic units (electron masses), while // PLUMED expects them in amu this->masses_amu_.resize(nat); - for (int i = 0; i < nat; ++i) - { + for (int i = 0; i < nat; ++i) { this->masses_amu_[i] = masses_au[i] * ELECTRON_MASS_TO_AMU; } this->plumed_ = plumed_create(); this->has_plumed_ = true; - if (this->plumed_.p == nullptr) - { + if (this->plumed_.p == nullptr) { ModuleBase::WARNING_QUIT("ModulePlumed::PlumedInterface::init", "could not create a PLUMED object, check the PLUMED installation"); } @@ -91,8 +82,7 @@ bool PlumedInterface::init(const std::string& plumed_file, const int nat, const #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; // when the MD run does not provide a virial, a zero buffer is registered @@ -102,19 +92,16 @@ void PlumedInterface::prepare_virial_(const double* cell, const double* virial) // Quantum ESPRESSO interface); ABACUS keeps the potential part of the // stress (Hartree/Bohr^3) in its MD virial matrix, so the volume is // multiplied back here. - for (int i = 0; i < 9; ++i) - { + for (int i = 0; i < 9; ++i) { this->virial_buf_[i] = 0.0; } - if (virial == nullptr || cell == nullptr) - { + if (virial == nullptr || cell == nullptr) { 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])); - for (int i = 0; i < 9; ++i) - { + 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; } } @@ -124,8 +111,7 @@ void PlumedInterface::compute_forces_(const int istep, double* force, const double* cell, const double potential, - const double* virial) -{ + const double* virial) { #ifdef USE_PLUMED int step = istep; double energy = potential; @@ -166,11 +152,9 @@ void PlumedInterface::move(const int istep, ModuleBase::Vector3* force, const double* cell, const double potential, - const double* virial) -{ + const double* virial) { #ifdef USE_PLUMED - if (!this->active_run_) - { + if (!this->active_run_) { return; } this->compute_forces_(istep, pos, reinterpret_cast(force), cell, potential, virial); @@ -184,11 +168,9 @@ void PlumedInterface::move(const int istep, #endif } -void PlumedInterface::finalize() -{ +void PlumedInterface::finalize() { #ifdef USE_PLUMED - if (this->has_plumed_) - { + if (this->has_plumed_) { plumed_finalize(this->plumed_); this->has_plumed_ = false; } diff --git a/source/source_md/module_plumed/plumed_interface.h b/source/source_md/module_plumed/plumed_interface.h index 5749781fe2c..c5019d7d166 100644 --- a/source/source_md/module_plumed/plumed_interface.h +++ b/source/source_md/module_plumed/plumed_interface.h @@ -9,13 +9,11 @@ /// 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 -{ +struct plumed_handle { void* p; }; -namespace ModulePlumed -{ +namespace ModulePlumed { /** * @brief interface between an ABACUS MD run and the PLUMED plugin @@ -51,8 +49,7 @@ namespace ModulePlumed * ``-DENABLE_PLUMED=ON``; requesting PLUMED in INPUT while ABACUS was built * without it stops the run with an error message. */ -class PlumedInterface -{ +class PlumedInterface { public: PlumedInterface() = default; ~PlumedInterface(); @@ -99,10 +96,7 @@ class PlumedInterface 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 diff --git a/source/source_md/msst.cpp b/source/source_md/msst.cpp index 03fb73aa1be..16b7247dd5e 100644 --- a/source/source_md/msst.cpp +++ b/source/source_md/msst.cpp @@ -12,8 +12,7 @@ MSST::MSST(const MD_para& mdp_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) -{ + : 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; msst_vis = mdp.msst_vis / ModuleBase::AU_to_MASS / ModuleBase::ANGSTROM_AU * ModuleBase::AU_to_FS; @@ -36,32 +35,26 @@ MSST::MSST(const MD_para& mdp_in, #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"); MD_base::setup(p_esolver, global_readin_dir, decomp); - if (mdcell.has_backing_unitcell()) - { + if (mdcell.has_backing_unitcell()) { mdcell.backing_unitcell().cell_parameter_updated = true; } int sd = mdp.msst_direction; - if (!mdp.md_restart) - { + if (!mdp.md_restart) { lag_pos = 0; v0 = mdcell.omega(); p0 = stress(sd, sd); e0 = potential + kinetic; - if (kinetic > 0 && mdp.msst_tscale > 0) - { + if (kinetic > 0 && mdp.msst_tscale > 0) { double fac1 = mdp.msst_tscale * totmass * 2.0 * kinetic / msst_qmass; omega[sd] = -1.0 * sqrt(fac1); double fac2 = omega[sd] / v0; @@ -81,8 +74,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"); @@ -98,8 +90,7 @@ void MSST::first_half(std::ofstream& ofs) /// save the velocities old_v.resize(mdcell.owned_atoms().size()); - for (int i = 0; i < mdcell.owned_atoms().size(); ++i) - { + for (int i = 0; i < mdcell.owned_atoms().size(); ++i) { old_v[static_cast(i)] = mdcell.owned_atoms()[static_cast(i)].vel; } @@ -109,8 +100,7 @@ void MSST::first_half(std::ofstream& ofs) vsum = vel_sum(); /// reset the velocities - for (int i = 0; i < mdcell.owned_atoms().size(); ++i) - { + for (int i = 0; i < mdcell.owned_atoms().size(); ++i) { mdcell.owned_atoms()[static_cast(i)].vel = old_v[static_cast(i)]; } @@ -137,8 +127,7 @@ void MSST::first_half(std::ofstream& ofs) return; } -void MSST::second_half() -{ +void MSST::second_half() { ModuleBase::TITLE("MSST", "second_half"); ModuleBase::timer::start("MSST", "second_half"); @@ -164,17 +153,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) -{ - if (!my_rank) - { +void MSST::write_restart(const std::string& global_out_dir) { + if (!my_rank) { std::stringstream ssc; ssc << global_out_dir << "Restart_md.txt"; std::ofstream file(ssc.str().c_str()); @@ -196,26 +182,21 @@ 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; - if (!my_rank) - { + if (!my_rank) { std::stringstream ssc; ssc << global_readin_dir << "Restart_md.txt"; std::ifstream file(ssc.str().c_str()); - if (!file) - { + if (!file) { ok = false; } - if (ok) - { + if (ok) { file >> step_rst_ >> md_tfirst >> omega[mdp.msst_direction] >> e0 >> v0 >> p0 >> lag_pos; - if (!file) - { + if (!file) { ok = false; } file.close(); @@ -226,8 +207,7 @@ void MSST::restart(const std::string& global_readin_dir) MPI_Bcast(&ok, 1, MPI_C_BOOL, 0, mdcell.communicator()); #endif - if (!ok) - { + if (!ok) { ModuleBase::WARNING_QUIT("mdrun", "no Restart_md.txt !"); } @@ -244,8 +224,7 @@ void MSST::restart(const std::string& global_readin_dir) return; } -double MSST::vel_sum() const -{ +double MSST::vel_sum() const { double vsum = 0; for (const LocalAtom& atom: mdcell.owned_atoms()) @@ -256,8 +235,7 @@ double MSST::vel_sum() const 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; assert(mdcell.omega() > 0.0); @@ -276,39 +254,32 @@ void MSST::rescale(std::ofstream& ofs, const double& volume) static_cast(ofs); } -void MSST::propagate_vel() -{ +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; const_D.set(fac / atom.mass, fac / atom.mass, fac / atom.mass); const_D[sd] -= 2 * omega[sd] / mdcell.omega(); - for (int k = 0; k < 3; ++k) - { - if (fabs(dthalf * const_D[k]) > 1e-6) - { + for (int k = 0; k < 3; ++k) { + if (fabs(dthalf * const_D[k]) > 1e-6) { double expd = exp(dthalf * const_D[k]); atom.vel[k] = expd * (const_C[k] + const_D[k] * atom.vel[k] - const_C[k] / expd) / const_D[k]; - } - 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; + } 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; } } } return; } -void MSST::propagate_voldot() -{ +void MSST::propagate_voldot() { const int sd = mdp.msst_direction; const double dthalf = 0.5 * md_dt; double p_current = stress(sd, sd); @@ -317,19 +288,15 @@ void MSST::propagate_voldot() double const_B = totmass * msst_vis / (msst_qmass * mdcell.omega()); /// prevent the increase of volume - if (mdcell.omega() > v0 && const_A > 0) - { + if (mdcell.omega() > v0 && const_A > 0) { const_A = -const_A; } /// avoid singularity at B = 0 with Taylor expansion double fac = const_B * dthalf; - if (fac > 1e-6) - { + if (fac > 1e-6) { omega[sd] = (omega[sd] + const_A * (exp(fac) - 1) / const_B) * exp(-fac); - } - else - { + } else { omega[sd] += (const_A - const_B * omega[sd]) * dthalf + 0.5 * (const_B * const_B * omega[sd] - const_A * const_B) * dthalf * dthalf; } diff --git a/source/source_md/msst.h b/source/source_md/msst.h index f230522141c..9845ecb6109 100644 --- a/source/source_md/msst.h +++ b/source/source_md/msst.h @@ -13,41 +13,25 @@ * The MSST varies the cell volume and temperature to restrain the system to the shock Hugoniot and the Rayleigh line. * These restraints correspond to the macroscopic conservation laws dictated by a shock front. */ -class MSST : public MD_base -{ +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(); /// @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); diff --git a/source/source_md/nhchain.cpp b/source/source_md/nhchain.cpp index 4c4bd3a8502..7fe7b4a1647 100644 --- a/source/source_md/nhchain.cpp +++ b/source/source_md/nhchain.cpp @@ -11,8 +11,7 @@ Nose_Hoover::Nose_Hoover(const MD_para& mdp_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) -{ + : 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; @@ -23,31 +22,24 @@ Nose_Hoover::Nose_Hoover(const MD_para& mdp_in, md_plast = mdp.md_plast / unit_transform; md_pfreq = mdp.md_pfreq * ModuleBase::AU_to_FS; - if (md_tfirst == 0) - { + if (md_tfirst == 0) { ModuleBase::WARNING_QUIT("Nose_Hoover", " md_tfirst must be larger than 0 in NHC"); } /// init NPT related variables - for (int i = 0; i < 6; ++i) - { + for (int i = 0; i < 6; ++i) { pstart[i] = pstop[i] = pfreq[i] = p_target[i] = pflag[i] = 0; } - if (mdp.md_type == "npt") - { + if (mdp.md_type == "npt") { /// determine the NPT methods - if (mdp.md_pmode == "iso") - { + if (mdp.md_pmode == "iso") { pstart[0] = pstart[1] = pstart[2] = md_pfirst; pstop[0] = pstop[1] = pstop[2] = md_plast; pfreq[0] = pfreq[1] = pfreq[2] = md_pfreq; pflag[0] = pflag[1] = pflag[2] = 1; - } - else if (mdp.md_pmode == "aniso") - { - if (mdp.md_pcouple == "xyz") - { + } else if (mdp.md_pmode == "aniso") { + if (mdp.md_pcouple == "xyz") { ModuleBase::WARNING_QUIT("Nose_Hoover", "md_pcouple==xyz will convert aniso to iso!"); } pstart[0] = pstart[1] = pstart[2] = md_pfirst; @@ -62,11 +54,9 @@ Nose_Hoover::Nose_Hoover(const MD_para& mdp_in, * e31 e32 e33 * Under Voigt notation, xx, yy, zz, yz, xz, xy. */ - else if (mdp.md_pmode == "tri") - { + else if (mdp.md_pmode == "tri") { const ModuleBase::Matrix3& latvec = mdcell.latvec(); - if (latvec.e12 || latvec.e13 || latvec.e23) - { + if (latvec.e12 || latvec.e13 || latvec.e23) { ModuleBase::WARNING_QUIT("Nose_Hoover", "the lattice must be lower-triangular when md_pmode == tri!"); } pstart[0] = pstart[1] = pstart[2] = md_pfirst; @@ -78,17 +68,14 @@ Nose_Hoover::Nose_Hoover(const MD_para& mdp_in, pstop[3] = pstop[4] = pstop[5] = 0; pfreq[3] = pfreq[4] = pfreq[5] = md_pfreq; pflag[3] = pflag[4] = pflag[5] = 1; - } - else - { + } else { ModuleBase::WARNING_QUIT("Nose_Hoover", "No such md_pmode yet!"); } } /// determine whether NPT ensemble npt_flag = 0; - for (int i = 0; i < 6; ++i) - { + for (int i = 0; i < 6; ++i) { npt_flag += pflag[i]; } pdim = pflag[0] + pflag[1] + pflag[2]; @@ -104,29 +91,24 @@ Nose_Hoover::Nose_Hoover(const MD_para& mdp_in, g_eta = new double[mdp.md_tchain]; v_eta[mdp.md_tchain] = 0; - for (int i = 0; i < mdp.md_tchain; ++i) - { + for (int i = 0; i < mdp.md_tchain; ++i) { eta[i] = v_eta[i] = g_eta[i] = 0; } /// allocate barostat and thermostats coupled with barostat - if (npt_flag) - { - for (int i = 0; i < 6; ++i) - { + if (npt_flag) { + for (int i = 0; i < 6; ++i) { v_omega[i] = mass_omega[i] = 0; } - if (mdp.md_pchain) - { + if (mdp.md_pchain) { mass_peta = new double[mdp.md_pchain]; peta = new double[mdp.md_pchain]; v_peta = new double[mdp.md_pchain + 1]; g_peta = new double[mdp.md_pchain]; v_peta[mdp.md_pchain] = 0; - for (int i = 0; i < mdp.md_pchain; ++i) - { + for (int i = 0; i < mdp.md_pchain; ++i) { peta[i] = v_peta[i] = g_peta[i] = 0; } } @@ -143,15 +125,13 @@ Nose_Hoover::Nose_Hoover(const MD_para& mdp_in, w[3] = 1 - w[0] - w[1] - w[2] - w[4] - w[5] - w[6]; } -Nose_Hoover::~Nose_Hoover() -{ +Nose_Hoover::~Nose_Hoover() { delete[] mass_eta; delete[] eta; delete[] v_eta; delete[] g_eta; - if (npt_flag && mdp.md_pchain) - { + if (npt_flag && mdp.md_pchain) { delete[] mass_peta; delete[] peta; delete[] v_peta; @@ -161,8 +141,7 @@ Nose_Hoover::~Nose_Hoover() void Nose_Hoover::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_readin_dir, - DomainDecomposition& decomp) -{ + DomainDecomposition& decomp) { ModuleBase::TITLE("Nose_Hoover", "setup"); ModuleBase::timer::start("Nose_Hoover", "setup"); @@ -173,15 +152,13 @@ void Nose_Hoover::setup(ModuleESolver::ESolver* p_esolver, /// init thermostats coupled with particles mass_eta[0] = tdof * t_target / md_tfreq / md_tfreq; - for (int m = 1; m < mdp.md_tchain; ++m) - { + for (int m = 1; m < mdp.md_tchain; ++m) { mass_eta[m] = t_target / md_tfreq / md_tfreq; g_eta[m] = (mass_eta[m - 1] * v_eta[m - 1] * v_eta[m - 1] - t_target) / mass_eta[m]; } /// NPT ensemble - if (npt_flag) - { + if (npt_flag) { /// determine target stress target_stress(); @@ -191,20 +168,16 @@ void Nose_Hoover::setup(ModuleESolver::ESolver* p_esolver, /// init barostat double nkt = (static_cast(mdcell.nat()) + 1.0) * t_target; - for (int i = 0; i < 6; ++i) - { - if (pflag[i]) - { + for (int i = 0; i < 6; ++i) { + if (pflag[i]) { mass_omega[i] = nkt / pfreq[i] / pfreq[i]; } } /// init thermostats coupled with barostat - if (mdp.md_pchain) - { + if (mdp.md_pchain) { mass_peta[0] = t_target / md_pfreq / md_pfreq; - for (int m = 1; m < mdp.md_pchain; ++m) - { + for (int m = 1; m < mdp.md_pchain; ++m) { mass_peta[m] = t_target / md_pfreq / md_pfreq; g_peta[m] = (mass_peta[m - 1] * v_peta[m - 1] * v_peta[m - 1] - t_target) / mass_peta[m]; } @@ -216,14 +189,12 @@ void Nose_Hoover::setup(ModuleESolver::ESolver* p_esolver, 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"); /// update thermostats coupled with barostat if NPT ensemble - if (npt_flag && mdp.md_pchain) - { + if (npt_flag && mdp.md_pchain) { baro_thermo(); } @@ -233,8 +204,7 @@ void Nose_Hoover::first_half(std::ofstream& ofs) /// update thermostats coupled with particles particle_thermo(); - if (npt_flag) - { + if (npt_flag) { /// 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); @@ -255,8 +225,7 @@ void Nose_Hoover::first_half(std::ofstream& ofs) /// perform half-step update of vel due to atomic force MD_base::update_vel(); - if (npt_flag) - { + if (npt_flag) { /// perform half-step update of volume update_volume(ofs); } @@ -264,8 +233,7 @@ void Nose_Hoover::first_half(std::ofstream& ofs) /// perform one step update of pos due to atomic velocity MD_base::update_pos(); - if (npt_flag) - { + if (npt_flag) { /// perform half-step update of volume update_volume(ofs); } @@ -275,16 +243,14 @@ void Nose_Hoover::first_half(std::ofstream& ofs) return; } -void Nose_Hoover::second_half() -{ +void Nose_Hoover::second_half() { ModuleBase::TITLE("Nose_Hoover", "second_half"); ModuleBase::timer::start("Nose_Hoover", "second_half"); /// perform half-step update of vel due to atomic force MD_base::update_vel(); - if (npt_flag) - { + if (npt_flag) { /// update vel due to barostat vel_baro(); } @@ -292,8 +258,7 @@ void Nose_Hoover::second_half() /// update temperature and kinetic energy due to velocity rescaling t_current = MD_func::current_temp(kinetic, mdcell, frozen_freedom_); - if (npt_flag) - { + if (npt_flag) { /// update stress due to velocity rescaling MD_func::compute_stress(mdcell, cal_stress, virial, stress); @@ -308,8 +273,7 @@ void Nose_Hoover::second_half() particle_thermo(); /// update thermostats coupled with barostat if NPT ensemble - if (npt_flag && mdp.md_pchain) - { + if (npt_flag && mdp.md_pchain) { baro_thermo(); } @@ -318,15 +282,10 @@ 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) -{ - if (!my_rank) - { +void Nose_Hoover::write_restart(const std::string& global_out_dir) { + if (!my_rank) { std::stringstream ssc; ssc << global_out_dir << "Restart_md.txt"; std::ofstream file(ssc.str().c_str()); @@ -334,34 +293,28 @@ void Nose_Hoover::write_restart(const std::string& global_out_dir) file << step_ + step_rst_ << std::endl; file << md_tfirst << std::endl; file << mdp.md_tchain << std::endl; - for (int i = 0; i < mdp.md_tchain; ++i) - { + for (int i = 0; i < mdp.md_tchain; ++i) { file << eta[i] << " "; } file << std::endl; - for (int i = 0; i < mdp.md_tchain; ++i) - { + for (int i = 0; i < mdp.md_tchain; ++i) { file << v_eta[i] << " "; } file << std::endl; /// npt - if (npt_flag) - { - for (int i = 0; i < 6; ++i) - { + if (npt_flag) { + for (int i = 0; i < 6; ++i) { file << v_omega[i] << " "; } file << std::endl; file << mdp.md_pchain << std::endl; - for (int i = 0; i < mdp.md_pchain; ++i) - { + for (int i = 0; i < mdp.md_pchain; ++i) { file << peta[i] << " "; } file << std::endl; - for (int i = 0; i < mdp.md_pchain; ++i) - { + for (int i = 0; i < mdp.md_pchain; ++i) { file << v_peta[i] << " "; } file << std::endl; @@ -373,67 +326,53 @@ 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; bool ok3 = true; - if (!my_rank) - { + if (!my_rank) { std::stringstream ssc; ssc << global_readin_dir << "Restart_md.txt"; std::ifstream file(ssc.str().c_str()); - if (!file) - { + if (!file) { ok = false; } - if (ok) - { + if (ok) { double Mnum = 0.0; file >> step_rst_ >> md_tfirst >> Mnum; - if (Mnum != mdp.md_tchain) - { + if (Mnum != mdp.md_tchain) { ok2 = false; } - if (ok2) - { - for (int i = 0; i < mdp.md_tchain; ++i) - { + if (ok2) { + for (int i = 0; i < mdp.md_tchain; ++i) { file >> eta[i]; } - for (int i = 0; i < mdp.md_tchain; ++i) - { + for (int i = 0; i < mdp.md_tchain; ++i) { file >> v_eta[i]; } } /// npt - if (npt_flag) - { - for (int i = 0; i < 6; ++i) - { + if (npt_flag) { + for (int i = 0; i < 6; ++i) { file >> v_omega[i]; } file >> Mnum; - if (Mnum != mdp.md_pchain) - { + if (Mnum != mdp.md_pchain) { ok3 = false; } - if (ok3) - { - for (int i = 0; i < mdp.md_pchain; ++i) - { + if (ok3) { + for (int i = 0; i < mdp.md_pchain; ++i) { file >> peta[i]; } - for (int i = 0; i < mdp.md_pchain; ++i) - { + for (int i = 0; i < mdp.md_pchain; ++i) { file >> v_peta[i]; } } @@ -449,16 +388,13 @@ void Nose_Hoover::restart(const std::string& global_readin_dir) MPI_Bcast(&ok3, 1, MPI_C_BOOL, 0, mdcell.communicator()); #endif - if (!ok) - { + if (!ok) { ModuleBase::WARNING_QUIT("Nose_Hoover", "no Restart_md.txt !"); } - if (!ok2) - { + if (!ok2) { ModuleBase::WARNING_QUIT("Nose_Hoover", "Num of thermostats coupled with particles is not the same !"); } - if (!ok3) - { + if (!ok3) { ModuleBase::WARNING_QUIT("Nose_Hoover", "Num of thermostats coupled with barostat is not the same !"); } @@ -467,8 +403,7 @@ void Nose_Hoover::restart(const std::string& global_readin_dir) 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) - { + 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()); @@ -476,22 +411,17 @@ void Nose_Hoover::restart(const std::string& global_readin_dir) #endif } -void Nose_Hoover::particle_thermo() -{ +void Nose_Hoover::particle_thermo() { /// update mass_eta mass_eta[0] = tdof * t_target / md_tfreq / md_tfreq; - for (int m = 1; m < mdp.md_tchain; ++m) - { + for (int m = 1; m < mdp.md_tchain; ++m) { mass_eta[m] = t_target / md_tfreq / md_tfreq; } /// propogate g_eta - if (mass_eta[0] > 0) - { + if (mass_eta[0] > 0) { g_eta[0] = (2 * kinetic - tdof * t_target) / mass_eta[0]; - } - else - { + } else { g_eta[0] = 0; } @@ -499,15 +429,12 @@ void Nose_Hoover::particle_thermo() double factor = 0.0; double scale = 1.0; double KE = kinetic; - for (int i = 0; i < nc_tchain; ++i) - { - for (int j = 0; j < nys; ++j) - { + for (int i = 0; i < nc_tchain; ++i) { + for (int j = 0; j < nys; ++j) { double delta = w[j] * md_dt / nc_tchain; /// propogate v_eta - for (int m = mdp.md_tchain - 1; m >= 0; --m) - { + for (int m = mdp.md_tchain - 1; m >= 0; --m) { factor = exp(-v_eta[m + 1] * delta / 8.0); v_eta[m] *= factor; v_eta[m] += g_eta[m] * delta / 4.0; @@ -515,26 +442,21 @@ void Nose_Hoover::particle_thermo() } /// propogate eta - for (int m = 0; m < mdp.md_tchain; ++m) - { + for (int m = 0; m < mdp.md_tchain; ++m) { eta[m] += v_eta[m] * delta / 2.0; } /// update rescale factor of particle velocity scale *= exp(-v_eta[0] * delta / 2.0); - if (!std::isfinite(scale)) - { + if (!std::isfinite(scale)) { ModuleBase::WARNING_QUIT("Nose_Hoover", "Please set a proper md_tfreq!"); } KE = kinetic * scale * scale; /// propogate g_eta - if (mass_eta[0] > 0) - { + if (mass_eta[0] > 0) { g_eta[0] = (2 * KE - tdof * t_target) / mass_eta[0]; - } - else - { + } else { g_eta[0] = 0; } @@ -543,8 +465,7 @@ void Nose_Hoover::particle_thermo() v_eta[0] += g_eta[0] * delta / 4.0; v_eta[0] *= factor; - for (int m = 1; m < mdp.md_tchain; ++m) - { + for (int m = 1; m < mdp.md_tchain; ++m) { factor = exp(-v_eta[m + 1] * delta / 8.0); v_eta[m] *= factor; g_eta[m] = (mass_eta[m - 1] * v_eta[m - 1] * v_eta[m - 1] - t_target) / mass_eta[m]; @@ -559,25 +480,21 @@ void Nose_Hoover::particle_thermo() atom.vel *= scale; } -void Nose_Hoover::baro_thermo() -{ +void Nose_Hoover::baro_thermo() { /// the freedom of lattice int pdof = npt_flag; /// update kenetic energy of lattice double ke_omega = 0; - for (int i = 0; i < 6; ++i) - { - if (pflag[i]) - { + for (int i = 0; i < 6; ++i) { + if (pflag[i]) { ke_omega += mass_omega[i] * v_omega[i] * v_omega[i]; } } /// update force double lkt_press = t_target; - if (mdp.md_pmode != "iso") - { + if (mdp.md_pmode != "iso") { lkt_press *= pdof; } g_peta[0] = (ke_omega - lkt_press) / mass_peta[0]; @@ -586,15 +503,12 @@ void Nose_Hoover::baro_thermo() double factor = 0.0; double scale = 1.0; double kecurrent = ke_omega; - for (int i = 0; i < nc_pchain; ++i) - { - for (int j = 0; j < nys; ++j) - { + for (int i = 0; i < nc_pchain; ++i) { + for (int j = 0; j < nys; ++j) { double delta = w[j] * md_dt / nc_pchain; /// propogate v_peta - for (int m = mdp.md_pchain - 1; m >= 0; --m) - { + for (int m = mdp.md_pchain - 1; m >= 0; --m) { factor = exp(-v_peta[m + 1] * delta / 8.0); v_peta[m] *= factor; v_peta[m] += g_peta[m] * delta / 4.0; @@ -602,8 +516,7 @@ void Nose_Hoover::baro_thermo() } /// propogate peta - for (int m = 0; m < mdp.md_pchain; ++m) - { + for (int m = 0; m < mdp.md_pchain; ++m) { peta[m] += v_peta[m] * delta / 2.0; } @@ -619,8 +532,7 @@ void Nose_Hoover::baro_thermo() v_peta[0] += g_peta[0] * delta / 4.0; v_peta[0] *= factor; - for (int m = 1; m < mdp.md_pchain; ++m) - { + for (int m = 1; m < mdp.md_pchain; ++m) { factor = exp(-v_peta[m + 1] * delta / 8.0); v_peta[m] *= factor; g_peta[m] = (mass_peta[m - 1] * v_peta[m - 1] * v_peta[m - 1] - t_target) / mass_peta[m]; @@ -631,24 +543,18 @@ void Nose_Hoover::baro_thermo() } /// rescale lattice due to thermostats - for (int i = 0; i < 6; ++i) - { - if (pflag[i]) - { + for (int i = 0; i < 6; ++i) { + if (pflag[i]) { v_omega[i] *= scale; } } } -void Nose_Hoover::update_baro() -{ +void Nose_Hoover::update_baro() { double term_one = 0; - if (mdp.md_pmode == "iso") - { + if (mdp.md_pmode == "iso") { term_one = tdof * t_current; - } - else - { + } else { ModuleBase::matrix t_vector(3, 3); for (const LocalAtom& atom: mdcell.owned_atoms()) for (int i = 0; i < 3; ++i) @@ -658,10 +564,8 @@ void Nose_Hoover::update_baro() 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) - { - if (pflag[i]) - { + for (int i = 0; i < 3; ++i) { + if (pflag[i]) { term_one += t_vector(i, i); } } @@ -670,10 +574,8 @@ void Nose_Hoover::update_baro() double g_omega = 0.0; double term_two = 0; - for (int i = 0; i < 3; ++i) - { - if (pflag[i]) - { + for (int i = 0; i < 3; ++i) { + if (pflag[i]) { g_omega = (p_current[i] - p_hydro) * mdcell.omega() / mass_omega[i] + term_one / mass_omega[i]; v_omega[i] += g_omega * md_dt / 2.0; term_two += v_omega[i]; @@ -681,10 +583,8 @@ void Nose_Hoover::update_baro() } term_two /= static_cast(pdim) * mdcell.nat(); - for (int i = 3; i < 6; ++i) - { - if (pflag[i]) - { + for (int i = 3; i < 6; ++i) { + if (pflag[i]) { g_omega = p_current[i] * mdcell.omega() / mass_omega[i]; v_omega[i] += g_omega * md_dt / 2.0; } @@ -693,70 +593,58 @@ void Nose_Hoover::update_baro() mtk_term = term_two; } -void Nose_Hoover::vel_baro() -{ +void Nose_Hoover::vel_baro() { double factor[3] = {0.0}; - for (int i = 0; i < 3; ++i) - { + for (int i = 0; i < 3; ++i) { factor[i] = exp(-(v_omega[i] + mtk_term) * md_dt / 4); } - for (LocalAtom& atom: mdcell.owned_atoms()) - { - for (int j = 0; j < 3; ++j) - { + for (LocalAtom& atom: mdcell.owned_atoms()) { + for (int j = 0; j < 3; ++j) { atom.vel[j] *= factor[j]; } /// Note: I am not sure whether fixed atoms should update here - if (atom.mbl[0]) - { + if (atom.mbl[0]) { atom.vel[0] -= (atom.vel[1] * v_omega[5] + atom.vel[2] * v_omega[4]) * md_dt / 2; } - if (atom.mbl[1]) - { + if (atom.mbl[1]) { atom.vel[1] -= atom.vel[2] * v_omega[3] * md_dt / 2; } - for (int j = 0; j < 3; ++j) - { + for (int j = 0; j < 3; ++j) { atom.vel[j] *= factor[j]; } } } -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(); /// tri mode, off-diagonal components, first half - if (pflag[4]) - { + if (pflag[4]) { factor = exp(v_omega[0] * md_dt / 16); latvec.e31 *= factor; latvec.e31 += (v_omega[5] * latvec.e32 + v_omega[4] * latvec.e33); latvec.e31 *= factor; } - if (pflag[3]) - { + if (pflag[3]) { factor = exp(v_omega[1] * md_dt / 8); latvec.e32 *= factor; latvec.e32 += (v_omega[3] * latvec.e33); latvec.e32 *= factor; } - if (pflag[5]) - { + if (pflag[5]) { factor = exp(v_omega[0] * md_dt / 8); latvec.e21 *= factor; latvec.e21 += (v_omega[5] * latvec.e22); latvec.e21 *= factor; } - if (pflag[4]) - { + if (pflag[4]) { factor = exp(v_omega[0] * md_dt / 16); latvec.e31 *= factor; latvec.e31 += (v_omega[5] * latvec.e32 + v_omega[4] * latvec.e33); @@ -764,51 +652,44 @@ void Nose_Hoover::update_volume(std::ofstream& ofs) } /// Diagonal components - if (pflag[0]) - { + if (pflag[0]) { factor = exp(v_omega[0] * md_dt / 2); latvec.e11 *= factor; } - if (pflag[1]) - { + if (pflag[1]) { factor = exp(v_omega[1] * md_dt / 2); latvec.e22 *= factor; } - if (pflag[2]) - { + if (pflag[2]) { factor = exp(v_omega[2] * md_dt / 2); latvec.e33 *= factor; } /// tri mode, off-diagonal components, second half - if (pflag[4]) - { + if (pflag[4]) { factor = exp(v_omega[0] * md_dt / 16); latvec.e31 *= factor; latvec.e31 += (v_omega[5] * latvec.e32 + v_omega[4] * latvec.e33); latvec.e31 *= factor; } - if (pflag[3]) - { + if (pflag[3]) { factor = exp(v_omega[1] * md_dt / 8); latvec.e32 *= factor; latvec.e32 += (v_omega[3] * latvec.e33); latvec.e32 *= factor; } - if (pflag[5]) - { + if (pflag[5]) { factor = exp(v_omega[0] * md_dt / 8); latvec.e21 *= factor; latvec.e21 += (v_omega[5] * latvec.e22); latvec.e21 *= factor; } - if (pflag[4]) - { + if (pflag[4]) { factor = exp(v_omega[0] * md_dt / 16); latvec.e31 *= factor; latvec.e31 += (v_omega[5] * latvec.e32 + v_omega[4] * latvec.e33); @@ -820,67 +701,50 @@ void Nose_Hoover::update_volume(std::ofstream& ofs) static_cast(ofs); } -void Nose_Hoover::target_stress() -{ +void Nose_Hoover::target_stress() { double delta = static_cast(step_ + step_rst_) / mdp.md_nstep; p_hydro = 0; - for (int i = 0; i < 3; ++i) - { - if (pflag[i]) - { + for (int i = 0; i < 3; ++i) { + if (pflag[i]) { p_target[i] = pstart[i] + delta * (pstop[i] - pstart[i]); p_hydro += p_target[i]; } } - if (pdim) - { + if (pdim) { p_hydro /= pdim; } - for (int i = 3; i < 6; ++i) - { - if (pflag[i]) - { + for (int i = 3; i < 6; ++i) { + if (pflag[i]) { p_target[i] = pstart[i] + delta * (pstop[i] - pstart[i]); } } } -void Nose_Hoover::couple_stress() -{ - if (mdp.md_pcouple == "xyz") - { +void Nose_Hoover::couple_stress() { + if (mdp.md_pcouple == "xyz") { double ave = (stress(0, 0) + stress(1, 1) + stress(2, 2)) / 3.0; p_current[0] = p_current[1] = p_current[2] = ave; - } - else if (mdp.md_pcouple == "xy") - { + } else if (mdp.md_pcouple == "xy") { double ave = (stress(0, 0) + stress(1, 1)) / 2.0; p_current[0] = p_current[1] = ave; p_current[2] = stress(2, 2); - } - else if (mdp.md_pcouple == "yz") - { + } else if (mdp.md_pcouple == "yz") { double ave = (stress(1, 1) + stress(2, 2)) / 2.0; p_current[1] = p_current[2] = ave; p_current[0] = stress(0, 0); - } - else if (mdp.md_pcouple == "xz") - { + } else if (mdp.md_pcouple == "xz") { double ave = (stress(0, 0) + stress(2, 2)) / 2.0; p_current[0] = p_current[2] = ave; p_current[1] = stress(1, 1); - } - else - { + } else { p_current[0] = stress(0, 0); p_current[1] = stress(1, 1); p_current[2] = stress(2, 2); } - if (!std::isfinite(p_current[0]) || !std::isfinite(p_current[1]) || !std::isfinite(p_current[2])) - { + if (!std::isfinite(p_current[0]) || !std::isfinite(p_current[1]) || !std::isfinite(p_current[2])) { ModuleBase::WARNING_QUIT("Nose_Hoover", "Non-numeric stress component!"); } @@ -888,8 +752,7 @@ void Nose_Hoover::couple_stress() p_current[4] = stress(0, 2); p_current[5] = stress(0, 1); - if (!std::isfinite(p_current[3]) || !std::isfinite(p_current[4]) || !std::isfinite(p_current[5])) - { + if (!std::isfinite(p_current[3]) || !std::isfinite(p_current[4]) || !std::isfinite(p_current[5])) { ModuleBase::WARNING_QUIT("Nose_Hoover", "Non-numeric stress component!"); } } diff --git a/source/source_md/nhchain.h b/source/source_md/nhchain.h index ee8330e7f2b..750734f8593 100644 --- a/source/source_md/nhchain.h +++ b/source/source_md/nhchain.h @@ -9,8 +9,7 @@ * Nose-Hoover style non-Hamiltonian equations of motion. [Phys. Rev. A 31, 1695 (1985)] * Designed to generate positions and velocities sampled from NVT and NPT ensemble. */ -class Nose_Hoover : public MD_base -{ +class Nose_Hoover : public MD_base { public: Nose_Hoover(const MD_para& mdp_in, const bool cal_stress_in, @@ -21,37 +20,22 @@ class Nose_Hoover : public MD_base /// @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); diff --git a/source/source_md/run_md.cpp b/source/source_md/run_md.cpp index c176aa8212e..02af0b8675c 100644 --- a/source/source_md/run_md.cpp +++ b/source/source_md/run_md.cpp @@ -29,15 +29,12 @@ #include #endif -namespace Run_MD -{ +namespace Run_MD { -void prepare_mdcell(MDCell& mdcell, const Parameter& param_in, DomainDecomposition& decomp) -{ +void prepare_mdcell(MDCell& mdcell, const Parameter& param_in, DomainDecomposition& decomp) { const Input_para& input = param_in.inp; std::vector effective_replicate = input.cell_replica; - if (input.mdp.md_restart) - { + if (input.mdp.md_restart) { effective_replicate = {1, 1, 1}; } @@ -52,16 +49,14 @@ void prepare_mdcell(MDCell& mdcell, const Parameter& param_in, DomainDecompositi 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); const std::vector owned_atoms = decomp.split_owned_atoms_from_ucell(ucell); std::vector type_labels; std::vector type_masses; std::vector type_atom_counts; - for (int it = 0; it < ucell.ntype; ++it) - { + for (int it = 0; it < ucell.ntype; ++it) { type_labels.push_back(ucell.atoms[it].label); type_masses.push_back(ucell.atoms[it].mass); type_atom_counts.push_back(ucell.atoms[it].na); @@ -81,56 +76,49 @@ void prepare_mdcell(MDCell& mdcell, UnitCell& ucell, DomainDecomposition& decomp mdcell.mutable_stru_meta() = unitcell::make_stru_meta(ucell); } -namespace -{ +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(); positions.resize(atoms.size()); - for (std::size_t i = 0; i < atoms.size(); ++i) - { + for (std::size_t i = 0; i < atoms.size(); ++i) { positions[i] = atoms[i].cart * lat0; } } /// @brief gather the atomic forces (in Hartree/Bohr) -void collect_plumed_forces(const MDCell& mdcell, std::vector>& forces) -{ +void collect_plumed_forces(const MDCell& mdcell, std::vector>& forces) { const std::vector& atoms = mdcell.owned_atoms(); forces.resize(atoms.size()); - for (std::size_t i = 0; i < atoms.size(); ++i) - { + for (std::size_t i = 0; i < atoms.size(); ++i) { forces[i] = atoms[i].force; } } /// @brief write the (PLUMED-biased) forces back into the cell -void apply_plumed_forces(MDCell& mdcell, const std::vector>& forces) -{ +void apply_plumed_forces(MDCell& mdcell, const std::vector>& forces) { std::vector& atoms = mdcell.owned_atoms(); - for (std::size_t i = 0; i < atoms.size(); ++i) - { + for (std::size_t i = 0; i < atoms.size(); ++i) { atoms[i].force = forces[i]; } } /// @brief 3x3 cell matrix (in Bohr, row-major) expected by PLUMED -void collect_plumed_cell(const MDCell& mdcell, double* cell) -{ +void collect_plumed_cell(const MDCell& mdcell, double* cell) { const double lat0 = mdcell.lat0(); const double* latvec = &mdcell.latvec().e11; // unitless lattice vectors - for (int i = 0; i < 9; ++i) - { + for (int i = 0; i < 9; ++i) { cell[i] = latvec[i] * lat0; } } } // namespace -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& param_in, + DomainDecomposition& decomp) { ModuleBase::TITLE("Run_MD", "md_line"); ModuleBase::timer::start("Run_MD", "md_line"); /// determine the md_type @@ -140,28 +128,18 @@ void md_line(MDCell& mdcell, ModuleESolver::ESolver* p_esolver, const Parameter& const bool cal_stress = param_in.inp.cal_stress; const bool init_vel = param_in.inp.init_vel; const int my_rank = param_in.globalv.myrank; - if (param_in.mdp.md_type == "fire") - { + if (param_in.mdp.md_type == "fire") { mdrun = new FIRE(mdp, cal_stress, init_vel, my_rank, param_in.inp.force_thr, mdcell); - } - else if ((param_in.mdp.md_type == "nvt" && param_in.mdp.md_thermostat == "nhc") || param_in.mdp.md_type == "npt") - { + } else if ((param_in.mdp.md_type == "nvt" && param_in.mdp.md_thermostat == "nhc") + || param_in.mdp.md_type == "npt") { mdrun = new Nose_Hoover(mdp, cal_stress, init_vel, my_rank, mdcell); - } - else if (param_in.mdp.md_type == "nve" || param_in.mdp.md_type == "nvt") - { + } else if (param_in.mdp.md_type == "nve" || param_in.mdp.md_type == "nvt") { mdrun = new Verlet(mdp, cal_stress, init_vel, my_rank, mdcell); - } - else if (param_in.mdp.md_type == "langevin") - { + } else if (param_in.mdp.md_type == "langevin") { mdrun = new Langevin(mdp, cal_stress, init_vel, my_rank, mdcell); - } - else if (param_in.mdp.md_type == "msst") - { + } else if (param_in.mdp.md_type == "msst") { mdrun = new MSST(mdp, cal_stress, init_vel, my_rank, mdcell); - } - else - { + } else { ModuleBase::WARNING_QUIT("md_line", "no such md_type!"); } @@ -186,36 +164,30 @@ void md_line(MDCell& mdcell, ModuleESolver::ESolver* p_esolver, const Parameter& }; /// md cycle, mohan update 2026-01-04, change '<=' to '<' - while ((mdrun->step_ + mdrun->step_rst_) < param_in.mdp.md_nstep && !mdrun->stop) - { + while ((mdrun->step_ + mdrun->step_rst_) < param_in.mdp.md_nstep && !mdrun->stop) { #ifdef __JSON // JSON output currently follows the UnitCell-backed electronic-structure path. // Start one output record before the solver appends SCF information for this MD step. - if (mdcell.has_backing_unitcell()) - { + if (mdcell.has_backing_unitcell()) { Json::init_output_array_obj(); } #endif - if (mdrun->step_ == 0) - { + if (mdrun->step_ == 0) { mdrun->setup(p_esolver, param_in.globalv.global_readin_dir, decomp); - if (param_in.mdp.plumed) - { + if (param_in.mdp.plumed) { int nproc = 1; #ifdef __MPI MPI_Comm_size(MPI_COMM_WORLD, &nproc); #endif - if (nproc > 1) - { + 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)"); } std::vector plumed_masses(mdcell.owned_atoms().size()); - for (std::size_t i = 0; i < mdcell.owned_atoms().size(); ++i) - { + 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, @@ -226,9 +198,7 @@ void md_line(MDCell& mdcell, ModuleESolver::ESolver* p_esolver, const Parameter& /// forces, so PLUMED can be called right away run_plumed_step(mdrun->step_ + mdrun->step_rst_); } - } - else - { + } else { // mohan add 2026-01-04 const int stress_step = 0; const int force_step = 0; @@ -246,8 +216,7 @@ void md_line(MDCell& mdcell, ModuleESolver::ESolver* p_esolver, const Parameter& mdrun->virial, param_in.mdp.md_out_force); - if (plumed.active()) - { + if (plumed.active()) { /// let PLUMED compute the collective variables and add the /// biasing forces; they are used by the second half of the /// current step and by the first half of the next one @@ -261,8 +230,7 @@ void md_line(MDCell& mdcell, ModuleESolver::ESolver* p_esolver, const Parameter& } 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, @@ -271,10 +239,8 @@ void md_line(MDCell& mdcell, ModuleESolver::ESolver* p_esolver, const Parameter& mdrun->virial); } - if (param_in.mdp.md_restartfreq > 0 && (mdrun->step_ + mdrun->step_rst_) % param_in.mdp.md_restartfreq == 0) - { - if (mdcell.has_backing_unitcell()) - { + if (param_in.mdp.md_restartfreq > 0 && (mdrun->step_ + mdrun->step_rst_) % param_in.mdp.md_restartfreq == 0) { + if (mdcell.has_backing_unitcell()) { mdcell.sync_backing_unitcell(); } std::stringstream file; diff --git a/source/source_md/run_md.h b/source/source_md/run_md.h index 2beaf7c93a4..29f689840bc 100644 --- a/source/source_md/run_md.h +++ b/source/source_md/run_md.h @@ -6,8 +6,7 @@ class DomainDecomposition; class UnitCell; struct Parameter; -namespace ModuleESolver -{ +namespace ModuleESolver { class ESolver; } @@ -15,8 +14,7 @@ class ESolver; * @brief the md loop line * */ -namespace Run_MD -{ +namespace Run_MD { void prepare_mdcell(MDCell& mdcell, const Parameter& param_in, DomainDecomposition& decomp); void prepare_mdcell(MDCell& mdcell, UnitCell& ucell, DomainDecomposition& decomp); diff --git a/source/source_md/test/fire_test.cpp b/source/source_md/test/fire_test.cpp index 05ae519c433..541c8e3a7fc 100644 --- a/source/source_md/test/fire_test.cpp +++ b/source/source_md/test/fire_test.cpp @@ -33,8 +33,7 @@ * - output MD information such as energy, temperature, and pressure */ -class FIREtest : public testing::Test -{ +class FIREtest : public testing::Test { protected: MD_base* mdrun; UnitCell ucell; @@ -46,8 +45,7 @@ class FIREtest : public testing::Test const std::string out_dir = "./"; const std::string readin_dir = "./"; - void SetUp() - { + void SetUp() { Setcell::setupcell(ucell); Setcell::parameters(inp); @@ -59,15 +57,13 @@ class FIREtest : public testing::Test mdrun->setup(p_esolver, readin_dir, decomp); } - void TearDown() - { + void TearDown() { delete mdrun; delete p_esolver; } }; -TEST_F(FIREtest, Setup) -{ +TEST_F(FIREtest, Setup) { EXPECT_NEAR(mdrun->t_current * ModuleBase::Hartree_to_K, 299.99999999999994, doublethreshold); EXPECT_NEAR(mdrun->stress(0, 0), 6.0100555286436806e-06, doublethreshold); EXPECT_NEAR(mdrun->stress(0, 1), -1.4746713013791574e-06, doublethreshold); @@ -80,8 +76,7 @@ TEST_F(FIREtest, Setup) EXPECT_NEAR(mdrun->stress(2, 2), 1.6060561926126463e-06, doublethreshold); } -TEST_F(FIREtest, FirstHalf) -{ +TEST_F(FIREtest, FirstHalf) { mdrun->first_half(GlobalV::ofs_running); EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, @@ -135,8 +130,7 @@ TEST_F(FIREtest, FirstHalf) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, -2.5498181033639999e-05, doublethreshold); } -TEST_F(FIREtest, SecondHalf) -{ +TEST_F(FIREtest, SecondHalf) { mdrun->first_half(GlobalV::ofs_running); mdrun->second_half(); @@ -191,8 +185,7 @@ TEST_F(FIREtest, SecondHalf) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, -2.5498181033639999e-05, doublethreshold); } -TEST_F(FIREtest, WriteRestart) -{ +TEST_F(FIREtest, WriteRestart) { mdrun->step_ = 1; mdrun->step_rst_ = 2; mdrun->write_restart(out_dir); @@ -214,8 +207,7 @@ TEST_F(FIREtest, WriteRestart) ifs.close(); } -TEST_F(FIREtest, Restart) -{ +TEST_F(FIREtest, Restart) { mdrun->restart(readin_dir); remove("Restart_md.txt"); @@ -227,8 +219,7 @@ TEST_F(FIREtest, Restart) EXPECT_EQ(fire->get_md_dt(), 41.3414); } -TEST_F(FIREtest, PrintMD) -{ +TEST_F(FIREtest, PrintMD) { std::ofstream ofs("running_fire.log"); mdrun->print_md(ofs, true); ofs.close(); diff --git a/source/source_md/test/langevin_test.cpp b/source/source_md/test/langevin_test.cpp index cbb62d1a355..13aabb5b625 100644 --- a/source/source_md/test/langevin_test.cpp +++ b/source/source_md/test/langevin_test.cpp @@ -33,8 +33,7 @@ * - output MD information such as energy, temperature, and pressure */ -class Langevin_test : public testing::Test -{ +class Langevin_test : public testing::Test { protected: MD_base* mdrun; UnitCell ucell; @@ -46,8 +45,7 @@ class Langevin_test : public testing::Test const std::string out_dir = "./"; const std::string readin_dir = "./"; - void SetUp() - { + void SetUp() { Setcell::setupcell(ucell); Setcell::parameters(inp); @@ -59,15 +57,13 @@ class Langevin_test : public testing::Test mdrun->setup(p_esolver, readin_dir, decomp); } - void TearDown() - { + void TearDown() { delete mdrun; delete p_esolver; } }; -TEST_F(Langevin_test, setup) -{ +TEST_F(Langevin_test, setup) { EXPECT_NEAR(mdrun->t_current * ModuleBase::Hartree_to_K, 299.99999999999994, doublethreshold); EXPECT_NEAR(mdrun->stress(0, 0), 6.0100555286436806e-06, doublethreshold); EXPECT_NEAR(mdrun->stress(0, 1), -1.4746713013791574e-06, doublethreshold); @@ -80,8 +76,7 @@ TEST_F(Langevin_test, setup) EXPECT_NEAR(mdrun->stress(2, 2), 1.6060561926126463e-06, doublethreshold); } -TEST_F(Langevin_test, first_half) -{ +TEST_F(Langevin_test, first_half) { mdrun->first_half(GlobalV::ofs_running); EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, @@ -135,8 +130,7 @@ TEST_F(Langevin_test, first_half) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, 2.189013211254232e-05, doublethreshold); } -TEST_F(Langevin_test, second_half) -{ +TEST_F(Langevin_test, second_half) { mdrun->first_half(GlobalV::ofs_running); mdrun->second_half(); ; @@ -192,8 +186,7 @@ TEST_F(Langevin_test, second_half) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, 9.6048301397465443e-05, doublethreshold); } -TEST_F(Langevin_test, write_restart) -{ +TEST_F(Langevin_test, write_restart) { mdrun->step_ = 1; mdrun->step_rst_ = 2; mdrun->write_restart(out_dir); @@ -205,16 +198,14 @@ TEST_F(Langevin_test, write_restart) ifs.close(); } -TEST_F(Langevin_test, restart) -{ +TEST_F(Langevin_test, restart) { mdrun->restart(readin_dir); remove("Restart_md.txt"); EXPECT_EQ(mdrun->step_rst_, 3); } -TEST_F(Langevin_test, print_md) -{ +TEST_F(Langevin_test, print_md) { std::ofstream ofs("running_langevin.log"); mdrun->print_md(ofs, true); ofs.close(); diff --git a/source/source_md/test/lj_pot_test.cpp b/source/source_md/test/lj_pot_test.cpp index 05d619658cf..33eaa8e4b54 100644 --- a/source/source_md/test/lj_pot_test.cpp +++ b/source/source_md/test/lj_pot_test.cpp @@ -17,8 +17,7 @@ * - calculate energy, force, virial for lj pot */ -class LJ_pot_test : public testing::Test -{ +class LJ_pot_test : public testing::Test { protected: ModuleBase::matrix stress; double potential; @@ -26,8 +25,7 @@ class LJ_pot_test : public testing::Test UnitCell ucell; Input_para inp; - void SetUp() - { + void SetUp() { Setcell::setupcell(ucell); natom = ucell.nat; @@ -36,13 +34,10 @@ class LJ_pot_test : public testing::Test Setcell::parameters(inp); } - void TearDown() - { - } + void TearDown() {} }; -TEST_F(LJ_pot_test, potential) -{ +TEST_F(LJ_pot_test, potential) { ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); MDCell mdcell; mdcell = Setcell::setup_mdcell(ucell); @@ -53,8 +48,7 @@ TEST_F(LJ_pot_test, potential) EXPECT_NEAR(potential, -0.011957818623534381, doublethreshold); } -TEST_F(LJ_pot_test, unitcell_compatibility) -{ +TEST_F(LJ_pot_test, unitcell_compatibility) { ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); p_esolver->before_all_runners(ucell, inp); p_esolver->runner(ucell, 0); @@ -66,8 +60,7 @@ TEST_F(LJ_pot_test, unitcell_compatibility) delete p_esolver; } -TEST_F(LJ_pot_test, force) -{ +TEST_F(LJ_pot_test, force) { ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); MDCell mdcell; mdcell = Setcell::setup_mdcell(ucell); @@ -90,8 +83,7 @@ TEST_F(LJ_pot_test, force) EXPECT_NEAR(atoms[3].force.z, 0.0, doublethreshold); } -TEST_F(LJ_pot_test, mdcell_cal_force) -{ +TEST_F(LJ_pot_test, mdcell_cal_force) { ModuleESolver::ESolver_LJ p_esolver; MDCell mdcell; mdcell = Setcell::setup_mdcell(ucell); @@ -103,8 +95,7 @@ TEST_F(LJ_pot_test, mdcell_cal_force) ModuleBase::matrix force; p_esolver.cal_force(mdcell, force); - for (int iat = 0; iat < mdcell.owned_atoms().size(); ++iat) - { + for (int iat = 0; iat < mdcell.owned_atoms().size(); ++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); @@ -112,8 +103,7 @@ TEST_F(LJ_pot_test, mdcell_cal_force) } } -TEST_F(LJ_pot_test, stress) -{ +TEST_F(LJ_pot_test, stress) { ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); MDCell mdcell; mdcell = Setcell::setup_mdcell(ucell); @@ -132,8 +122,7 @@ TEST_F(LJ_pot_test, stress) EXPECT_NEAR(stress(2, 2), 6.4275429572682057e-07, doublethreshold); } -TEST_F(LJ_pot_test, mdcell_stress_includes_external_pressure) -{ +TEST_F(LJ_pot_test, mdcell_stress_includes_external_pressure) { ModuleESolver::ESolver_LJ p_esolver; MDCell mdcell; mdcell = Setcell::setup_mdcell(ucell); @@ -163,17 +152,14 @@ TEST_F(LJ_pot_test, mdcell_stress_includes_external_pressure) input.press3 = saved_press3; } -TEST_F(LJ_pot_test, RcutSearchRadius) -{ +TEST_F(LJ_pot_test, RcutSearchRadius) { 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); - for (int i = 0; i < ucell.ntype; i++) - { - for (int j = 0; j < ucell.ntype; j++) - { + for (int i = 0; i < ucell.ntype; i++) { + for (int j = 0; j < ucell.ntype; j++) { EXPECT_NEAR(p_esolver->get_lj_rcut()(i, j), 3.0 * ModuleBase::ANGSTROM_AU, doublethreshold); } } @@ -188,8 +174,7 @@ TEST_F(LJ_pot_test, RcutSearchRadius) EXPECT_NEAR(p_esolver->get_search_radius(), 5.0 * ModuleBase::ANGSTROM_AU + 0.01, doublethreshold); } -TEST_F(LJ_pot_test, SetC6C12) -{ +TEST_F(LJ_pot_test, SetC6C12) { ModuleESolver::ESolver_LJ* p_esolver = new ModuleESolver::ESolver_LJ(); ucell.ntype = 2; @@ -200,10 +185,8 @@ TEST_F(LJ_pot_test, SetC6C12) 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++) - { + for (int i = 0; i < ucell.ntype; i++) { + for (int j = 0; j <= i; j++) { 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), @@ -222,14 +205,12 @@ TEST_F(LJ_pot_test, SetC6C12) p_esolver->set_c6_c12_for_testing(ucell.ntype, rule, lj_epsilon, lj_sigma); - for (int i = 0; i < ucell.ntype; i++) - { + for (int i = 0; i < ucell.ntype; i++) { double temp = pow(lj_sigma[i] * ModuleBase::ANGSTROM_AU, 6); EXPECT_NEAR(p_esolver->get_lj_c6()(i, i), 4.0 * lj_epsilon[i] * temp / ModuleBase::Ry_to_eV, doublethreshold); EXPECT_NEAR(p_esolver->get_lj_c12()(i, i), p_esolver->get_lj_c6()(i, i) * temp, doublethreshold); - for (int j = 0; j < i; j++) - { + 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); @@ -248,10 +229,8 @@ TEST_F(LJ_pot_test, SetC6C12) 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++) - { + 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, @@ -263,8 +242,7 @@ TEST_F(LJ_pot_test, SetC6C12) } } -TEST_F(LJ_pot_test, CalEnShift) -{ +TEST_F(LJ_pot_test, CalEnShift) { ModuleESolver::ESolver_LJ* p_esolver = new ModuleESolver::ESolver_LJ(); ucell.ntype = 2; @@ -278,10 +256,8 @@ TEST_F(LJ_pot_test, CalEnShift) // false p_esolver->cal_en_shift_for_testing(ucell.ntype, false); - for (int i = 0; i < ucell.ntype; i++) - { - for (int j = 0; j < ucell.ntype; j++) - { + for (int i = 0; i < ucell.ntype; i++) { + for (int j = 0; j < ucell.ntype; j++) { EXPECT_DOUBLE_EQ(p_esolver->get_en_shift()(i, j), 0.0); } } diff --git a/source/source_md/test/md_func_test.cpp b/source/source_md/test/md_func_test.cpp index 0ae45e9b06b..8edc3c1a4ee 100644 --- a/source/source_md/test/md_func_test.cpp +++ b/source/source_md/test/md_func_test.cpp @@ -8,18 +8,15 @@ #define doublethreshold 1e-12 -namespace -{ -void initialize_mdcell_from_ucell(MDCell& mdcell, UnitCell& ucell, DomainDecomposition& decomp) -{ +namespace { +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); const std::vector owned_atoms = decomp.split_owned_atoms_from_ucell(ucell); std::vector type_labels; std::vector type_masses; std::vector type_atom_counts; - for (int it = 0; it < ucell.ntype; ++it) - { + for (int it = 0; it < ucell.ntype; ++it) { type_labels.push_back(ucell.atoms[it].label); type_masses.push_back(ucell.atoms[it].mass); type_atom_counts.push_back(ucell.atoms[it].na); @@ -76,19 +73,15 @@ void initialize_mdcell_from_ucell(MDCell& mdcell, UnitCell& ucell, DomainDecompo * - test the current_md_info function with an incorrect file path */ -class MD_func_test : public MdFuncTestFixture -{ -}; +class MD_func_test : public MdFuncTestFixture {}; -TEST_F(MD_func_test, gaussrand) -{ +TEST_F(MD_func_test, gaussrand) { EXPECT_DOUBLE_EQ(MD_func::gaussrand(), 1.1122716058967226); EXPECT_DOUBLE_EQ(MD_func::gaussrand(), -0.34532367182326629); EXPECT_DOUBLE_EQ(MD_func::gaussrand(), 0.60805637857480721); } -TEST_F(MD_func_test, readvel) -{ +TEST_F(MD_func_test, readvel) { MD_func::read_vel(ucell, vel); EXPECT_DOUBLE_EQ(vel[0].x, -0.0001320807363640); @@ -105,11 +98,9 @@ TEST_F(MD_func_test, readvel) EXPECT_DOUBLE_EQ(vel[3].z, -2.83313122596e-05); } -TEST_F(MD_func_test, RescaleVel) -{ +TEST_F(MD_func_test, RescaleVel) { int frozen_freedom = 3; - for (int i = 0; i < natom; ++i) - { + for (int i = 0; i < natom; ++i) { allmass[i] = 39.948 / ModuleBase::AU_to_MASS; vel[i].x = 0.1; vel[i].y = 0.2; @@ -132,14 +123,12 @@ TEST_F(MD_func_test, RescaleVel) EXPECT_DOUBLE_EQ(vel[3].z, 0.00013737032325207373); } -TEST_F(MD_func_test, compute_stress) -{ +TEST_F(MD_func_test, compute_stress) { const ModuleBase::Vector3 test_velocity(0.1, 0.2, 0.3); 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); @@ -154,13 +143,11 @@ TEST_F(MD_func_test, compute_stress) EXPECT_DOUBLE_EQ(stress(2, 2), 26.215470600896506); } -TEST_F(MD_func_test, dump_info) -{ +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); @@ -298,8 +285,7 @@ TEST_F(MD_func_test, dump_info) remove("MD_dump"); } -TEST_F(MD_func_test, print_stress) -{ +TEST_F(MD_func_test, print_stress) { GlobalV::ofs_running.open("running.log"); MD_func::print_stress(GlobalV::ofs_running, virial, stress); @@ -337,8 +323,7 @@ TEST_F(MD_func_test, print_stress) remove("running.log"); } -TEST_F(MD_func_test, current_md_info_mdcell_accepts_step_only_restart) -{ +TEST_F(MD_func_test, current_md_info_mdcell_accepts_step_only_restart) { std::ofstream file("Restart_md.txt"); file << 123; file.close(); @@ -355,22 +340,19 @@ TEST_F(MD_func_test, current_md_info_mdcell_accepts_step_only_restart) remove("Restart_md.txt"); } -TEST_F(MD_func_test, global_dof_mdcell) -{ +TEST_F(MD_func_test, global_dof_mdcell) { MDCell mdcell; DomainDecomposition decomp; 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; } EXPECT_EQ(MD_func::global_dof(mdcell), 6); } -TEST_F(MD_func_test, current_step_warning) -{ +TEST_F(MD_func_test, current_step_warning) { MDCell mdcell; DomainDecomposition decomp; initialize_mdcell_from_ucell(mdcell, ucell, decomp); diff --git a/source/source_md/test/md_test_fixture.h b/source/source_md/test/md_test_fixture.h index dac429f978a..b406ed28d3b 100644 --- a/source/source_md/test/md_test_fixture.h +++ b/source/source_md/test/md_test_fixture.h @@ -12,15 +12,13 @@ #include #include -class MdTestBase : public testing::Test -{ +class MdTestBase : public testing::Test { protected: UnitCell ucell; Input_para inp; std::unique_ptr p_esolver; - void SetUp() override - { + void SetUp() override { Setcell::setupcell(ucell); Setcell::parameters(inp); @@ -30,8 +28,7 @@ class MdTestBase : public testing::Test }; template -class MdIntegratorFixture : public MdTestBase -{ +class MdIntegratorFixture : public MdTestBase { protected: MDCell mdcell; DomainDecomposition decomp; @@ -39,8 +36,7 @@ class MdIntegratorFixture : public MdTestBase /// setup() takes the directory explicitly, so the fixture owns it const std::string readin_dir = "./"; - void SetUp() override - { + void SetUp() override { MdTestBase::SetUp(); mdcell = Setcell::setup_mdcell(ucell); decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0); @@ -50,8 +46,7 @@ class MdIntegratorFixture : public MdTestBase } }; -class MdFuncTestFixture : public testing::Test -{ +class MdFuncTestFixture : public testing::Test { protected: UnitCell ucell; std::vector allmass_store; @@ -74,8 +69,7 @@ class MdFuncTestFixture : public testing::Test /// the md writers take the directory explicitly, so the fixture owns it const std::string out_dir = "./"; - void SetUp() override - { + void SetUp() override { Setcell::setupcell(ucell); Setcell::parameters(inp); natom = ucell.nat; @@ -95,8 +89,7 @@ class MdFuncTestFixture : public testing::Test } }; -class LjPotTestFixture : public testing::Test -{ +class LjPotTestFixture : public testing::Test { protected: std::vector> force_store; ModuleBase::Vector3* force = nullptr; @@ -106,8 +99,7 @@ class LjPotTestFixture : public testing::Test UnitCell ucell; Input_para input; - void SetUp() override - { + void SetUp() override { Setcell::setupcell(ucell); natom = ucell.nat; diff --git a/source/source_md/test/msst_test.cpp b/source/source_md/test/msst_test.cpp index b037304de20..da08302856f 100644 --- a/source/source_md/test/msst_test.cpp +++ b/source/source_md/test/msst_test.cpp @@ -33,8 +33,7 @@ * - output MD information such as energy, temperature, and pressure */ -class MSST_test : public testing::Test -{ +class MSST_test : public testing::Test { protected: MD_base* mdrun; UnitCell ucell; @@ -46,8 +45,7 @@ class MSST_test : public testing::Test const std::string out_dir = "./"; const std::string readin_dir = "./"; - void SetUp() - { + void SetUp() { Setcell::setupcell(ucell); Setcell::parameters(inp); @@ -59,15 +57,13 @@ class MSST_test : public testing::Test mdrun->setup(p_esolver, readin_dir, decomp); } - void TearDown() - { + void TearDown() { delete mdrun; delete p_esolver; } }; -TEST_F(MSST_test, setup) -{ +TEST_F(MSST_test, setup) { EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.x, -0.0001314186733659715, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.y, 7.0985331994796372e-05, doublethreshold); EXPECT_NEAR(mdcell.owned_atoms()[static_cast(0)].vel.z, -1.3947731701005279e-05, doublethreshold); @@ -92,8 +88,7 @@ TEST_F(MSST_test, setup) EXPECT_NEAR(mdrun->stress(2, 2), 1.5964231736437884e-06, doublethreshold); } -TEST_F(MSST_test, first_half) -{ +TEST_F(MSST_test, first_half) { mdrun->first_half(GlobalV::ofs_running); EXPECT_NEAR(mdcell.lat0(), 1.0, doublethreshold); @@ -160,8 +155,7 @@ TEST_F(MSST_test, first_half) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, -2.8189297471809918e-05, doublethreshold); } -TEST_F(MSST_test, second_half) -{ +TEST_F(MSST_test, second_half) { mdrun->first_half(GlobalV::ofs_running); mdrun->second_half(); ; @@ -230,8 +224,7 @@ TEST_F(MSST_test, second_half) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, -2.8200698165338931e-05, doublethreshold); } -TEST_F(MSST_test, write_restart) -{ +TEST_F(MSST_test, write_restart) { mdrun->step_ = 1; mdrun->step_rst_ = 2; mdrun->write_restart(out_dir); @@ -255,8 +248,7 @@ TEST_F(MSST_test, write_restart) ifs.close(); } -TEST_F(MSST_test, restart) -{ +TEST_F(MSST_test, restart) { mdrun->restart(readin_dir); remove("Restart_md.txt"); @@ -269,8 +261,7 @@ TEST_F(MSST_test, restart) EXPECT_EQ(msst->get_lag_pos(), 0); } -TEST_F(MSST_test, print_md) -{ +TEST_F(MSST_test, print_md) { std::ofstream ofs("running_msst.log"); mdrun->print_md(ofs, true); ofs.close(); diff --git a/source/source_md/test/nhchain_test.cpp b/source/source_md/test/nhchain_test.cpp index e0756fe9675..9d43f5c5d81 100644 --- a/source/source_md/test/nhchain_test.cpp +++ b/source/source_md/test/nhchain_test.cpp @@ -31,8 +31,7 @@ * - Nose_Hoover::print_md * - output MD information such as energy, temperature, and pressure */ -class NHC_test : public testing::Test -{ +class NHC_test : public testing::Test { protected: MD_base* mdrun; UnitCell ucell; @@ -44,8 +43,7 @@ class NHC_test : public testing::Test const std::string out_dir = "./"; const std::string readin_dir = "./"; - void SetUp() - { + void SetUp() { Setcell::setupcell(ucell); Setcell::parameters(inp); @@ -61,15 +59,13 @@ class NHC_test : public testing::Test mdrun->setup(p_esolver, readin_dir, decomp); } - void TearDown() - { + void TearDown() { delete mdrun; delete p_esolver; } }; -TEST_F(NHC_test, setup) -{ +TEST_F(NHC_test, setup) { EXPECT_NEAR(mdrun->t_current * ModuleBase::Hartree_to_K, 299.99999999999994, doublethreshold); EXPECT_NEAR(mdrun->stress(0, 0), 6.0100555286436806e-06, doublethreshold); EXPECT_NEAR(mdrun->stress(0, 1), -1.4746713013791574e-06, doublethreshold); @@ -82,8 +78,7 @@ TEST_F(NHC_test, setup) EXPECT_NEAR(mdrun->stress(2, 2), 1.6060561926126463e-06, doublethreshold); } -TEST_F(NHC_test, first_half) -{ +TEST_F(NHC_test, first_half) { mdrun->first_half(GlobalV::ofs_running); EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, @@ -137,8 +132,7 @@ TEST_F(NHC_test, first_half) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, -2.9596755163433345e-05, doublethreshold); } -TEST_F(NHC_test, second_half) -{ +TEST_F(NHC_test, second_half) { mdrun->first_half(GlobalV::ofs_running); mdrun->second_half(); @@ -193,8 +187,7 @@ TEST_F(NHC_test, second_half) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, -3.2917171190756263e-05, doublethreshold); } -TEST_F(NHC_test, write_restart) -{ +TEST_F(NHC_test, write_restart) { mdrun->first_half(GlobalV::ofs_running); mdrun->second_half(); @@ -226,8 +219,7 @@ TEST_F(NHC_test, write_restart) ifs.close(); } -TEST_F(NHC_test, restart) -{ +TEST_F(NHC_test, restart) { mdrun->restart(readin_dir); remove("Restart_md.txt"); @@ -259,8 +251,7 @@ TEST_F(NHC_test, restart) EXPECT_EQ(nhc->get_v_peta()[3], 226.197); } -TEST_F(NHC_test, print_md) -{ +TEST_F(NHC_test, print_md) { std::ofstream ofs("running_nhchain.log"); mdrun->print_md(ofs, true); ofs.close(); diff --git a/source/source_md/test/run_md_test.cpp b/source/source_md/test/run_md_test.cpp index 52b4e4f14f0..e4a22da7fb2 100644 --- a/source/source_md/test/run_md_test.cpp +++ b/source/source_md/test/run_md_test.cpp @@ -6,8 +6,7 @@ #include "gtest/gtest.h" -TEST(RunMDTest, prepare_mdcell_from_unitcell) -{ +TEST(RunMDTest, prepare_mdcell_from_unitcell) { UnitCell ucell; ucell.ntype = 1; ucell.nat = 1; diff --git a/source/source_md/test/setcell.h b/source/source_md/test/setcell.h index b6027cdb293..2a39f8ffa27 100644 --- a/source/source_md/test/setcell.h +++ b/source/source_md/test/setcell.h @@ -17,20 +17,15 @@ #include #include -Magnetism::Magnetism() -{ +Magnetism::Magnetism() { this->tot_mag = 0.0; this->abs_mag = 0.0; } -Magnetism::~Magnetism() -{ -} +Magnetism::~Magnetism() {} -class Setcell -{ +class Setcell { public: - static void setupcell(UnitCell& ucell) - { + static void setupcell(UnitCell& ucell) { ucell.ntype = 1; ucell.atoms = new Atom[ucell.ntype]; @@ -82,8 +77,7 @@ class Setcell ucell.atoms[0].vel[1].set(0.000153039878532, -0.000146533266608, 9.64491480698e-05); ucell.atoms[0].vel[2].set(-0.000133789480226, -3.0451038112e-06, -5.40998380137e-05); ucell.atoms[0].vel[3].set(0.000112830338059, 7.82354274358e-05, -2.83313122596e-05); - for (int ia = 0; ia < 4; ++ia) - { + for (int ia = 0; ia < 4; ++ia) { ucell.atoms[0].tau[ia] = ucell.atoms[0].taud[ia] * ucell.latvec; ucell.atoms[0].mbl[ia].set(1, 1, 1); } @@ -104,8 +98,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,19 +133,16 @@ 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; std::vector type_masses; std::vector type_atom_counts; - for (int it = 0; it < ucell.ntype; ++it) - { + for (int it = 0; it < ucell.ntype; ++it) { type_labels.push_back(ucell.atoms[it].label); type_masses.push_back(ucell.atoms[it].mass); type_atom_counts.push_back(ucell.atoms[it].na); - for (int ia = 0; ia < ucell.atoms[it].na; ++ia) - { + 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], @@ -180,8 +170,7 @@ class Setcell 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), ModuleBase::Vector3(0.52, 0.52, 0.0), ModuleBase::Vector3(0.51, 0.0, 0.5), diff --git a/source/source_md/test/test_md_func.cpp b/source/source_md/test/test_md_func.cpp index c4d36889641..789ee5e3302 100644 --- a/source/source_md/test/test_md_func.cpp +++ b/source/source_md/test/test_md_func.cpp @@ -6,49 +6,34 @@ #include -namespace -{ -class GhostForceSolver : public ModuleESolver::ESolver -{ +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 - { + 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); 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 - { + double cal_energy() override { return 6.0; } + 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; } }; -TEST(MdForceVirialTest, ReturnsGhostForcesBeforeConvertingUnitsExactlyOnce) -{ +TEST(MdForceVirialTest, ReturnsGhostForcesBeforeConvertingUnitsExactlyOnce) { MDCell cell; DomainDecomposition decomp; const ModuleBase::CommunicationDomain domain = ModuleBase::world_comm_domain(); @@ -75,15 +60,13 @@ TEST(MdForceVirialTest, ReturnsGhostForcesBeforeConvertingUnitsExactlyOnce) GhostForceSolver solver; double potential = 0.0; ModuleBase::matrix virial(3, 3); - for (int step = 0; step < 2; ++step) - { + for (int step = 0; step < 2; ++step) { 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 bb285efb4c6..82a978172aa 100644 --- a/source/source_md/test/verlet_test.cpp +++ b/source/source_md/test/verlet_test.cpp @@ -34,8 +34,7 @@ * - output MD information such as energy, temperature, and pressure */ -class Verlet_test : public testing::Test -{ +class Verlet_test : public testing::Test { protected: MD_base* mdrun; UnitCell ucell; @@ -47,8 +46,7 @@ class Verlet_test : public testing::Test const std::string out_dir = "./"; const std::string readin_dir = "./"; - void SetUp() - { + void SetUp() { Setcell::setupcell(ucell); Setcell::parameters(inp); @@ -60,15 +58,13 @@ class Verlet_test : public testing::Test mdrun->setup(p_esolver, readin_dir, decomp); } - void TearDown() - { + void TearDown() { delete mdrun; delete p_esolver; } }; -TEST_F(Verlet_test, setup) -{ +TEST_F(Verlet_test, setup) { EXPECT_NEAR(mdrun->t_current * ModuleBase::Hartree_to_K, 299.99999999999994, doublethreshold); EXPECT_NEAR(mdrun->stress(0, 0), 6.0100555286436806e-06, doublethreshold); EXPECT_NEAR(mdrun->stress(0, 1), -1.4746713013791574e-06, doublethreshold); @@ -81,8 +77,7 @@ TEST_F(Verlet_test, setup) EXPECT_NEAR(mdrun->stress(2, 2), 1.6060561926126463e-06, doublethreshold); } -TEST_F(Verlet_test, first_half) -{ +TEST_F(Verlet_test, first_half) { mdrun->first_half(GlobalV::ofs_running); EXPECT_NEAR(Setcell::fractional_displacement(mdcell.owned_atoms()[static_cast(0)]).x, @@ -136,8 +131,7 @@ TEST_F(Verlet_test, first_half) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, -2.83313122596e-05, doublethreshold); } -TEST_F(Verlet_test, NVE) -{ +TEST_F(Verlet_test, NVE) { mdrun->first_half(GlobalV::ofs_running); inp.mdp.md_type = "nve"; mdrun->second_half(); @@ -194,8 +188,7 @@ TEST_F(Verlet_test, NVE) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, -2.83313122596e-05, doublethreshold); } -TEST_F(Verlet_test, Anderson) -{ +TEST_F(Verlet_test, Anderson) { mdrun->first_half(GlobalV::ofs_running); inp.mdp.md_type = "nvt"; inp.mdp.md_thermostat = "anderson"; @@ -253,8 +246,7 @@ TEST_F(Verlet_test, Anderson) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, -2.83313122596e-05, doublethreshold); } -TEST_F(Verlet_test, Berendsen) -{ +TEST_F(Verlet_test, Berendsen) { mdrun->first_half(GlobalV::ofs_running); inp.mdp.md_type = "nvt"; inp.mdp.md_thermostat = "berendsen"; @@ -312,8 +304,7 @@ TEST_F(Verlet_test, Berendsen) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, -2.8329987777389342e-05, doublethreshold); } -TEST_F(Verlet_test, rescaling) -{ +TEST_F(Verlet_test, rescaling) { mdrun->first_half(GlobalV::ofs_running); inp.mdp.md_type = "nvt"; inp.mdp.md_thermostat = "rescaling"; @@ -371,8 +362,7 @@ TEST_F(Verlet_test, rescaling) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, -2.8328663233253657e-05, doublethreshold); } -TEST_F(Verlet_test, rescale_v) -{ +TEST_F(Verlet_test, rescale_v) { mdrun->first_half(GlobalV::ofs_running); inp.mdp.md_type = "nvt"; inp.mdp.md_thermostat = "rescale_v"; @@ -430,8 +420,7 @@ TEST_F(Verlet_test, rescale_v) EXPECT_NEAR(mdcell.owned_atoms()[static_cast(3)].vel.z, -2.8328663233253657e-05, doublethreshold); } -TEST_F(Verlet_test, CSVR) -{ +TEST_F(Verlet_test, CSVR) { std::ofstream ofs; mdrun->first_half(ofs); inp.mdp.md_type = "nvt"; @@ -457,8 +446,7 @@ TEST_F(Verlet_test, CSVR) EXPECT_LT(temp, 1000.0); } -TEST_F(Verlet_test, write_restart) -{ +TEST_F(Verlet_test, write_restart) { mdrun->step_ = 1; mdrun->step_rst_ = 2; mdrun->write_restart(out_dir); @@ -470,16 +458,14 @@ TEST_F(Verlet_test, write_restart) ifs.close(); } -TEST_F(Verlet_test, restart) -{ +TEST_F(Verlet_test, restart) { mdrun->restart(readin_dir); remove("Restart_md.txt"); EXPECT_EQ(mdrun->step_rst_, 3); } -TEST_F(Verlet_test, print_md) -{ +TEST_F(Verlet_test, print_md) { std::ofstream ofs("running_verlet.log"); mdrun->print_md(ofs, true); ofs.close(); diff --git a/source/source_md/verlet.cpp b/source/source_md/verlet.cpp index bf4bc11a768..4584940bca5 100644 --- a/source/source_md/verlet.cpp +++ b/source/source_md/verlet.cpp @@ -12,16 +12,13 @@ Verlet::Verlet(const MD_para& mdp_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) -{ -} + : 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"); @@ -30,8 +27,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"); @@ -41,8 +37,7 @@ void Verlet::first_half(std::ofstream& ofs) ModuleBase::timer::end("Verlet", "first_half"); } -void Verlet::second_half() -{ +void Verlet::second_half() { ModuleBase::TITLE("Verlet", "second_half"); ModuleBase::timer::start("Verlet", "second_half"); @@ -52,72 +47,48 @@ void Verlet::second_half() ModuleBase::timer::end("Verlet", "second_half"); } -void Verlet::apply_thermostat(void) -{ +void Verlet::apply_thermostat(void) { double t_target = 0.0; t_current = MD_func::current_temp(kinetic, mdcell, frozen_freedom_); - if (mdp.md_type == "nve") - { - } - else if (mdp.md_thermostat == "rescaling") - { + if (mdp.md_type == "nve") { + } else if (mdp.md_thermostat == "rescaling") { 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) - { + if (std::abs(t_target - t_current) * ModuleBase::Hartree_to_K > mdp.md_tolerance) { thermalize(0, t_current, t_target); } - } - else if (mdp.md_thermostat == "rescale_v") - { - if ((step_ + step_rst_) % mdp.md_nraise == 0) - { + } else if (mdp.md_thermostat == "rescale_v") { + if ((step_ + step_rst_) % mdp.md_nraise == 0) { 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()) - { - if (static_cast(std::rand()) / RAND_MAX <= 1.0 / mdp.md_nraise) - { + } else if (mdp.md_thermostat == "anderson") { + for (LocalAtom& atom: mdcell.owned_atoms()) { + if (static_cast(std::rand()) / RAND_MAX <= 1.0 / mdp.md_nraise) { const double deviation = sqrt(md_tlast / atom.mass); - for (int k = 0; k < 3; ++k) - { - if (atom.mbl[k]) - { + for (int k = 0; k < 3; ++k) { + if (atom.mbl[k]) { atom.vel[k] = deviation * MD_func::gaussrand(); } } } } - } - else if (mdp.md_thermostat == "berendsen") - { + } else if (mdp.md_thermostat == "berendsen") { 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") - { + } else if (mdp.md_thermostat == "csvr") { t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); apply_csvr(t_current, t_target); - } - else - { + } else { ModuleBase::WARNING_QUIT("Verlet", "No such thermostat!"); } } -void Verlet::thermalize(const int& nraise, const double& current_temp, const double& target_temp) -{ +void Verlet::thermalize(const int& nraise, const double& current_temp, const double& target_temp) { double fac = 0.0; - if (nraise > 0 && current_temp > 0 && target_temp > 0) - { + if (nraise > 0 && current_temp > 0 && target_temp > 0) { fac = sqrt(1 + (target_temp / current_temp - 1) / nraise); - } - else if (nraise == 0 && current_temp > 0 && target_temp > 0) - { + } else if (nraise == 0 && current_temp > 0 && target_temp > 0) { fac = sqrt(target_temp / current_temp); } @@ -125,13 +96,11 @@ void Verlet::thermalize(const int& nraise, const double& current_temp, const dou atom.vel *= fac; } -void Verlet::apply_csvr(const double& current_temp, const double& target_temp) -{ +void Verlet::apply_csvr(const double& current_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) - if (current_temp <= 0.0 || target_temp <= 0.0) - { + if (current_temp <= 0.0 || target_temp <= 0.0) { return; } @@ -147,44 +116,38 @@ void Verlet::apply_csvr(const double& current_temp, const double& target_temp) // Calculate decay factor double factor = 0.0; - if (taut > 0.1) - { + if (taut > 0.1) { factor = exp(-1.0 / taut); } double scale = 1.0; #ifdef __MPI - if (mdcell.mpi_size() > 1) - { - if (mdcell.mpi_rank() == 0) - { + if (mdcell.mpi_size() > 1) { + if (mdcell.mpi_rank() == 0) { const double rr = MD_func::gaussrand(); double sumnoises = 0.0; - for (int i = 0; i < ndeg - 1; ++i) - { + for (int i = 0; i < ndeg - 1; ++i) { const double random_value = MD_func::gaussrand(); 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()); - } - else + } else #endif { const double rr = MD_func::gaussrand(); double sumnoises = 0.0; - for (int i = 0; i < ndeg - 1; ++i) - { + for (int i = 0; i < ndeg - 1; ++i) { const double random_value = MD_func::gaussrand(); 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); } @@ -193,20 +156,17 @@ void Verlet::apply_csvr(const double& current_temp, const double& target_temp) 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 670688251ef..a798cdc89bc 100644 --- a/source/source_md/verlet.h +++ b/source/source_md/verlet.h @@ -7,8 +7,7 @@ * @brief the md methods based on the velocity-Verlet equation * */ -class Verlet : public MD_base -{ +class Verlet : public MD_base { public: Verlet(const MD_para& mdp_in, const bool cal_stress_in,