From 57675b18f0284ff4d3c60e2c1c4c4792abc7a54c Mon Sep 17 00:00:00 2001 From: SY Wang Date: Thu, 15 Oct 2026 22:51:36 +0800 Subject: [PATCH 1/5] Prettify module_vdw --- .pre-commit-config.yaml | 11 +- .../module_vdw/test/vdw_test.cpp | 406 +++++++++--------- .../module_vdw/test/vdwd3_evaluator_test.cpp | 165 +++---- source/source_hamilt/module_vdw/vdw.cpp | 53 ++- source/source_hamilt/module_vdw/vdw.h | 20 +- .../source_hamilt/module_vdw/vdw_parameters.h | 10 +- .../source_hamilt/module_vdw/vdw_xcname.cpp | 7 +- source/source_hamilt/module_vdw/vdwd2.cpp | 4 +- source/source_hamilt/module_vdw/vdwd2.h | 17 +- .../module_vdw/vdwd2_parameters.cpp | 12 +- .../module_vdw/vdwd2_parameters.h | 45 +- source/source_hamilt/module_vdw/vdwd3.cpp | 34 +- source/source_hamilt/module_vdw/vdwd3.h | 7 +- .../source_hamilt/module_vdw/vdwd3_data.cpp | 3 +- .../module_vdw/vdwd3_evaluator.cpp | 190 ++++---- .../module_vdw/vdwd3_parameters.cpp | 25 +- source/source_hamilt/module_vdw/vdwd3_types.h | 4 +- source/source_hamilt/module_vdw/vdwd4.cpp | 24 +- source/source_hamilt/module_vdw/vdwd4.h | 12 +- 19 files changed, 485 insertions(+), 564 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 8873791f2ce..5a76d574d54 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -9,15 +9,12 @@ repos: - repo: https://github.com/pocc/pre-commit-hooks rev: v1.3.5 hooks: - # - id: clang-format - # args: [-i] - - id: clang-tidy - args: [-p=build, --fix-errors] - # - id: oclint - # - id: uncrustify + - id: clang-format + args: [-i, --no-diff] + # - id: clang-tidy + # args: [-p=build, --fix-errors] # - id: cppcheck # - id: cpplint - # - id: include-what-you-use - repo: local hooks: - id: abacus-agent-governance diff --git a/source/source_hamilt/module_vdw/test/vdw_test.cpp b/source/source_hamilt/module_vdw/test/vdw_test.cpp index f03f4fb1588..71e3682ce62 100644 --- a/source/source_hamilt/module_vdw/test/vdw_test.cpp +++ b/source/source_hamilt/module_vdw/test/vdw_test.cpp @@ -1,32 +1,32 @@ -#include "source_cell/unitcell.h" - +#include "mpi.h" #include "source_base/mathzone.h" #include "source_base/vector3.h" -#include"gtest/gtest.h" -#include"gmock/gmock.h" -#include "mpi.h" -#include "source_hamilt/module_vdw/vdwd2_parameters.h" +#include "source_cell/unitcell.h" #include "source_hamilt/module_vdw/vdwd2.h" +#include "source_hamilt/module_vdw/vdwd2_parameters.h" #include "source_hamilt/module_vdw/vdwd3.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" #ifdef __DFTD4 #include "source_hamilt/module_vdw/vdwd4.h" #endif #include "source_hamilt/module_vdw/vdw.h" /************************************************ -* unit test of class VDW and related functions -***********************************************/ + * unit test of class VDW and related functions + ***********************************************/ /** -* - Tested functions: -* - vdw::make_vdw(): -* Based on the value of INPUT.vdw_method, construct -* Vdwd2 or Vdwd3 class, and do the initialization. -* - vdw::Vdw::evaluate(): -* Calculate the requested vdW energy, force and stress in one evaluation. -* - Vdwd2Parameters::initial_parameters() -* - native s-dftd3-compatible D3 adapter -*/ + * - Tested functions: + * - vdw::make_vdw(): + * Based on the value of INPUT.vdw_method, construct + * Vdwd2 or Vdwd3 class, and do the initialization. + * - vdw::Vdw::evaluate(): + * Calculate the requested vdW energy, force and stress in one evaluation. + * - Vdwd2Parameters::initial_parameters() + * - native s-dftd3-compatible D3 adapter + */ pseudo::pseudo() { @@ -59,10 +59,18 @@ Magnetism::~Magnetism() { } -SepPot::SepPot(){} -SepPot::~SepPot(){} -Sep_Cell::Sep_Cell() noexcept {} -Sep_Cell::~Sep_Cell() noexcept {} +SepPot::SepPot() +{ +} +SepPot::~SepPot() +{ +} +Sep_Cell::Sep_Cell() noexcept +{ +} +Sep_Cell::~Sep_Cell() noexcept +{ +} struct atomtype_ { @@ -76,15 +84,21 @@ struct stru_ std::vector all_type; }; -void construct_ucell(stru_ &stru, UnitCell &ucell) +void construct_ucell(stru_& stru, UnitCell& ucell) { std::vector coord = stru.all_type; ucell.a1 = ModuleBase::Vector3(stru.cell[0], stru.cell[1], stru.cell[2]); ucell.a2 = ModuleBase::Vector3(stru.cell[3], stru.cell[4], stru.cell[5]); ucell.a3 = ModuleBase::Vector3(stru.cell[6], stru.cell[7], stru.cell[8]); - ucell.latvec.e11 = stru.cell[0]; ucell.latvec.e12 = stru.cell[1]; ucell.latvec.e13 = stru.cell[2]; - ucell.latvec.e21 = stru.cell[3]; ucell.latvec.e22 = stru.cell[4]; ucell.latvec.e23 = stru.cell[5]; - ucell.latvec.e31 = stru.cell[6]; ucell.latvec.e32 = stru.cell[7]; ucell.latvec.e33 = stru.cell[8]; + ucell.latvec.e11 = stru.cell[0]; + ucell.latvec.e12 = stru.cell[1]; + ucell.latvec.e13 = stru.cell[2]; + ucell.latvec.e21 = stru.cell[3]; + ucell.latvec.e22 = stru.cell[4]; + ucell.latvec.e23 = stru.cell[5]; + ucell.latvec.e31 = stru.cell[6]; + ucell.latvec.e32 = stru.cell[7]; + ucell.latvec.e33 = stru.cell[8]; ucell.lat0 = 10.2; ucell.ntype = stru.all_type.size(); ucell.atoms = new Atom[ucell.ntype]; @@ -126,30 +140,28 @@ void construct_ucell(stru_ &stru, UnitCell &ucell) ucell.omega = std::abs(ucell.latvec.Det()) * ucell.lat0 * ucell.lat0 * ucell.lat0; ucell.itia2iat.create(ucell.ntype, nmax); - int iat=0; - for(int it = 0;it < ucell.ntype;it++) - { - for(int ia=0; ia{0.5, 0.5, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.5}, - std::vector{atomtype_{"Si", - std::vector>{{0., 0., 0.}, - {0.3, 0.25, 0.25}}}, - atomtype_{"Si_empty", - std::vector>{{0.1, 0.1, 0.1}}}}}; + stru_ structure{ + std::vector{0.5, 0.5, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.5}, + std::vector{atomtype_{"Si", std::vector>{{0., 0., 0.}, {0.3, 0.25, 0.25}}}, + atomtype_{"Si_empty", std::vector>{{0.1, 0.1, 0.1}}}}}; construct_ucell(structure, ucell); ucell.atoms[0].ncpp.zv = 4.0; ucell.atoms[1].flag_empty_element = true; @@ -157,9 +169,7 @@ void construct_two_si_with_ghost(UnitCell& ucell) ucell.atoms[1].ncpp.zv = 4.0; } -void expect_empty_atoms_excluded(const vdw::VdwResult& reference, - const vdw::VdwResult& ghost, - int reference_nat) +void expect_empty_atoms_excluded(const vdw::VdwResult& reference, const vdw::VdwResult& ghost, int reference_nat) { EXPECT_NEAR(ghost.energy, reference.energy, 1E-12); EXPECT_EQ(ghost.has_force, reference.has_force); @@ -197,24 +207,22 @@ void expect_empty_atoms_excluded(const vdw::VdwResult& reference, } } -class vdwd2Test: public testing::Test +class vdwd2Test : public testing::Test { - protected: + protected: UnitCell ucell; Input_para input; - void SetUp(){ - stru_ structure{std::vector{0.5,0.5,0.0,0.5,0.0,0.5,0.0,0.5,0.5}, - std::vector{atomtype_{"Si", - std::vector>{ - {0., 0., 0.}, - {0.3, 0.25, 0.25} - }}}}; - construct_ucell(structure,ucell); + void SetUp() + { + stru_ structure{ + std::vector{0.5, 0.5, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.5}, + std::vector{atomtype_{"Si", std::vector>{{0., 0., 0.}, {0.3, 0.25, 0.25}}}}}; + construct_ucell(structure, ucell); input.vdw_method = "d2"; input.vdw_s6 = "0.75"; - input.vdw_d = 20; + input.vdw_d = 20; input.vdw_C6_file = "default"; input.vdw_C6_unit = "Jnm6/mol"; input.vdw_R0_file = "default"; @@ -224,7 +232,8 @@ class vdwd2Test: public testing::Test input.vdw_radius_unit = "Bohr"; } - void TearDown(){ + void TearDown() + { ClearUcell(ucell); } }; @@ -238,9 +247,10 @@ TEST_F(vdwd2Test, D2Default) EXPECT_EQ(vdwd2_test.parameter().damping(), input.vdw_d); EXPECT_EQ(vdwd2_test.parameter().model(), input.vdw_cutoff_type); EXPECT_EQ(vdwd2_test.parameter().radius(), 56.6918); - double Si_C6 = 9.23*1e6 / (ModuleBase::ELECTRONVOLT_SI * ModuleBase::NA) / pow(ModuleBase::BOHR_TO_A, 6)/ ModuleBase::Ry_to_eV; - EXPECT_NEAR(vdwd2_test.parameter().C6().at("Si"), Si_C6,1e-13); - EXPECT_EQ(vdwd2_test.parameter().R0().at("Si"), 1.716/ModuleBase::BOHR_TO_A); + double Si_C6 = 9.23 * 1e6 / (ModuleBase::ELECTRONVOLT_SI * ModuleBase::NA) / pow(ModuleBase::BOHR_TO_A, 6) + / ModuleBase::Ry_to_eV; + EXPECT_NEAR(vdwd2_test.parameter().C6().at("Si"), Si_C6, 1e-13); + EXPECT_EQ(vdwd2_test.parameter().R0().at("Si"), 1.716 / ModuleBase::BOHR_TO_A); EXPECT_EQ(vdwd2_test.parameter().period().x, 2 * ceil(56.6918 / ucell.lat0 / sqrt(ucell.a1.norm2())) + 1); EXPECT_EQ(vdwd2_test.parameter().period().y, 2 * ceil(56.6918 / ucell.lat0 / sqrt(ucell.a2.norm2())) + 1); EXPECT_EQ(vdwd2_test.parameter().period().z, 2 * ceil(56.6918 / ucell.lat0 / sqrt(ucell.a3.norm2())) + 1); @@ -250,11 +260,10 @@ TEST_F(vdwd2Test, WrongVdwType) { input.vdw_method = "d2d3"; testing::internal::CaptureStdout(); - EXPECT_EXIT( vdw::make_vdw(ucell, input); ,::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(vdw::make_vdw(ucell, input);, ::testing::ExitedWithCode(1), ""); std::string output = testing::internal::GetCapturedStdout(); } - TEST_F(vdwd2Test, D2ReadFile) { input.vdw_C6_file = "c6.txt"; @@ -262,9 +271,10 @@ TEST_F(vdwd2Test, D2ReadFile) vdw::Vdwd2 vdwd2_test(ucell); vdwd2_test.parameter().initial_parameters(input); - double Si_C6 = 9.13*1e6 / (ModuleBase::ELECTRONVOLT_SI * ModuleBase::NA) / pow(ModuleBase::BOHR_TO_A, 6)/ ModuleBase::Ry_to_eV; - EXPECT_NEAR(vdwd2_test.parameter().C6().at("Si"), Si_C6,1e-13); - EXPECT_EQ(vdwd2_test.parameter().R0().at("Si"), 1.626/ModuleBase::BOHR_TO_A); + double Si_C6 = 9.13 * 1e6 / (ModuleBase::ELECTRONVOLT_SI * ModuleBase::NA) / pow(ModuleBase::BOHR_TO_A, 6) + / ModuleBase::Ry_to_eV; + EXPECT_NEAR(vdwd2_test.parameter().C6().at("Si"), Si_C6, 1e-13); + EXPECT_EQ(vdwd2_test.parameter().R0().at("Si"), 1.626 / ModuleBase::BOHR_TO_A); } TEST_F(vdwd2Test, D2ReadFileError) @@ -274,8 +284,12 @@ TEST_F(vdwd2Test, D2ReadFileError) vdw::Vdwd2 vdwd2_test(ucell); testing::internal::CaptureStdout(); - EXPECT_EXIT(vdwd2_test.parameter().C6_input(input.vdw_C6_file, input.vdw_C6_unit), ::testing::ExitedWithCode(1), ""); - EXPECT_EXIT(vdwd2_test.parameter().R0_input(input.vdw_R0_file, input.vdw_R0_unit), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(vdwd2_test.parameter().C6_input(input.vdw_C6_file, input.vdw_C6_unit), + ::testing::ExitedWithCode(1), + ""); + EXPECT_EXIT(vdwd2_test.parameter().R0_input(input.vdw_R0_file, input.vdw_R0_unit), + ::testing::ExitedWithCode(1), + ""); std::string output = testing::internal::GetCapturedStdout(); } @@ -286,7 +300,7 @@ TEST_F(vdwd2Test, D2c6UniteVA6) vdwd2_test.parameter().initial_parameters(input); double Si_C6 = 9.23 / pow(ModuleBase::BOHR_TO_A, 6) * ModuleBase::Ry_to_eV; - EXPECT_NEAR(vdwd2_test.parameter().C6().at("Si"), Si_C6,1e-13); + EXPECT_NEAR(vdwd2_test.parameter().C6().at("Si"), Si_C6, 1e-13); } TEST_F(vdwd2Test, D2r0UnitBohr) @@ -305,8 +319,12 @@ TEST_F(vdwd2Test, D2WrongUnit) vdw::Vdwd2 vdwd2_test(ucell); testing::internal::CaptureStdout(); - EXPECT_EXIT(vdwd2_test.parameter().C6_input(input.vdw_C6_file, input.vdw_C6_unit), ::testing::ExitedWithCode(1), ""); - EXPECT_EXIT(vdwd2_test.parameter().R0_input(input.vdw_R0_file, input.vdw_R0_unit), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(vdwd2_test.parameter().C6_input(input.vdw_C6_file, input.vdw_C6_unit), + ::testing::ExitedWithCode(1), + ""); + EXPECT_EXIT(vdwd2_test.parameter().R0_input(input.vdw_R0_file, input.vdw_R0_unit), + ::testing::ExitedWithCode(1), + ""); std::string output = testing::internal::GetCapturedStdout(); } @@ -317,13 +335,13 @@ TEST_F(vdwd2Test, D2RadiusUnitAngstrom) vdw::Vdwd2 vdwd2_test(ucell); vdwd2_test.parameter().initial_parameters(input); - EXPECT_EQ(vdwd2_test.parameter().radius(), 56.6918/ModuleBase::BOHR_TO_A); + EXPECT_EQ(vdwd2_test.parameter().radius(), 56.6918 / ModuleBase::BOHR_TO_A); } TEST_F(vdwd2Test, D2CutoffTypePeriod) { input.vdw_cutoff_type = "period"; - input.vdw_cutoff_period = {3,3,3}; + input.vdw_cutoff_period = {3, 3, 3}; vdw::Vdwd2 vdwd2_test(ucell); vdwd2_test.parameter().initial_parameters(input); @@ -348,7 +366,7 @@ TEST_F(vdwd2Test, D2GetEnergy) auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false)); const double ene = result.energy; - EXPECT_NEAR(ene,-0.034526673470525196,1E-10); + EXPECT_NEAR(ene, -0.034526673470525196, 1E-10); } TEST_F(vdwd2Test, D2GetForce) @@ -359,12 +377,12 @@ TEST_F(vdwd2Test, D2GetForce) ASSERT_TRUE(result.has_force); EXPECT_FALSE(result.has_stress); const std::vector>& force = result.force; - EXPECT_NEAR(force[0].x, -0.00078824525563651242,1e-12); - EXPECT_NEAR(force[0].y, 2.6299822052061785e-08,1e-12); - EXPECT_NEAR(force[0].z, 2.6299822050796364e-08,1e-12); - EXPECT_NEAR(force[1].x, 0.00078824525563651025,1e-12); - EXPECT_NEAR(force[1].y, -2.6299822051961511e-08,1e-12); - EXPECT_NEAR(force[1].z, -2.6299822050777911e-08,1e-12); + EXPECT_NEAR(force[0].x, -0.00078824525563651242, 1e-12); + EXPECT_NEAR(force[0].y, 2.6299822052061785e-08, 1e-12); + EXPECT_NEAR(force[0].z, 2.6299822050796364e-08, 1e-12); + EXPECT_NEAR(force[1].x, 0.00078824525563651025, 1e-12); + EXPECT_NEAR(force[1].y, -2.6299822051961511e-08, 1e-12); + EXPECT_NEAR(force[1].z, -2.6299822050777911e-08, 1e-12); } TEST_F(vdwd2Test, D2GetStress) @@ -375,47 +393,41 @@ TEST_F(vdwd2Test, D2GetStress) ASSERT_TRUE(result.has_force); ASSERT_TRUE(result.has_stress); const ModuleBase::Matrix3& stress = result.stress; - EXPECT_NEAR(stress.e11, -0.00020532319044269705,1e-12); - EXPECT_NEAR(stress.e12, -3.5642821939401251e-08,1e-12); - EXPECT_NEAR(stress.e13, -3.5642821939437223e-08,1e-12); - EXPECT_NEAR(stress.e21, -3.5642821939401251e-08,1e-12); - EXPECT_NEAR(stress.e22, -0.00020491799872060745,1e-12); - EXPECT_NEAR(stress.e23, -4.4369449138539291e-06,1e-12); - EXPECT_NEAR(stress.e31, -3.5642821939437223e-08,1e-12); - EXPECT_NEAR(stress.e32, -4.4369449138539291e-06,1e-12); - EXPECT_NEAR(stress.e33, -0.00020491799872060734,1e-12); + EXPECT_NEAR(stress.e11, -0.00020532319044269705, 1e-12); + EXPECT_NEAR(stress.e12, -3.5642821939401251e-08, 1e-12); + EXPECT_NEAR(stress.e13, -3.5642821939437223e-08, 1e-12); + EXPECT_NEAR(stress.e21, -3.5642821939401251e-08, 1e-12); + EXPECT_NEAR(stress.e22, -0.00020491799872060745, 1e-12); + EXPECT_NEAR(stress.e23, -4.4369449138539291e-06, 1e-12); + EXPECT_NEAR(stress.e31, -3.5642821939437223e-08, 1e-12); + EXPECT_NEAR(stress.e32, -4.4369449138539291e-06, 1e-12); + EXPECT_NEAR(stress.e33, -0.00020491799872060734, 1e-12); } TEST_F(vdwd2Test, EmptyAtomsAreExcluded) { - const vdw::VdwResult reference = - vdw::make_vdw(ucell, input)->evaluate(vdw::VdwRequest(true, true)); + const vdw::VdwResult reference = vdw::make_vdw(ucell, input)->evaluate(vdw::VdwRequest(true, true)); UnitCell ghost_ucell; construct_two_si_with_ghost(ghost_ucell); - const vdw::VdwResult ghost = - vdw::make_vdw(ghost_ucell, input)->evaluate(vdw::VdwRequest(true, true)); + const vdw::VdwResult ghost = vdw::make_vdw(ghost_ucell, input)->evaluate(vdw::VdwRequest(true, true)); expect_empty_atoms_excluded(reference, ghost, ucell.nat); ClearUcell(ghost_ucell); } - - -class vdwd3Test: public testing::Test +class vdwd3Test : public testing::Test { - protected: + protected: UnitCell ucell; Input_para input; - void SetUp(){ - stru_ structure{std::vector{0.5,0.5,0.0,0.5,0.0,0.5,0.0,0.5,0.5}, - std::vector{atomtype_{"Si", - std::vector>{ - {0., 0., 0.}, - {0.3, 0.25, 0.25} - }}}}; - construct_ucell(structure,ucell); + void SetUp() + { + stru_ structure{ + std::vector{0.5, 0.5, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.5}, + std::vector{atomtype_{"Si", std::vector>{{0., 0., 0.}, {0.3, 0.25, 0.25}}}}}; + construct_ucell(structure, ucell); input.dft_functional = "pbe"; input.vdw_method = "d3_0"; @@ -427,7 +439,8 @@ class vdwd3Test: public testing::Test input.vdw_cn_thr = 40; } - void TearDown(){ + void TearDown() + { ClearUcell(ucell); } }; @@ -483,12 +496,12 @@ TEST_F(vdwd3Test, D30GetForce) ASSERT_TRUE(result.has_force); EXPECT_FALSE(result.has_stress); const std::vector>& force = result.force; - EXPECT_NEAR(force[0].x, -0.00079560409324726114,1e-12); - EXPECT_NEAR(force[0].y, 0.0,1e-12); - EXPECT_NEAR(force[0].z, 0.0,1e-12); - EXPECT_NEAR(force[1].x, 0.00079560409324726114,1e-12); - EXPECT_NEAR(force[1].y, 0.0,1e-12); - EXPECT_NEAR(force[1].z, 0.0,1e-12); + EXPECT_NEAR(force[0].x, -0.00079560409324726114, 1e-12); + EXPECT_NEAR(force[0].y, 0.0, 1e-12); + EXPECT_NEAR(force[0].z, 0.0, 1e-12); + EXPECT_NEAR(force[1].x, 0.00079560409324726114, 1e-12); + EXPECT_NEAR(force[1].y, 0.0, 1e-12); + EXPECT_NEAR(force[1].z, 0.0, 1e-12); } TEST_F(vdwd3Test, D30GetStress) @@ -499,27 +512,25 @@ TEST_F(vdwd3Test, D30GetStress) ASSERT_TRUE(result.has_force); ASSERT_TRUE(result.has_stress); const ModuleBase::Matrix3& stress = result.stress; - EXPECT_NEAR(stress.e11, 3.719638454389264e-05,1e-12); - EXPECT_NEAR(stress.e12, 0.0,1e-12); - EXPECT_NEAR(stress.e13, 0.0,1e-12); - EXPECT_NEAR(stress.e21, 0.0,1e-12); - EXPECT_NEAR(stress.e22, 4.5226077638555907e-05,1e-12); - EXPECT_NEAR(stress.e23, -1.2732393110133044e-05,1e-12); - EXPECT_NEAR(stress.e31, 0.0,1e-12); - EXPECT_NEAR(stress.e32, -1.2732393110133044e-05,1e-12); - EXPECT_NEAR(stress.e33, 4.5226077638555961e-05,1e-12); + EXPECT_NEAR(stress.e11, 3.719638454389264e-05, 1e-12); + EXPECT_NEAR(stress.e12, 0.0, 1e-12); + EXPECT_NEAR(stress.e13, 0.0, 1e-12); + EXPECT_NEAR(stress.e21, 0.0, 1e-12); + EXPECT_NEAR(stress.e22, 4.5226077638555907e-05, 1e-12); + EXPECT_NEAR(stress.e23, -1.2732393110133044e-05, 1e-12); + EXPECT_NEAR(stress.e31, 0.0, 1e-12); + EXPECT_NEAR(stress.e32, -1.2732393110133044e-05, 1e-12); + EXPECT_NEAR(stress.e33, 4.5226077638555961e-05, 1e-12); } TEST_F(vdwd3Test, EmptyAtomsAreExcluded) { - const vdw::VdwResult reference = - vdw::make_vdw(ucell, input)->evaluate(vdw::VdwRequest(true, true)); + const vdw::VdwResult reference = vdw::make_vdw(ucell, input)->evaluate(vdw::VdwRequest(true, true)); UnitCell ghost_ucell; construct_two_si_with_ghost(ghost_ucell); - const vdw::VdwResult ghost = - vdw::make_vdw(ghost_ucell, input)->evaluate(vdw::VdwRequest(true, true)); + const vdw::VdwResult ghost = vdw::make_vdw(ghost_ucell, input)->evaluate(vdw::VdwRequest(true, true)); expect_empty_atoms_excluded(reference, ghost, ucell.nat); ClearUcell(ghost_ucell); @@ -543,12 +554,12 @@ TEST_F(vdwd3Test, D3bjGetForce) ASSERT_TRUE(result.has_force); EXPECT_FALSE(result.has_stress); const std::vector>& force = result.force; - EXPECT_NEAR(force[0].x, -0.0025992693513955879,1e-12); - EXPECT_NEAR(force[0].y, 0.0,1e-12); - EXPECT_NEAR(force[0].z, 0.0,1e-12); - EXPECT_NEAR(force[1].x, 0.0025992693513955879,1e-12); - EXPECT_NEAR(force[1].y, 0.0,1e-12); - EXPECT_NEAR(force[1].z, 0.0,1e-12); + EXPECT_NEAR(force[0].x, -0.0025992693513955879, 1e-12); + EXPECT_NEAR(force[0].y, 0.0, 1e-12); + EXPECT_NEAR(force[0].z, 0.0, 1e-12); + EXPECT_NEAR(force[1].x, 0.0025992693513955879, 1e-12); + EXPECT_NEAR(force[1].y, 0.0, 1e-12); + EXPECT_NEAR(force[1].z, 0.0, 1e-12); } TEST_F(vdwd3Test, D3bjGetStress) @@ -560,30 +571,30 @@ TEST_F(vdwd3Test, D3bjGetStress) ASSERT_TRUE(result.has_force); ASSERT_TRUE(result.has_stress); const ModuleBase::Matrix3& stress = result.stress; - EXPECT_NEAR(stress.e11, 0.00014355975918208151,1e-12); - EXPECT_NEAR(stress.e12, 0.0,1e-12); - EXPECT_NEAR(stress.e13, 0.0,1e-12); - EXPECT_NEAR(stress.e21, 0.0,1e-12); - EXPECT_NEAR(stress.e22, 0.00015345761663797616,1e-12); - EXPECT_NEAR(stress.e23, -1.8098500191039634e-05,1e-12); - EXPECT_NEAR(stress.e31, 0.0,1e-12); - EXPECT_NEAR(stress.e32, -1.8098500191039634e-05,1e-12); - EXPECT_NEAR(stress.e33, 0.00015345761663797621,1e-12); + EXPECT_NEAR(stress.e11, 0.00014355975918208151, 1e-12); + EXPECT_NEAR(stress.e12, 0.0, 1e-12); + EXPECT_NEAR(stress.e13, 0.0, 1e-12); + EXPECT_NEAR(stress.e21, 0.0, 1e-12); + EXPECT_NEAR(stress.e22, 0.00015345761663797616, 1e-12); + EXPECT_NEAR(stress.e23, -1.8098500191039634e-05, 1e-12); + EXPECT_NEAR(stress.e31, 0.0, 1e-12); + EXPECT_NEAR(stress.e32, -1.8098500191039634e-05, 1e-12); + EXPECT_NEAR(stress.e33, 0.00015345761663797621, 1e-12); } - -class vdwd3abcTest: public testing::Test +class vdwd3abcTest : public testing::Test { - protected: + protected: UnitCell ucell; Input_para input; - void SetUp(){ + void SetUp() + { stru_ structure{ std::vector{0.75, 0.75, 0.0, 0.75, 0.0, 0.75, 0.0, 0.75, 0.75}, std::vector{atomtype_{"Si", std::vector>{{0., 0., 0.}, {0.3, 0.25, 0.25}}}, atomtype_{"C", std::vector>{{0.5, 0.5, 0.5}}}}}; - construct_ucell(structure,ucell); + construct_ucell(structure, ucell); input.dft_functional = "pbe"; input.vdw_method = "d3_0"; @@ -595,21 +606,19 @@ class vdwd3abcTest: public testing::Test input.vdw_cn_thr = 40; } - void TearDown(){ + void TearDown() + { ClearUcell(ucell); } }; - TEST_F(vdwd3abcTest, InconsistentPseudopotentialXcIsRejected) { input.dft_functional = "default"; ucell.atoms[0].ncpp.xc_func = "PBE"; ucell.atoms[1].ncpp.xc_func = "LDA"; - EXPECT_EXIT(vdw::make_vdw(ucell, input), - ::testing::ExitedWithCode(1), - "XC name automatic inference failed"); + EXPECT_EXIT(vdw::make_vdw(ucell, input), ::testing::ExitedWithCode(1), "XC name automatic inference failed"); } TEST_F(vdwd3abcTest, D30GetEnergy) @@ -628,12 +637,12 @@ TEST_F(vdwd3abcTest, D30GetForce) ASSERT_TRUE(result.has_force); EXPECT_FALSE(result.has_stress); const std::vector>& force = result.force; - EXPECT_NEAR(force[0].x, 0.00016103341708334517,1e-12); - EXPECT_NEAR(force[0].y, 0.00010195935695294143,1e-12); - EXPECT_NEAR(force[0].z, 0.00010195935695294057,1e-12); - EXPECT_NEAR(force[1].x, -0.00062849646016895316,1e-12); - EXPECT_NEAR(force[1].y, -0.00064921932894278223,1e-12); - EXPECT_NEAR(force[1].z, -0.00064921932894278212,1e-12); + EXPECT_NEAR(force[0].x, 0.00016103341708334517, 1e-12); + EXPECT_NEAR(force[0].y, 0.00010195935695294143, 1e-12); + EXPECT_NEAR(force[0].z, 0.00010195935695294057, 1e-12); + EXPECT_NEAR(force[1].x, -0.00062849646016895316, 1e-12); + EXPECT_NEAR(force[1].y, -0.00064921932894278223, 1e-12); + EXPECT_NEAR(force[1].z, -0.00064921932894278212, 1e-12); } TEST_F(vdwd3abcTest, D30GetStress) @@ -644,15 +653,15 @@ TEST_F(vdwd3abcTest, D30GetStress) ASSERT_TRUE(result.has_force); ASSERT_TRUE(result.has_stress); const ModuleBase::Matrix3& stress = result.stress; - EXPECT_NEAR(stress.e11, 1.4721914168426608e-05,1e-12); - EXPECT_NEAR(stress.e12, -1.0755730453993765e-06,1e-12); - EXPECT_NEAR(stress.e13, -1.0755730453993653e-06,1e-12); - EXPECT_NEAR(stress.e21, -1.075573045399377e-06,1e-12); - EXPECT_NEAR(stress.e22, 1.4594782301730673e-05,1e-12); - EXPECT_NEAR(stress.e23, 2.0893396123596713e-07,1e-12); - EXPECT_NEAR(stress.e31, -1.0755730453993645e-06,1e-12); - EXPECT_NEAR(stress.e32, 2.089339612359676e-07,1e-12); - EXPECT_NEAR(stress.e33, 1.4594782301731136e-05,1e-12); + EXPECT_NEAR(stress.e11, 1.4721914168426608e-05, 1e-12); + EXPECT_NEAR(stress.e12, -1.0755730453993765e-06, 1e-12); + EXPECT_NEAR(stress.e13, -1.0755730453993653e-06, 1e-12); + EXPECT_NEAR(stress.e21, -1.075573045399377e-06, 1e-12); + EXPECT_NEAR(stress.e22, 1.4594782301730673e-05, 1e-12); + EXPECT_NEAR(stress.e23, 2.0893396123596713e-07, 1e-12); + EXPECT_NEAR(stress.e31, -1.0755730453993645e-06, 1e-12); + EXPECT_NEAR(stress.e32, 2.089339612359676e-07, 1e-12); + EXPECT_NEAR(stress.e33, 1.4594782301731136e-05, 1e-12); } TEST_F(vdwd3abcTest, D3bjGetEnergy) @@ -673,12 +682,12 @@ TEST_F(vdwd3abcTest, D3bjGetForce) ASSERT_TRUE(result.has_force); EXPECT_FALSE(result.has_stress); const std::vector>& force = result.force; - EXPECT_NEAR(force[0].x, -0.0010624232313481282,1e-12); - EXPECT_NEAR(force[0].y, -0.0010027552655765149,1e-12); - EXPECT_NEAR(force[0].z, -0.0010027552655765149,1e-12); - EXPECT_NEAR(force[1].x, 0.00015483314976063873,1e-12); - EXPECT_NEAR(force[1].y, -0.00010882072675980092,1e-12); - EXPECT_NEAR(force[1].z, -0.00010882072675979906,1e-12); + EXPECT_NEAR(force[0].x, -0.0010624232313481282, 1e-12); + EXPECT_NEAR(force[0].y, -0.0010027552655765149, 1e-12); + EXPECT_NEAR(force[0].z, -0.0010027552655765149, 1e-12); + EXPECT_NEAR(force[1].x, 0.00015483314976063873, 1e-12); + EXPECT_NEAR(force[1].y, -0.00010882072675980092, 1e-12); + EXPECT_NEAR(force[1].z, -0.00010882072675979906, 1e-12); } TEST_F(vdwd3abcTest, D3bjGetStress) @@ -690,25 +699,24 @@ TEST_F(vdwd3abcTest, D3bjGetStress) ASSERT_TRUE(result.has_force); ASSERT_TRUE(result.has_stress); const ModuleBase::Matrix3& stress = result.stress; - EXPECT_NEAR(stress.e11, 3.3773268956804165e-05,1e-12); - EXPECT_NEAR(stress.e12, -5.1343672869104567e-06,1e-12); - EXPECT_NEAR(stress.e13, -5.1343672869104609e-06,1e-12); - EXPECT_NEAR(stress.e21, -5.1343672869104559e-06,1e-12); - EXPECT_NEAR(stress.e22, 3.4293409389506296e-05,1e-12); - EXPECT_NEAR(stress.e23, -4.3785438987372696e-06,1e-12); - EXPECT_NEAR(stress.e31, -5.1343672869104626e-06,1e-12); - EXPECT_NEAR(stress.e32, -4.3785438987372696e-06,1e-12); - EXPECT_NEAR(stress.e33, 3.4293409389506207e-05,1e-12); + EXPECT_NEAR(stress.e11, 3.3773268956804165e-05, 1e-12); + EXPECT_NEAR(stress.e12, -5.1343672869104567e-06, 1e-12); + EXPECT_NEAR(stress.e13, -5.1343672869104609e-06, 1e-12); + EXPECT_NEAR(stress.e21, -5.1343672869104559e-06, 1e-12); + EXPECT_NEAR(stress.e22, 3.4293409389506296e-05, 1e-12); + EXPECT_NEAR(stress.e23, -4.3785438987372696e-06, 1e-12); + EXPECT_NEAR(stress.e31, -5.1343672869104626e-06, 1e-12); + EXPECT_NEAR(stress.e32, -4.3785438987372696e-06, 1e-12); + EXPECT_NEAR(stress.e33, 3.4293409389506207e-05, 1e-12); } TEST_F(vdwd3abcTest, EmptyAtomsAreExcludedWithThreeBodyTerm) { const std::vector methods{"d3_0", "d3_bj"}; - for (const std::string& method : methods) + for (const std::string& method: methods) { input.vdw_method = method; - const vdw::VdwResult reference = - vdw::make_vdw(ucell, input)->evaluate(vdw::VdwRequest(true, true)); + const vdw::VdwResult reference = vdw::make_vdw(ucell, input)->evaluate(vdw::VdwRequest(true, true)); UnitCell ghost_ucell; stru_ structure{ @@ -720,8 +728,7 @@ TEST_F(vdwd3abcTest, EmptyAtomsAreExcludedWithThreeBodyTerm) ghost_ucell.atoms[2].flag_empty_element = true; ghost_ucell.atoms[2].ncpp.psd = "C"; - const vdw::VdwResult ghost = - vdw::make_vdw(ghost_ucell, input)->evaluate(vdw::VdwRequest(true, true)); + const vdw::VdwResult ghost = vdw::make_vdw(ghost_ucell, input)->evaluate(vdw::VdwRequest(true, true)); expect_empty_atoms_excluded(reference, ghost, ucell.nat); ClearUcell(ghost_ucell); } @@ -729,20 +736,18 @@ TEST_F(vdwd3abcTest, EmptyAtomsAreExcludedWithThreeBodyTerm) #ifdef __DFTD4 -class vdwd4Test: public testing::Test +class vdwd4Test : public testing::Test { - protected: + protected: UnitCell ucell; Input_para input; - void SetUp(){ - stru_ structure{std::vector{0.5,0.5,0.0,0.5,0.0,0.5,0.0,0.5,0.5}, - std::vector{atomtype_{"Si", - std::vector>{ - {0., 0., 0.}, - {0.3, 0.25, 0.25} - }}}}; - construct_ucell(structure,ucell); + void SetUp() + { + stru_ structure{ + std::vector{0.5, 0.5, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.5}, + std::vector{atomtype_{"Si", std::vector>{{0., 0., 0.}, {0.3, 0.25, 0.25}}}}}; + construct_ucell(structure, ucell); ucell.atoms[0].ncpp.zv = 4.0; input.nelec = 8.0; @@ -756,7 +761,8 @@ class vdwd4Test: public testing::Test input.vdw_cn_thr = 30; } - void TearDown(){ + void TearDown() + { ClearUcell(ucell); } }; @@ -826,13 +832,11 @@ TEST_F(vdwd4Test, D4GetStress) TEST_F(vdwd4Test, EmptyAtomsAreExcluded) { - const vdw::VdwResult reference = - vdw::make_vdw(ucell, input)->evaluate(vdw::VdwRequest(true, true)); + const vdw::VdwResult reference = vdw::make_vdw(ucell, input)->evaluate(vdw::VdwRequest(true, true)); UnitCell ghost_ucell; construct_two_si_with_ghost(ghost_ucell); - const vdw::VdwResult ghost = - vdw::make_vdw(ghost_ucell, input)->evaluate(vdw::VdwRequest(true, true)); + const vdw::VdwResult ghost = vdw::make_vdw(ghost_ucell, input)->evaluate(vdw::VdwRequest(true, true)); expect_empty_atoms_excluded(reference, ghost, ucell.nat); ClearUcell(ghost_ucell); @@ -887,13 +891,11 @@ TEST_F(vdwd4Test, D4SGetStress) TEST_F(vdwd4Test, EmptyAtomsAreExcludedForD4S) { input.vdw_d4_model = "d4s"; - const vdw::VdwResult reference = - vdw::make_vdw(ucell, input)->evaluate(vdw::VdwRequest(true, true)); + const vdw::VdwResult reference = vdw::make_vdw(ucell, input)->evaluate(vdw::VdwRequest(true, true)); UnitCell ghost_ucell; construct_two_si_with_ghost(ghost_ucell); - const vdw::VdwResult ghost = - vdw::make_vdw(ghost_ucell, input)->evaluate(vdw::VdwRequest(true, true)); + const vdw::VdwResult ghost = vdw::make_vdw(ghost_ucell, input)->evaluate(vdw::VdwRequest(true, true)); expect_empty_atoms_excluded(reference, ghost, ucell.nat); ClearUcell(ghost_ucell); @@ -901,7 +903,7 @@ TEST_F(vdwd4Test, EmptyAtomsAreExcludedForD4S) #endif // __DFTD4 -int main(int argc, char **argv) +int main(int argc, char** argv) { MPI_Init(&argc, &argv); testing::InitGoogleTest(&argc, argv); diff --git a/source/source_hamilt/module_vdw/test/vdwd3_evaluator_test.cpp b/source/source_hamilt/module_vdw/test/vdwd3_evaluator_test.cpp index fd7cd498073..6523746f134 100644 --- a/source/source_hamilt/module_vdw/test/vdwd3_evaluator_test.cpp +++ b/source/source_hamilt/module_vdw/test/vdwd3_evaluator_test.cpp @@ -1,10 +1,10 @@ -#include "source_hamilt/module_vdw/vdwd3_data.h" #include "source_hamilt/module_vdw/vdwd3_evaluator.h" -#include "source_hamilt/module_vdw/vdwd3_parameters.h" + #include "source_hamilt/module_vdw/vdw_xcname.h" +#include "source_hamilt/module_vdw/vdwd3_data.h" +#include "source_hamilt/module_vdw/vdwd3_parameters.h" #include "gtest/gtest.h" - #include #include #include @@ -131,8 +131,7 @@ Structure mindless09() Structure actinides() { Structure structure; - structure.atomic_numbers = {87, 88, 89, 90, 91, 92, 93, 94, 95, - 96, 97, 98, 99, 100, 101, 102, 103}; + structure.atomic_numbers = {87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103}; structure.positions = { {0.98692316414074, 6.12727238368797, -6.67861597188102}, {3.63898862390869, 5.12109301182962, 3.01908613326278}, @@ -159,8 +158,7 @@ Structure x23_acetic() { Structure structure; structure.atomic_numbers = { - 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, - 6, 6, 6, 6, 6, 6, 6, 6, 8, 8, 8, 8, 8, 8, 8, 8, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 6, 6, 6, 6, 6, 6, 6, 6, 8, 8, 8, 8, 8, 8, 8, 8, }; structure.positions = { {7.90377778530531, 7.12936818903081, 9.56314593193451}, @@ -196,9 +194,7 @@ Structure x23_acetic() {18.50211546160091, 6.44377570116260, 1.93885857925242}, {6.36100571503078, 0.93994956853816, 7.33629315436038}, }; - structure.lattice = {{{24.86304558760325, 0.0, 0.0}, - {0.0, 7.38360999643241, 0.0}, - {0.0, 0.0, 10.79482606446971}}}; + structure.lattice = {{{24.86304558760325, 0.0, 0.0}, {0.0, 7.38360999643241, 0.0}, {0.0, 0.0, 10.79482606446971}}}; structure.periodic = {{true, true, true}}; return structure; } @@ -209,12 +205,9 @@ TEST(VdwXcName, NormalizesAbacusAndLibxcSpellings) { EXPECT_EQ(vdw::normalize_xc_name("PBE"), "pbe"); EXPECT_EQ(vdw::normalize_xc_name(" XC_HYB_GGA_XC_PBEH "), "hyb_gga_xc_pbeh"); - EXPECT_EQ(vdw::normalize_xc_name("GGA_X_PBE+GGA_C_PBE"), - "gga_x_pbe:gga_c_pbe"); - EXPECT_EQ(vdw::normalize_xc_name("XC_GGA_X_PBE + XC_GGA_C_PBE"), - "gga_x_pbe:gga_c_pbe"); - EXPECT_EQ(vdw::normalize_xc_name("XC_GGA_X_PBE_R:XC_GGA_C_PBE"), - "gga_x_pbe_r:gga_c_pbe"); + EXPECT_EQ(vdw::normalize_xc_name("GGA_X_PBE+GGA_C_PBE"), "gga_x_pbe:gga_c_pbe"); + EXPECT_EQ(vdw::normalize_xc_name("XC_GGA_X_PBE + XC_GGA_C_PBE"), "gga_x_pbe:gga_c_pbe"); + EXPECT_EQ(vdw::normalize_xc_name("XC_GGA_X_PBE_R:XC_GGA_C_PBE"), "gga_x_pbe_r:gga_c_pbe"); } TEST(D3Parameters, PbeAndLibxcAliases) @@ -229,10 +222,7 @@ TEST(D3Parameters, PbeAndLibxcAliases) EXPECT_DOUBLE_EQ(zero.rs8, 1.0); Parameters rational; - ASSERT_TRUE(vdw::d3::lookup_parameters("GGA_X_PBE+GGA_C_PBE", - Damping::Rational, - rational, - canonical)); + ASSERT_TRUE(vdw::d3::lookup_parameters("GGA_X_PBE+GGA_C_PBE", Damping::Rational, rational, canonical)); EXPECT_EQ(canonical, "pbe"); EXPECT_DOUBLE_EQ(rational.s8, 0.7875); EXPECT_DOUBLE_EQ(rational.a1, 0.4289); @@ -242,52 +232,31 @@ TEST(D3Parameters, PbeAndLibxcAliases) EXPECT_EQ(canonical, "hse06"); EXPECT_DOUBLE_EQ(rational.a1, 0.383); - ASSERT_TRUE(vdw::d3::lookup_parameters("XC_GGA_X_B88+XC_GGA_C_OP_B88", - Damping::Rational, - rational, - canonical)); + ASSERT_TRUE(vdw::d3::lookup_parameters("XC_GGA_X_B88+XC_GGA_C_OP_B88", Damping::Rational, rational, canonical)); EXPECT_EQ(canonical, "bop"); EXPECT_DOUBLE_EQ(rational.a1, 0.487); - ASSERT_TRUE(vdw::d3::lookup_parameters("XC_GGA_X_PBE_R+XC_GGA_C_PBE", - Damping::Rational, - rational, - canonical)); + ASSERT_TRUE(vdw::d3::lookup_parameters("XC_GGA_X_PBE_R+XC_GGA_C_PBE", Damping::Rational, rational, canonical)); EXPECT_EQ(canonical, "revpbe"); // Legacy ABACUS LibXC spellings: the numerical parameter set is stored // under `wb97x` by s-dftd3, but only for the corresponding damping form. - ASSERT_TRUE(vdw::d3::lookup_parameters("XC_HYB_GGA_XC_WB97X_V", - Damping::Rational, - rational, - canonical)); + ASSERT_TRUE(vdw::d3::lookup_parameters("XC_HYB_GGA_XC_WB97X_V", Damping::Rational, rational, canonical)); EXPECT_EQ(canonical, "wb97xv"); EXPECT_DOUBLE_EQ(rational.s8, 0.2641); EXPECT_DOUBLE_EQ(rational.a1, 0.0); EXPECT_DOUBLE_EQ(rational.a2, 5.4959); - ASSERT_TRUE(vdw::d3::lookup_parameters("XC_HYB_GGA_XC_WB97X_D3", - Damping::Zero, - zero, - canonical)); + ASSERT_TRUE(vdw::d3::lookup_parameters("XC_HYB_GGA_XC_WB97X_D3", Damping::Zero, zero, canonical)); EXPECT_EQ(canonical, "wb97xd3"); EXPECT_DOUBLE_EQ(zero.s8, 1.0); EXPECT_DOUBLE_EQ(zero.rs6, 1.281); EXPECT_DOUBLE_EQ(zero.rs8, 1.094); - EXPECT_FALSE(vdw::d3::lookup_parameters("XC_HYB_GGA_XC_WB97X_V", - Damping::Zero, - zero, - canonical)); - EXPECT_FALSE(vdw::d3::lookup_parameters("XC_HYB_GGA_XC_WB97X_D3", - Damping::Rational, - rational, - canonical)); - - EXPECT_FALSE(vdw::d3::lookup_parameters("definitely-not-an-xc", - Damping::Rational, - rational, - canonical)); + EXPECT_FALSE(vdw::d3::lookup_parameters("XC_HYB_GGA_XC_WB97X_V", Damping::Zero, zero, canonical)); + EXPECT_FALSE(vdw::d3::lookup_parameters("XC_HYB_GGA_XC_WB97X_D3", Damping::Rational, rational, canonical)); + + EXPECT_FALSE(vdw::d3::lookup_parameters("definitely-not-an-xc", Damping::Rational, rational, canonical)); } TEST(D3Data, ReferencePackingAndSymmetry) @@ -318,17 +287,13 @@ TEST(D3Data, ReferencePackingAndSymmetry) { for (int reference_j = 0; reference_j < references_j; ++reference_j) { - const double c6 = vdw::d3::data::reference_c6(atomic_number_i, - reference_i, - atomic_number_j, - reference_j); + const double c6 + = vdw::d3::data::reference_c6(atomic_number_i, reference_i, atomic_number_j, reference_j); EXPECT_TRUE(std::isfinite(c6)); EXPECT_GE(c6, 0.0); - EXPECT_DOUBLE_EQ(c6, - vdw::d3::data::reference_c6(atomic_number_j, - reference_j, - atomic_number_i, - reference_i)); + EXPECT_DOUBLE_EQ( + c6, + vdw::d3::data::reference_c6(atomic_number_j, reference_j, atomic_number_i, reference_i)); } } } @@ -345,8 +310,7 @@ TEST(D3Evaluator, MatchesSdftd3Mindless01) parameters.s9 = 0.0; // Exact value from s-dftd3 v1.5.0 test/unit/test_dftd3.f90. - EXPECT_DOUBLE_EQ(energy(mindless01(), parameters, Cutoffs()), - -1.7882220155186028e-2); + EXPECT_DOUBLE_EQ(energy(mindless01(), parameters, Cutoffs()), -1.7882220155186028e-2); } TEST(D3Evaluator, MatchesSdftd3Actinides) @@ -357,8 +321,7 @@ TEST(D3Evaluator, MatchesSdftd3Actinides) parameters.s9 = 0.0; // Exact value from s-dftd3 v1.5.0 test/unit/test_dftd3.f90. - EXPECT_DOUBLE_EQ(energy(actinides(), parameters, Cutoffs()), - -1.4131143363689097e-1); + EXPECT_DOUBLE_EQ(energy(actinides(), parameters, Cutoffs()), -1.4131143363689097e-1); } TEST(D3Evaluator, MatchesSdftd3PeriodicX23Acetic) @@ -374,23 +337,18 @@ TEST(D3Evaluator, MatchesSdftd3PeriodicX23Acetic) cutoffs.disp3 = 15.0; // Exact value from s-dftd3 v1.5.0 test/unit/test_periodic_3d.f90. - EXPECT_DOUBLE_EQ(energy(x23_acetic(), parameters, cutoffs), - -6.6732836815486210e-2); + EXPECT_DOUBLE_EQ(energy(x23_acetic(), parameters, cutoffs), -6.6732836815486210e-2); } TEST(D3Evaluator, MatchesSdftd3AtmMindless17) { Parameters parameters; std::string canonical; - ASSERT_TRUE(vdw::d3::lookup_parameters("dsd-blyp", - Damping::Rational, - parameters, - canonical)); + ASSERT_TRUE(vdw::d3::lookup_parameters("dsd-blyp", Damping::Rational, parameters, canonical)); parameters.s9 = 1.0; // Exact value from s-dftd3 v1.5.0 test/unit/test_dftd3.f90. - EXPECT_DOUBLE_EQ(energy(mindless17(), parameters, Cutoffs()), - -1.3592755832923201e-2); + EXPECT_DOUBLE_EQ(energy(mindless17(), parameters, Cutoffs()), -1.3592755832923201e-2); } TEST(D3Evaluator, MatchesSdftd3ZeroDampingMindless09) @@ -401,19 +359,17 @@ TEST(D3Evaluator, MatchesSdftd3ZeroDampingMindless09) parameters.s9 = 0.0; // Exact value from s-dftd3 v1.5.0 test/unit/test_dftd3.f90. - EXPECT_DOUBLE_EQ(energy(mindless09(), parameters, Cutoffs()), - -2.0178760785797962e-2); + EXPECT_DOUBLE_EQ(energy(mindless09(), parameters, Cutoffs()), -2.0178760785797962e-2); } TEST(D3Evaluator, HighCoordinationMgGradientAndVirial) { - // Two Mg atoms in a slightly distorted FCC supercell. Their D3 CN is near 12, - // well above the largest Mg reference CN (1.9496). + // Two Mg atoms in a slightly distorted FCC supercell. Their D3 CN is near + // 12, well above the largest Mg reference CN (1.9496). const double half_diagonal = 5.3 / std::sqrt(2.0); Structure structure; structure.atomic_numbers = {12, 12}; - structure.positions = {{0.0, 0.0, 0.0}, - {0.06, half_diagonal - 0.09, half_diagonal + 0.13}}; + structure.positions = {{0.0, 0.0, 0.0}, {0.06, half_diagonal - 0.09, half_diagonal + 0.13}}; structure.lattice = {{{0.0, 2.0 * half_diagonal, 2.0 * half_diagonal}, {half_diagonal, 0.0, half_diagonal}, {half_diagonal, half_diagonal, 0.0}}}; @@ -432,11 +388,9 @@ TEST(D3Evaluator, HighCoordinationMgGradientAndVirial) std::string error; ASSERT_TRUE(vdw::d3::evaluate(structure, parameters, cutoffs, true, analytic, error)) << error; ASSERT_TRUE(std::isfinite(analytic.energy)); - for (const Vec3& gradient : analytic.gradient) + for (const Vec3& gradient: analytic.gradient) { - const double magnitude = std::sqrt(gradient.x * gradient.x - + gradient.y * gradient.y - + gradient.z * gradient.z); + const double magnitude = std::sqrt(gradient.x * gradient.x + gradient.y * gradient.y + gradient.z * gradient.z); EXPECT_TRUE(std::isfinite(magnitude)); EXPECT_LT(magnitude, 1.0); } @@ -450,8 +404,8 @@ TEST(D3Evaluator, HighCoordinationMgGradientAndVirial) Structure minus = structure; coordinate(plus.positions[atom], component) += step; coordinate(minus.positions[atom], component) -= step; - const double numerical = (energy(plus, parameters, cutoffs) - - energy(minus, parameters, cutoffs)) / (2.0 * step); + const double numerical + = (energy(plus, parameters, cutoffs) - energy(minus, parameters, cutoffs)) / (2.0 * step); EXPECT_NEAR(coordinate(analytic.gradient[atom], component), numerical, 1.0e-8); } } @@ -470,8 +424,8 @@ TEST(D3Evaluator, HighCoordinationMgGradientAndVirial) plus.lattice[vector].x += displacement; minus.lattice[vector].x -= displacement; } - const double numerical_virial = (energy(plus, parameters, cutoffs) - - energy(minus, parameters, cutoffs)) / (2.0 * step); + const double numerical_virial + = (energy(plus, parameters, cutoffs) - energy(minus, parameters, cutoffs)) / (2.0 * step); EXPECT_NEAR(analytic.virial.value[0][0], numerical_virial, 1.0e-8); } @@ -507,12 +461,10 @@ TEST(D3Evaluator, SmoothCutoffGradientAndVirial) Structure minus = structure; coordinate(plus.positions[atom], component) += step; coordinate(minus.positions[atom], component) -= step; - const double numerical = (energy(plus, parameters, cutoffs) - - energy(minus, parameters, cutoffs)) - / (2.0 * step); - maximum_gradient_error = std::max( - maximum_gradient_error, - std::abs(numerical - coordinate(analytic.gradient[atom], component))); + const double numerical + = (energy(plus, parameters, cutoffs) - energy(minus, parameters, cutoffs)) / (2.0 * step); + maximum_gradient_error = std::max(maximum_gradient_error, + std::abs(numerical - coordinate(analytic.gradient[atom], component))); } } EXPECT_LT(maximum_gradient_error, 1.0e-9); @@ -530,12 +482,10 @@ TEST(D3Evaluator, SmoothCutoffGradientAndVirial) coordinate(plus.positions[atom], strain_row) += displacement; coordinate(minus.positions[atom], strain_row) -= displacement; } - const double numerical = (energy(plus, parameters, cutoffs) - - energy(minus, parameters, cutoffs)) - / (2.0 * step); - maximum_virial_error = std::max( - maximum_virial_error, - std::abs(numerical - analytic.virial.value[strain_row][strain_column])); + const double numerical + = (energy(plus, parameters, cutoffs) - energy(minus, parameters, cutoffs)) / (2.0 * step); + maximum_virial_error = std::max(maximum_virial_error, + std::abs(numerical - analytic.virial.value[strain_row][strain_column])); } } EXPECT_LT(maximum_virial_error, 1.0e-8); @@ -575,12 +525,10 @@ TEST(D3Evaluator, PeriodicSmoothCutoffVirial) Structure minus = structure; coordinate(plus.positions[atom], component) += step; coordinate(minus.positions[atom], component) -= step; - const double numerical = (energy(plus, parameters, cutoffs) - - energy(minus, parameters, cutoffs)) - / (2.0 * step); - maximum_gradient_error = std::max( - maximum_gradient_error, - std::abs(numerical - coordinate(analytic.gradient[atom], component))); + const double numerical + = (energy(plus, parameters, cutoffs) - energy(minus, parameters, cutoffs)) / (2.0 * step); + maximum_gradient_error = std::max(maximum_gradient_error, + std::abs(numerical - coordinate(analytic.gradient[atom], component))); } } EXPECT_LT(maximum_gradient_error, 1.0e-8); @@ -604,12 +552,10 @@ TEST(D3Evaluator, PeriodicSmoothCutoffVirial) coordinate(plus.lattice[vector], strain_row) += displacement; coordinate(minus.lattice[vector], strain_row) -= displacement; } - const double numerical = (energy(plus, parameters, cutoffs) - - energy(minus, parameters, cutoffs)) - / (2.0 * step); - maximum_virial_error = std::max( - maximum_virial_error, - std::abs(numerical - analytic.virial.value[strain_row][strain_column])); + const double numerical + = (energy(plus, parameters, cutoffs) - energy(minus, parameters, cutoffs)) / (2.0 * step); + maximum_virial_error = std::max(maximum_virial_error, + std::abs(numerical - analytic.virial.value[strain_row][strain_column])); } } EXPECT_LT(maximum_virial_error, 1.0e-8); @@ -636,8 +582,7 @@ TEST(D3Evaluator, SharpCutoffCompatibilityAndWidthValidation) Cutoffs equal_width = sharp; equal_width.width2 = equal_width.disp2; equal_width.width3 = equal_width.disp3; - EXPECT_DOUBLE_EQ(energy(structure, parameters, sharp), - energy(structure, parameters, equal_width)); + EXPECT_DOUBLE_EQ(energy(structure, parameters, sharp), energy(structure, parameters, equal_width)); Cutoffs invalid = sharp; invalid.width2 = invalid.disp2 + 1.0; diff --git a/source/source_hamilt/module_vdw/vdw.cpp b/source/source_hamilt/module_vdw/vdw.cpp index de65497d245..2bcad9de68d 100644 --- a/source/source_hamilt/module_vdw/vdw.cpp +++ b/source/source_hamilt/module_vdw/vdw.cpp @@ -1,17 +1,17 @@ -#include -#include - #include "vdw.h" + +#include "source_base/tool_quit.h" #include "vdwd2.h" #include "vdwd3.h" -#include "source_base/tool_quit.h" + +#include +#include #ifdef __DFTD4 #include "vdwd4.h" #endif -std::string parse_xcname(const std::string &xc_input, - const std::vector &xc_psp) +std::string parse_xcname(const std::string& xc_input, const std::vector& xc_psp) { if (xc_input != "default") { @@ -20,41 +20,41 @@ std::string parse_xcname(const std::string &xc_input, if (xc_psp.size() <= 0) { - ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::parse_xcname", - "XC name automatic inference failed: no pseudopotential files are found"); + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::parse_xcname", + "XC name automatic inference failed: no " + "pseudopotential files are found"); } std::vector xc_psp_uniq = xc_psp; std::sort(xc_psp_uniq.begin(), xc_psp_uniq.end()); auto last = std::unique(xc_psp_uniq.begin(), xc_psp_uniq.end()); xc_psp_uniq.erase(last, xc_psp_uniq.end()); - + if (xc_psp_uniq.size() > 1) { - ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::parse_xcname", - "XC name automatic inference failed: inconsistency in XC names is found" - " in the pseudopotential files"); + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::parse_xcname", + "XC name automatic inference failed: " + "inconsistency in XC names is found" + " in the pseudopotential files"); } const std::string xc = xc_psp_uniq[0]; std::cout << " ***WARNING*** ModuleHamiltGeneral::ModuleVDW::parse_xcname: " - << "XC name is automatically inferred from pseudopotential as `" - << xc << "`" << std::endl; + << "XC name is automatically inferred from pseudopotential as `" << xc << "`" << std::endl; return xc; } namespace vdw { -std::unique_ptr make_vdw(const UnitCell &ucell, - const Input_para &input, - std::ofstream* plog) +std::unique_ptr make_vdw(const UnitCell& ucell, const Input_para& input, std::ofstream* plog) { - // NOTE: the following lines are incorrect! + // NOTE: the following lines are incorrect! // NOTE: VDW interaction exists between images even if there is only one - // NOTE: atom in cell. See issue#5401: https://github.com/deepmodeling/abacus-develop/issues/5401 - // if (ucell.nat < 2 && input.vdw_method != "none") + // NOTE: atom in cell. See issue#5401: + // https://github.com/deepmodeling/abacus-develop/issues/5401 if (ucell.nat + // < 2 && input.vdw_method != "none") // { - // ModuleBase::WARNING("VDW", "Only one atom in this system, and will not do the calculation of VDW"); - // return nullptr; + // ModuleBase::WARNING("VDW", "Only one atom in this system, and will + // not do the calculation of VDW"); return nullptr; // } if (input.vdw_method == "d2") { @@ -70,10 +70,7 @@ std::unique_ptr make_vdw(const UnitCell &ucell, { xc_psp[it] = ucell.atoms[it].ncpp.xc_func; } - return make_unique_compat(ucell, - parse_xcname(input.dft_functional, xc_psp), - input, - plog); + return make_unique_compat(ucell, parse_xcname(input.dft_functional, xc_psp), input, plog); } else if (input.vdw_method == "d4") { @@ -99,8 +96,8 @@ std::unique_ptr make_vdw(const UnitCell &ucell, } else if (input.vdw_method != "none") { - ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::make_vdw", - "Unrecognized Van der Waals correction method: " + input.vdw_method); + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::make_vdw", + "Unrecognized Van der Waals correction method: " + input.vdw_method); return nullptr; } return nullptr; // "none" method diff --git a/source/source_hamilt/module_vdw/vdw.h b/source/source_hamilt/module_vdw/vdw.h index 48025cfea47..0bc8e5146e2 100644 --- a/source/source_hamilt/module_vdw/vdw.h +++ b/source/source_hamilt/module_vdw/vdw.h @@ -1,14 +1,14 @@ #ifndef VDW_H #define VDW_H -#include +#include "source_cell/unitcell.h" +#include "source_io/module_parameter/input_parameter.h" + #include +#include #include #include -#include "source_cell/unitcell.h" -#include "source_io/module_parameter/input_parameter.h" - namespace vdw { @@ -20,7 +20,9 @@ std::unique_ptr make_unique_compat(Args&&... args) struct VdwRequest { - VdwRequest(const bool force_in, const bool stress_in) : force(force_in), stress(stress_in) {} + VdwRequest(const bool force_in, const bool stress_in) : force(force_in), stress(stress_in) + { + } bool force; bool stress; @@ -43,7 +45,9 @@ struct VdwResult class Vdw { public: - Vdw(const UnitCell& unit_in) : ucell_(unit_in) {} + Vdw(const UnitCell& unit_in) : ucell_(unit_in) + { + } virtual ~Vdw() = default; @@ -68,9 +72,7 @@ class Vdw * @param plog optional, for logging the parameter setting process * @return std::unique_ptr */ -std::unique_ptr make_vdw(const UnitCell& ucell, - const Input_para& input, - std::ofstream* plog = nullptr); +std::unique_ptr make_vdw(const UnitCell& ucell, const Input_para& input, std::ofstream* plog = nullptr); } // namespace vdw diff --git a/source/source_hamilt/module_vdw/vdw_parameters.h b/source/source_hamilt/module_vdw/vdw_parameters.h index 630ff3d9e7f..4972efd8801 100644 --- a/source/source_hamilt/module_vdw/vdw_parameters.h +++ b/source/source_hamilt/module_vdw/vdw_parameters.h @@ -14,8 +14,14 @@ class VdwParameters VdwParameters() = default; virtual ~VdwParameters() = default; - inline const std::string &model() const { return model_; } - inline const ModuleBase::Vector3 &period() const { return period_; }; + inline const std::string& model() const + { + return model_; + } + inline const ModuleBase::Vector3& period() const + { + return period_; + }; protected: std::string model_; diff --git a/source/source_hamilt/module_vdw/vdw_xcname.cpp b/source/source_hamilt/module_vdw/vdw_xcname.cpp index 93b624949b8..5741f3bb2b6 100644 --- a/source/source_hamilt/module_vdw/vdw_xcname.cpp +++ b/source/source_hamilt/module_vdw/vdw_xcname.cpp @@ -20,8 +20,8 @@ std::string trim_ascii(const std::string& value) } const auto last = std::find_if_not(value.rbegin(), value.rend(), [](unsigned char character) { - return std::isspace(character) != 0; - }).base(); + return std::isspace(character) != 0; + }).base(); return std::string(first, last); } @@ -52,8 +52,7 @@ std::string normalize_xc_name(const std::string& input) return normalize_component(input); } - return normalize_component(input.substr(0, separator)) + ":" - + normalize_component(input.substr(separator + 1)); + return normalize_component(input.substr(0, separator)) + ":" + normalize_component(input.substr(separator + 1)); } } // namespace vdw diff --git a/source/source_hamilt/module_vdw/vdwd2.cpp b/source/source_hamilt/module_vdw/vdwd2.cpp index 9489ea62a75..d2f9ad47df3 100644 --- a/source/source_hamilt/module_vdw/vdwd2.cpp +++ b/source/source_hamilt/module_vdw/vdwd2.cpp @@ -5,6 +5,7 @@ //========================================================== #include "vdwd2.h" + #include "source_base/timer.h" #include @@ -49,8 +50,7 @@ void Vdwd2::evaluate_impl(const VdwRequest& request, VdwResult& result) } const double tmp_factor = C6_product / pow(r_sqr, 3) / r / tmp_damp_recip - * (-6.0 / r - + tmp_exp / tmp_damp_recip * para_.damping() / R0_sum); + * (-6.0 / r + tmp_exp / tmp_damp_recip * para_.damping() / R0_sum); if (request.force) { diff --git a/source/source_hamilt/module_vdw/vdwd2.h b/source/source_hamilt/module_vdw/vdwd2.h index 43d93615db5..6671ebcbe5b 100644 --- a/source/source_hamilt/module_vdw/vdwd2.h +++ b/source/source_hamilt/module_vdw/vdwd2.h @@ -17,12 +17,20 @@ class Vdwd2 : public Vdw { public: - Vdwd2(const UnitCell& unit_in) : Vdw(unit_in) {} + Vdwd2(const UnitCell& unit_in) : Vdw(unit_in) + { + } ~Vdwd2() = default; - Vdwd2Parameters& parameter() { return para_; } - const Vdwd2Parameters& parameter() const { return para_; } + Vdwd2Parameters& parameter() + { + return para_; + } + const Vdwd2Parameters& parameter() const + { + return para_; + } private: Vdwd2Parameters para_; @@ -70,8 +78,7 @@ class Vdwd2 : public Vdw { for (ilat_loop.z = -zidx; ilat_loop.z <= zidx; ++ilat_loop.z) { - if ((!(ilat_loop.x || ilat_loop.y || ilat_loop.z)) && (it1 == it2) - && (ia1 == ia2)) + if ((!(ilat_loop.x || ilat_loop.y || ilat_loop.z)) && (it1 == it2) && (ia1 == ia2)) { continue; } diff --git a/source/source_hamilt/module_vdw/vdwd2_parameters.cpp b/source/source_hamilt/module_vdw/vdwd2_parameters.cpp index 3ac48e18199..bff55dc5b9b 100644 --- a/source/source_hamilt/module_vdw/vdwd2_parameters.cpp +++ b/source/source_hamilt/module_vdw/vdwd2_parameters.cpp @@ -9,7 +9,7 @@ namespace vdw { -void Vdwd2Parameters::initial_parameters(const Input_para &input, std::ofstream* plog) +void Vdwd2Parameters::initial_parameters(const Input_para& input, std::ofstream* plog) { scaling_ = std::stod(input.vdw_s6); damping_ = input.vdw_d; @@ -30,7 +30,7 @@ void Vdwd2Parameters::initial_parameters(const Input_para &input, std::ofstream* } } -void Vdwd2Parameters::initset(const UnitCell &ucell) +void Vdwd2Parameters::initset(const UnitCell& ucell) { if (model_ == "radius") { @@ -40,7 +40,7 @@ void Vdwd2Parameters::initset(const UnitCell &ucell) } } -void Vdwd2Parameters::C6_input(const std::string &file, const std::string &unit) +void Vdwd2Parameters::C6_input(const std::string& file, const std::string& unit) { C6_ = C6_default_; if (file != "default") @@ -55,7 +55,7 @@ void Vdwd2Parameters::C6_input(const std::string &file, const std::string &unit) C6_[element] = value; ifs.close(); } - for (auto &c6: C6_) + for (auto& c6: C6_) { if (unit == "Jnm6/mol") c6.second *= 1e6 / (ModuleBase::ELECTRONVOLT_SI * ModuleBase::NA) / pow(ModuleBase::BOHR_TO_A, 6) @@ -67,7 +67,7 @@ void Vdwd2Parameters::C6_input(const std::string &file, const std::string &unit) } } -void Vdwd2Parameters::R0_input(const std::string &file, const std::string &unit) +void Vdwd2Parameters::R0_input(const std::string& file, const std::string& unit) { R0_ = R0_default_; if (file != "default") @@ -82,7 +82,7 @@ void Vdwd2Parameters::R0_input(const std::string &file, const std::string &unit) R0_[element] = value; ifs.close(); } - for (auto &r0: R0_) + for (auto& r0: R0_) { if (unit == "A") r0.second /= ModuleBase::BOHR_TO_A; diff --git a/source/source_hamilt/module_vdw/vdwd2_parameters.h b/source/source_hamilt/module_vdw/vdwd2_parameters.h index a22fe350320..87d93273337 100644 --- a/source/source_hamilt/module_vdw/vdwd2_parameters.h +++ b/source/source_hamilt/module_vdw/vdwd2_parameters.h @@ -7,11 +7,12 @@ #ifndef VDWD2_PARAMETERS_H #define VDWD2_PARAMETERS_H -#include -#include "source_io/module_parameter/parameter.h" #include "source_cell/unitcell.h" +#include "source_io/module_parameter/parameter.h" #include "vdw_parameters.h" +#include + namespace vdw { @@ -30,26 +31,42 @@ class Vdwd2Parameters : public VdwParameters ~Vdwd2Parameters() = default; - void C6_input(const std::string &file, const std::string &unit); - void R0_input(const std::string &file, const std::string &unit); + void C6_input(const std::string& file, const std::string& unit); + void R0_input(const std::string& file, const std::string& unit); - void initset(const UnitCell &ucell); // init sets of vdwd2 once this correction is called + void initset(const UnitCell& ucell); // init sets of vdwd2 once this + // correction is called /** * @brief initial parameters of Vdwd2 with INPUT file - * + * * @param input Parameter instance - * @param plog optional, for logging the parameter setting process (not implemented) + * @param plog optional, for logging the parameter setting process (not + * implemented) */ - void initial_parameters(const Input_para &input, - std::ofstream* plog = nullptr); + void initial_parameters(const Input_para& input, std::ofstream* plog = nullptr); - inline const std::map C6() const { return C6_; } - inline const std::map R0() const { return R0_; } - inline double damping() const { return damping_; } - inline double scaling() const { return scaling_; } + inline const std::map C6() const + { + return C6_; + } + inline const std::map R0() const + { + return R0_; + } + inline double damping() const + { + return damping_; + } + inline double scaling() const + { + return scaling_; + } /// @brief cutoff radius of the D2 summation, in Bohr - inline double radius() const { return radius_; } + inline double radius() const + { + return radius_; + } private: double scaling_; diff --git a/source/source_hamilt/module_vdw/vdwd3.cpp b/source/source_hamilt/module_vdw/vdwd3.cpp index 923ea619300..523d5231ac8 100644 --- a/source/source_hamilt/module_vdw/vdwd3.cpp +++ b/source/source_hamilt/module_vdw/vdwd3.cpp @@ -1,11 +1,11 @@ #include "vdwd3.h" -#include "vdwd3_evaluator.h" -#include "vdwd3_parameters.h" #include "source_base/constants.h" #include "source_base/element_name.h" #include "source_base/timer.h" #include "source_base/tool_quit.h" +#include "vdwd3_evaluator.h" +#include "vdwd3_parameters.h" #include #include @@ -60,10 +60,7 @@ const char* damping_name(d3::Damping damping) } // namespace -Vdwd3::Vdwd3(const UnitCell& unit_in, - const std::string& xc_name, - const Input_para& input, - std::ofstream* plog) +Vdwd3::Vdwd3(const UnitCell& unit_in, const std::string& xc_name, const Input_para& input, std::ofstream* plog) : Vdw(unit_in) { if (input.vdw_cutoff_type != "radius") @@ -86,23 +83,18 @@ Vdwd3::Vdwd3(const UnitCell& unit_in, ModuleBase::WARNING_QUIT("Vdwd3::Vdwd3", "Unsupported DFT-D3 method: " + input.vdw_method); } - const bool all_custom = input.vdw_s6 != "default" - && input.vdw_s8 != "default" - && input.vdw_a1 != "default" + const bool all_custom = input.vdw_s6 != "default" && input.vdw_s8 != "default" && input.vdw_a1 != "default" && input.vdw_a2 != "default"; if (all_custom) { canonical_method_ = d3::canonicalize_method_name(xc_name); } - else if (!d3::lookup_parameters(xc_name, - parameters_.damping, - parameters_, - canonical_method_)) + else if (!d3::lookup_parameters(xc_name, parameters_.damping, parameters_, canonical_method_)) { ModuleBase::WARNING_QUIT("Vdwd3::Vdwd3", - "No s-dftd3 damping parameters found for XC functional '" - + xc_name + "'. Define vdw_s6, vdw_s8, vdw_a1 and vdw_a2 " - "together to use a fully custom parameter set."); + "No s-dftd3 damping parameters found for XC functional '" + xc_name + + "'. Define vdw_s6, vdw_s8, vdw_a1 and vdw_a2 " + "together to use a fully custom parameter set."); } if (input.vdw_s6 != "default") @@ -169,9 +161,7 @@ void Vdwd3::write_parameters(std::ofstream* plog) const *plog << "\nDFT-D3 parameters (s-dftd3 v1.5.0 numerical specification)\n" << "XC functional: " << canonical_method_ << '\n' << "damping: " << damping_name(parameters_.damping) << '\n' - << std::setprecision(10) - << "s6=" << parameters_.s6 << " s8=" << parameters_.s8 - << " s9=" << parameters_.s9; + << std::setprecision(10) << "s6=" << parameters_.s6 << " s8=" << parameters_.s8 << " s9=" << parameters_.s9; if (parameters_.damping == d3::Damping::Rational) { *plog << " a1=" << parameters_.a1 << " a2=" << parameters_.a2; @@ -180,10 +170,8 @@ void Vdwd3::write_parameters(std::ofstream* plog) const { *plog << " rs6=" << parameters_.rs6 << " rs8=" << parameters_.rs8; } - *plog << "\ncutoffs (Bohr): disp2=" << cutoffs_.disp2 - << " disp3=" << cutoffs_.disp3 << " CN=" << cutoffs_.cn - << " smooth_width_2b=" << cutoffs_.width2 - << " smooth_width_3b=" << cutoffs_.width3 << std::endl; + *plog << "\ncutoffs (Bohr): disp2=" << cutoffs_.disp2 << " disp3=" << cutoffs_.disp3 << " CN=" << cutoffs_.cn + << " smooth_width_2b=" << cutoffs_.width2 << " smooth_width_3b=" << cutoffs_.width3 << std::endl; } d3::Structure Vdwd3::build_structure(std::vector& atom_indices) const diff --git a/source/source_hamilt/module_vdw/vdwd3.h b/source/source_hamilt/module_vdw/vdwd3.h index e8a1e50a9e7..170e6760cef 100644 --- a/source/source_hamilt/module_vdw/vdwd3.h +++ b/source/source_hamilt/module_vdw/vdwd3.h @@ -1,8 +1,8 @@ #ifndef VDWD3_H #define VDWD3_H -#include "vdwd3_types.h" #include "vdw.h" +#include "vdwd3_types.h" #include #include @@ -13,10 +13,7 @@ namespace vdw class Vdwd3 : public Vdw { public: - Vdwd3(const UnitCell& unit_in, - const std::string& xc_name, - const Input_para& input, - std::ofstream* plog = nullptr); + Vdwd3(const UnitCell& unit_in, const std::string& xc_name, const Input_para& input, std::ofstream* plog = nullptr); ~Vdwd3() override = default; diff --git a/source/source_hamilt/module_vdw/vdwd3_data.cpp b/source/source_hamilt/module_vdw/vdwd3_data.cpp index 28a64706d18..8360be19731 100644 --- a/source/source_hamilt/module_vdw/vdwd3_data.cpp +++ b/source/source_hamilt/module_vdw/vdwd3_data.cpp @@ -61,8 +61,7 @@ double reference_c6(int atomic_number_i, int reference_i, int atomic_number_j, i } const std::size_t pair = pair_index(high, low); - const std::size_t offset = kReferenceC6Offset[pair] - + high_reference * reference_count(low) + low_reference; + const std::size_t offset = kReferenceC6Offset[pair] + high_reference * reference_count(low) + low_reference; return kReferenceC6[offset]; } diff --git a/source/source_hamilt/module_vdw/vdwd3_evaluator.cpp b/source/source_hamilt/module_vdw/vdwd3_evaluator.cpp index f130971b37b..de90543175b 100644 --- a/source/source_hamilt/module_vdw/vdwd3_evaluator.cpp +++ b/source/source_hamilt/module_vdw/vdwd3_evaluator.cpp @@ -29,9 +29,7 @@ double dot(const Vec3& lhs, const Vec3& rhs) Vec3 cross(const Vec3& lhs, const Vec3& rhs) { - return Vec3(lhs.y * rhs.z - rhs.y * lhs.z, - lhs.z * rhs.x - rhs.z * lhs.x, - lhs.x * rhs.y - rhs.x * lhs.y); + return Vec3(lhs.y * rhs.z - rhs.y * lhs.z, lhs.z * rhs.x - rhs.z * lhs.x, lhs.x * rhs.y - rhs.x * lhs.y); } double norm(const Vec3& value) @@ -59,7 +57,7 @@ bool validate(const Structure& structure, const Cutoffs& cutoffs, std::string& e error = "atomic-number and position arrays have different sizes"; return false; } - for (int atomic_number : structure.atomic_numbers) + for (int atomic_number: structure.atomic_numbers) { if (atomic_number <= 0 || atomic_number > data::max_element) { @@ -67,21 +65,18 @@ bool validate(const Structure& structure, const Cutoffs& cutoffs, std::string& e return false; } } - if (!std::isfinite(cutoffs.disp2) || !std::isfinite(cutoffs.disp3) - || !std::isfinite(cutoffs.cn) || !(cutoffs.disp2 > 0.0) - || !(cutoffs.disp3 > 0.0) || !(cutoffs.cn > 0.0)) + if (!std::isfinite(cutoffs.disp2) || !std::isfinite(cutoffs.disp3) || !std::isfinite(cutoffs.cn) + || !(cutoffs.disp2 > 0.0) || !(cutoffs.disp3 > 0.0) || !(cutoffs.cn > 0.0)) { error = "DFT-D3 real-space cutoffs must be positive"; return false; } - if (!std::isfinite(cutoffs.width2) || cutoffs.width2 < 0.0 - || cutoffs.width2 > cutoffs.disp2) + if (!std::isfinite(cutoffs.width2) || cutoffs.width2 < 0.0 || cutoffs.width2 > cutoffs.disp2) { error = "two-body smooth width must satisfy 0 <= width <= cutoff"; return false; } - if (!std::isfinite(cutoffs.width3) || cutoffs.width3 < 0.0 - || cutoffs.width3 > cutoffs.disp3) + if (!std::isfinite(cutoffs.width3) || cutoffs.width3 < 0.0 || cutoffs.width3 > cutoffs.disp3) { error = "three-body smooth width must satisfy 0 <= width <= cutoff"; return false; @@ -103,22 +98,19 @@ std::vector lattice_points(const Structure& structure, double cutoff) for (int direction = 0; direction < 3; ++direction) { normal[direction] = normal[direction] / norm(normal[direction]); - repeat[direction] = static_cast( - std::ceil(std::abs(cutoff / dot(normal[direction], structure.lattice[direction])))); + repeat[direction] + = static_cast(std::ceil(std::abs(cutoff / dot(normal[direction], structure.lattice[direction])))); } std::vector translations; - translations.reserve(static_cast((2 * repeat[0] + 1) - * (2 * repeat[1] + 1) - * (2 * repeat[2] + 1))); + translations.reserve(static_cast((2 * repeat[0] + 1) * (2 * repeat[1] + 1) * (2 * repeat[2] + 1))); for (int ix = -repeat[0]; ix <= repeat[0]; ++ix) { for (int iy = -repeat[1]; iy <= repeat[1]; ++iy) { for (int iz = -repeat[2]; iz <= repeat[2]; ++iz) { - translations.push_back(structure.lattice[0] * ix - + structure.lattice[1] * iy + translations.push_back(structure.lattice[0] * ix + structure.lattice[1] * iy + structure.lattice[2] * iz); } } @@ -156,19 +148,16 @@ void smooth_cutoff(double distance, double cutoff, double width, double& value, double coordination_count(int atomic_number_i, int atomic_number_j, double distance) { - const double radius = data::covalent_radius(atomic_number_i) - + data::covalent_radius(atomic_number_j); + const double radius = data::covalent_radius(atomic_number_i) + data::covalent_radius(atomic_number_j); return 1.0 / (1.0 + std::exp(-kCoordinationSteepness * (radius / distance - 1.0))); } double coordination_derivative(int atomic_number_i, int atomic_number_j, double distance) { - const double radius = data::covalent_radius(atomic_number_i) - + data::covalent_radius(atomic_number_j); + const double radius = data::covalent_radius(atomic_number_i) + data::covalent_radius(atomic_number_j); const double exponential = std::exp(-kCoordinationSteepness * (radius / distance - 1.0)); const double denominator = exponential + 1.0; - return -kCoordinationSteepness * radius * exponential - / (distance * distance * denominator * denominator); + return -kCoordinationSteepness * radius * exponential / (distance * distance * denominator * denominator); } std::vector coordination_numbers(const Structure& structure, @@ -182,7 +171,7 @@ std::vector coordination_numbers(const Structure& structure, { for (std::size_t j = 0; j <= i; ++j) { - for (const Vec3& translation : translations) + for (const Vec3& translation: translations) { const Vec3 vector = structure.positions[i] - (structure.positions[j] + translation); const double distance2 = dot(vector, vector); @@ -230,12 +219,14 @@ void reference_weights(const Structure& structure, max_exponent = std::max(max_exponent, raw[reference]); } - // Scale the exponentials to avoid underflow at high coordination numbers. + // Scale the exponentials to avoid underflow at high coordination + // numbers. double weight_sum = 0.0; double weighted_reference_cn = 0.0; for (int reference = 0; reference < references; ++reference) { - // Shifting by the maximum exponent prevents all weights from underflowing. + // Shifting by the maximum exponent prevents all weights from + // underflowing. raw[reference] = std::exp(raw[reference] - max_exponent); weight_sum += raw[reference]; if (derivatives) @@ -255,8 +246,7 @@ void reference_weights(const Structure& structure, if (derivatives) { const double reference_cn = data::reference_cn(atomic_number, reference); - weight_derivatives[index] = 2.0 * kReferenceWeight * weight - * (reference_cn - weighted_reference_cn); + weight_derivatives[index] = 2.0 * kReferenceWeight * weight * (reference_cn - weighted_reference_cn); } } } @@ -290,19 +280,18 @@ void atomic_c6(const Structure& structure, for (int reference_j = 0; reference_j < references_j; ++reference_j) { const double reference_c6 = data::reference_c6(structure.atomic_numbers[i], - reference_i, - structure.atomic_numbers[j], - reference_j); + reference_i, + structure.atomic_numbers[j], + reference_j); const double weight_i = weights[i * data::max_reference + reference_i]; const double weight_j = weights[j * data::max_reference + reference_j]; coefficient += weight_i * weight_j * reference_c6; if (derivatives) { - derivative_i += weight_derivatives[i * data::max_reference + reference_i] - * weight_j * reference_c6; - derivative_j += weight_i - * weight_derivatives[j * data::max_reference + reference_j] - * reference_c6; + derivative_i + += weight_derivatives[i * data::max_reference + reference_i] * weight_j * reference_c6; + derivative_j + += weight_i * weight_derivatives[j * data::max_reference + reference_j] * reference_c6; } } } @@ -356,7 +345,7 @@ void pairwise_dispersion(const Structure& structure, r0 = data::vdw_radius(atomic_number_i, atomic_number_j); } - for (const Vec3& translation : translations) + for (const Vec3& translation: translations) { const Vec3 vector = structure.positions[i] - (structure.positions[j] + translation); const double distance2 = dot(vector, vector); @@ -367,11 +356,7 @@ void pairwise_dispersion(const Structure& structure, const double distance = std::sqrt(distance2); double switching = 1.0; double switching_derivative = 0.0; - smooth_cutoff(distance, - cutoffs.disp2, - cutoffs.width2, - switching, - switching_derivative); + smooth_cutoff(distance, cutoffs.disp2, cutoffs.width2, switching, switching_derivative); double bare_energy = 0.0; double bare_radial = 0.0; @@ -381,12 +366,9 @@ void pairwise_dispersion(const Structure& structure, const double inverse_sixth = 1.0 / (fourth * distance2 + r0_sixth); const double inverse_eighth = 1.0 / (fourth * fourth + r0_eighth); const double derivative_sixth = -6.0 * fourth * inverse_sixth * inverse_sixth; - const double derivative_eighth = -8.0 * fourth * distance2 - * inverse_eighth * inverse_eighth; - bare_energy = parameters.s6 * inverse_sixth - + parameters.s8 * rrij * inverse_eighth; - bare_radial = parameters.s6 * derivative_sixth - + parameters.s8 * rrij * derivative_eighth; + const double derivative_eighth = -8.0 * fourth * distance2 * inverse_eighth * inverse_eighth; + bare_energy = parameters.s6 * inverse_sixth + parameters.s8 * rrij * inverse_eighth; + bare_radial = parameters.s6 * derivative_sixth + parameters.s8 * rrij * derivative_eighth; } else { @@ -398,14 +380,10 @@ void pairwise_dispersion(const Structure& structure, const double t8 = std::pow(parameters.rs8 * r0 / distance, exponent8); const double f6 = 1.0 / (1.0 + 6.0 * t6); const double f8 = 1.0 / (1.0 + 6.0 * t8); - const double derivative6 = -6.0 * f6 / distance2 - + 6.0 * exponent6 * t6 * f6 * f6 / distance2; - const double derivative8 = -8.0 * f8 / distance2 - + 6.0 * exponent8 * t8 * f8 * f8 / distance2; - bare_energy = parameters.s6 * f6 / sixth - + parameters.s8 * rrij * f8 / eighth; - bare_radial = parameters.s6 * derivative6 / sixth - + parameters.s8 * rrij * derivative8 / eighth; + const double derivative6 = -6.0 * f6 / distance2 + 6.0 * exponent6 * t6 * f6 * f6 / distance2; + const double derivative8 = -8.0 * f8 / distance2 + 6.0 * exponent8 * t8 * f8 * f8 / distance2; + bare_energy = parameters.s6 * f6 / sixth + parameters.s8 * rrij * f8 / eighth; + bare_radial = parameters.s6 * derivative6 / sixth + parameters.s8 * rrij * derivative8 / eighth; } const double dispersion = switching * bare_energy; @@ -418,8 +396,7 @@ void pairwise_dispersion(const Structure& structure, if (derivatives) { - const double radial = switching * bare_radial - + switching_derivative * bare_energy / distance; + const double radial = switching * bare_radial + switching_derivative * bare_energy / distance; const Vec3 gradient = vector * (-coefficient * radial); dEdcn[i] -= dc6dcn[i * atoms + j] * dispersion; add_outer(result.virial, gradient, vector, 0.5); @@ -475,7 +452,7 @@ void atm_dispersion(const Structure& structure, const int atomic_number_j = structure.atomic_numbers[j]; const double c6ij = c6[j * atoms + i]; const double r0ij = kAtmRadiusScale * data::vdw_radius(atomic_number_j, atomic_number_i); - for (const Vec3& translation_j : translations) + for (const Vec3& translation_j: translations) { const Vec3 vij = structure.positions[j] + translation_j - structure.positions[i]; const double r2ij = dot(vij, vij); @@ -494,14 +471,12 @@ void atm_dispersion(const Structure& structure, const double c6ik = c6[k * atoms + i]; const double c6jk = c6[k * atoms + j]; const double c9 = -parameters.s9 * std::sqrt(std::abs(c6ij * c6ik * c6jk)); - const double r0ik = kAtmRadiusScale - * data::vdw_radius(atomic_number_k, atomic_number_i); - const double r0jk = kAtmRadiusScale - * data::vdw_radius(atomic_number_k, atomic_number_j); + const double r0ik = kAtmRadiusScale * data::vdw_radius(atomic_number_k, atomic_number_i); + const double r0jk = kAtmRadiusScale * data::vdw_radius(atomic_number_k, atomic_number_j); const double r0 = r0ij * r0ik * r0jk; const double scale = triple_scale(i, j, k); - for (const Vec3& translation_k : translations) + for (const Vec3& translation_k: translations) { const Vec3 vik = structure.positions[k] + translation_k - structure.positions[i]; const double r2ik = dot(vik, vik); @@ -532,10 +507,9 @@ void atm_dispersion(const Structure& structure, const double r5 = r3 * r2; const double damping_power = std::pow(r0 / r1, alpha_third); const double damping = 1.0 / (1.0 + 6.0 * damping_power); - const double angular = 0.375 * (r2ij + r2jk - r2ik) - * (r2ij - r2jk + r2ik) - * (-r2ij + r2jk + r2ik) / r5 - + 1.0 / r3; + const double angular + = 0.375 * (r2ij + r2jk - r2ik) * (r2ij - r2jk + r2ik) * (-r2ij + r2jk + r2ik) / r5 + + 1.0 / r3; const double rr = angular * damping; const double base_energy = rr * c9; const double energy = base_energy * scale * switching; @@ -548,45 +522,39 @@ void atm_dispersion(const Structure& structure, continue; } - const double damping_derivative = -2.0 * alpha * damping_power - * damping * damping; - double angular_derivative = -0.375 - * (r2ij * r2ij * r2ij - + r2ij * r2ij * (r2jk + r2ik) - + r2ij * (3.0 * r2jk * r2jk + 2.0 * r2jk * r2ik - + 3.0 * r2ik * r2ik) - - 5.0 * (r2jk - r2ik) * (r2jk - r2ik) * (r2jk + r2ik)) - / r5; + const double damping_derivative = -2.0 * alpha * damping_power * damping * damping; + double angular_derivative + = -0.375 + * (r2ij * r2ij * r2ij + r2ij * r2ij * (r2jk + r2ik) + + r2ij * (3.0 * r2jk * r2jk + 2.0 * r2jk * r2ik + 3.0 * r2ik * r2ik) + - 5.0 * (r2jk - r2ik) * (r2jk - r2ik) * (r2jk + r2ik)) + / r5; const Vec3 gradient_ij - = vij * (switching * c9 - * (-angular_derivative * damping - + angular * damping_derivative) / r2ij) + = vij + * (switching * c9 * (-angular_derivative * damping + angular * damping_derivative) + / r2ij) - vij * (base_energy * dswij / rij * swik * swjk); angular_derivative = -0.375 - * (r2ik * r2ik * r2ik - + r2ik * r2ik * (r2jk + r2ij) - + r2ik * (3.0 * r2jk * r2jk + 2.0 * r2jk * r2ij - + 3.0 * r2ij * r2ij) - - 5.0 * (r2jk - r2ij) * (r2jk - r2ij) * (r2jk + r2ij)) - / r5; + * (r2ik * r2ik * r2ik + r2ik * r2ik * (r2jk + r2ij) + + r2ik * (3.0 * r2jk * r2jk + 2.0 * r2jk * r2ij + 3.0 * r2ij * r2ij) + - 5.0 * (r2jk - r2ij) * (r2jk - r2ij) * (r2jk + r2ij)) + / r5; const Vec3 gradient_ik - = vik * (switching * c9 - * (-angular_derivative * damping - + angular * damping_derivative) / r2ik) + = vik + * (switching * c9 * (-angular_derivative * damping + angular * damping_derivative) + / r2ik) - vik * (base_energy * dswik / rik * swij * swjk); angular_derivative = -0.375 - * (r2jk * r2jk * r2jk - + r2jk * r2jk * (r2ik + r2ij) - + r2jk * (3.0 * r2ik * r2ik + 2.0 * r2ik * r2ij - + 3.0 * r2ij * r2ij) - - 5.0 * (r2ik - r2ij) * (r2ik - r2ij) * (r2ik + r2ij)) - / r5; + * (r2jk * r2jk * r2jk + r2jk * r2jk * (r2ik + r2ij) + + r2jk * (3.0 * r2ik * r2ik + 2.0 * r2ik * r2ij + 3.0 * r2ij * r2ij) + - 5.0 * (r2ik - r2ij) * (r2ik - r2ij) * (r2ik + r2ij)) + / r5; const Vec3 gradient_jk - = vjk * (switching * c9 - * (-angular_derivative * damping - + angular * damping_derivative) / r2jk) + = vjk + * (switching * c9 * (-angular_derivative * damping + angular * damping_derivative) + / r2jk) - vjk * (base_energy * dswjk / rjk * swij * swik); result.gradient[i] -= (gradient_ij + gradient_ik) * scale; @@ -596,15 +564,9 @@ void atm_dispersion(const Structure& structure, add_outer(result.virial, gradient_ik, vik, scale); add_outer(result.virial, gradient_jk, vjk, scale); - dEdcn[i] -= energy * 0.5 - * (dc6dcn[i * atoms + j] / c6ij - + dc6dcn[i * atoms + k] / c6ik); - dEdcn[j] -= energy * 0.5 - * (dc6dcn[j * atoms + i] / c6ij - + dc6dcn[j * atoms + k] / c6jk); - dEdcn[k] -= energy * 0.5 - * (dc6dcn[k * atoms + i] / c6ik - + dc6dcn[k * atoms + j] / c6jk); + dEdcn[i] -= energy * 0.5 * (dc6dcn[i * atoms + j] / c6ij + dc6dcn[i * atoms + k] / c6ik); + dEdcn[j] -= energy * 0.5 * (dc6dcn[j * atoms + i] / c6ij + dc6dcn[j * atoms + k] / c6jk); + dEdcn[k] -= energy * 0.5 * (dc6dcn[k * atoms + i] / c6ik + dc6dcn[k * atoms + j] / c6jk); } } } @@ -624,7 +586,7 @@ void add_coordination_derivatives(const Structure& structure, { for (std::size_t j = 0; j <= i; ++j) { - for (const Vec3& translation : translations) + for (const Vec3& translation: translations) { const Vec3 vector = structure.positions[i] - (structure.positions[j] + translation); const double distance2 = dot(vector, vector); @@ -634,10 +596,9 @@ void add_coordination_derivatives(const Structure& structure, } const double distance = std::sqrt(distance2); const Vec3 count_derivative - = vector * (coordination_derivative(structure.atomic_numbers[i], - structure.atomic_numbers[j], - distance) - / distance); + = vector + * (coordination_derivative(structure.atomic_numbers[i], structure.atomic_numbers[j], distance) + / distance); const double pair_derivative = dEdcn[i] + dEdcn[j]; result.gradient[i] += count_derivative * pair_derivative; result.gradient[j] -= count_derivative * pair_derivative; @@ -675,8 +636,7 @@ bool evaluate(const Structure& structure, } const std::vector cn_translations = lattice_points(structure, cutoffs.cn); - const std::vector coordination - = coordination_numbers(structure, cn_translations, cutoffs.cn); + const std::vector coordination = coordination_numbers(structure, cn_translations, cutoffs.cn); std::vector weights; std::vector weight_derivatives; diff --git a/source/source_hamilt/module_vdw/vdwd3_parameters.cpp b/source/source_hamilt/module_vdw/vdwd3_parameters.cpp index 5f6a4bffbf1..b215a44d66a 100644 --- a/source/source_hamilt/module_vdw/vdwd3_parameters.cpp +++ b/source/source_hamilt/module_vdw/vdwd3_parameters.cpp @@ -1,4 +1,5 @@ #include "vdwd3_parameters.h" + #include "vdw_xcname.h" #include @@ -50,7 +51,7 @@ std::string normalized(std::string value) value = lower_ascii(value); std::string result; result.reserve(value.size()); - for (unsigned char character : value) + for (unsigned char character: value) { if (std::isalnum(character) || character == '/') { @@ -93,11 +94,16 @@ std::string canonicalize_libxc_pair(const std::string& input) // libxc_component(). Only mixed exchange/correlation names need an // explicit mapping; identical component names are handled above. static const std::map special = { - {"b88:p86", "bp86"}, {"b88:lyp", "blyp"}, - {"b88:pbe", "bpbe"}, {"pber:pbe", "revpbe"}, - {"rpbe:pbe", "rpbe"}, {"optx:lyp", "olyp"}, - {"mpw91:pw91", "mpwpw"}, {"ms2:regtpss", "ms2"}, - {"ms2h:regtpss", "ms2h"}, {"pbe:oppbe", "pbeop"}, + {"b88:p86", "bp86"}, + {"b88:lyp", "blyp"}, + {"b88:pbe", "bpbe"}, + {"pber:pbe", "revpbe"}, + {"rpbe:pbe", "rpbe"}, + {"optx:lyp", "olyp"}, + {"mpw91:pw91", "mpwpw"}, + {"ms2:regtpss", "ms2"}, + {"ms2h:regtpss", "ms2h"}, + {"pbe:oppbe", "pbeop"}, {"b88:opb88", "bop"}, }; const auto found = special.find(exchange + ":" + correlation); @@ -122,8 +128,7 @@ std::string canonicalize_method_name(const std::string& input) if (canonical.empty()) { const std::size_t xc = xc_name.find("_xc_"); - canonical = xc == std::string::npos ? normalized(xc_name) - : normalized(xc_name.substr(xc + 4)); + canonical = xc == std::string::npos ? normalized(xc_name) : normalized(xc_name.substr(xc + 4)); } // ABACUS calls its HSE06 implementation "HSE" in user input and @@ -159,7 +164,7 @@ bool lookup_parameters(const std::string& method, } const char* method_id = nullptr; - for (const MethodAlias& alias : kMethodAliases) + for (const MethodAlias& alias: kMethodAliases) { if (normalized(alias.alias) == lookup_method) { @@ -172,7 +177,7 @@ bool lookup_parameters(const std::string& method, return false; } - for (const DampingParameterRecord& record : kDampingParameters) + for (const DampingParameterRecord& record: kDampingParameters) { if (record.damping == damping && std::string(record.method_id) == method_id) { diff --git a/source/source_hamilt/module_vdw/vdwd3_types.h b/source/source_hamilt/module_vdw/vdwd3_types.h index 803ec8ca69d..67134dcd82b 100644 --- a/source/source_hamilt/module_vdw/vdwd3_types.h +++ b/source/source_hamilt/module_vdw/vdwd3_types.h @@ -16,7 +16,9 @@ struct Vec3 double z = 0.0; Vec3() = default; - Vec3(double x_in, double y_in, double z_in) : x(x_in), y(y_in), z(z_in) {} + Vec3(double x_in, double y_in, double z_in) : x(x_in), y(y_in), z(z_in) + { + } Vec3& operator+=(const Vec3& rhs) { diff --git a/source/source_hamilt/module_vdw/vdwd4.cpp b/source/source_hamilt/module_vdw/vdwd4.cpp index ad4a93e9925..13ed6419337 100644 --- a/source/source_hamilt/module_vdw/vdwd4.cpp +++ b/source/source_hamilt/module_vdw/vdwd4.cpp @@ -1,10 +1,11 @@ #include "vdwd4.h" -#include "vdw_xcname.h" #include "source_base/constants.h" #include "source_base/element_name.h" +#include "vdw_xcname.h" #ifdef __MPI #include "source_base/parallel_reduce.h" + #include #endif #include "source_base/timer.h" @@ -96,13 +97,11 @@ Vdwd4::Vdwd4(const UnitCell& unit_in, const std::string& xc_name, const Input_pa if (smooth_width_2b_ < 0.0 || smooth_width_2b_ > cutoff_disp2_) { - ModuleBase::WARNING_QUIT("Vdwd4::Vdwd4", - "vdw_cutoff_width2 must satisfy 0 <= width <= two-body cutoff"); + ModuleBase::WARNING_QUIT("Vdwd4::Vdwd4", "vdw_cutoff_width2 must satisfy 0 <= width <= two-body cutoff"); } if (smooth_width_3b_ < 0.0 || smooth_width_3b_ > cutoff_disp3_) { - ModuleBase::WARNING_QUIT("Vdwd4::Vdwd4", - "vdw_cutoff_width3 must satisfy 0 <= width <= three-body cutoff"); + ModuleBase::WARNING_QUIT("Vdwd4::Vdwd4", "vdw_cutoff_width3 must satisfy 0 <= width <= three-body cutoff"); } double valence_charge = 0.0; @@ -190,11 +189,9 @@ void Vdwd4::compute(double& energy_ha, build_structure(numbers, positions, lattice, periodic, atom_indices); - if (gradient_ha_bohr != nullptr - && gradient_ha_bohr->size() < static_cast(3 * numbers.size())) + if (gradient_ha_bohr != nullptr && gradient_ha_bohr->size() < static_cast(3 * numbers.size())) { - ModuleBase::WARNING_QUIT("Vdwd4::compute", - "gradient_ha_bohr is too small for the filtered atom set."); + ModuleBase::WARNING_QUIT("Vdwd4::compute", "gradient_ha_bohr is too small for the filtered atom set."); } // These vectors own all arrays passed to DFT-D4. Their data() pointers @@ -248,7 +245,8 @@ void Vdwd4::compute(double& energy_ha, method.push_back('\0'); // Use the DFT-D4 library's internal rational damping parameters. - // The final boolean selects the three-body-specific parameterization in the C API. + // The final boolean selects the three-body-specific parameterization in the + // C API. const bool atm = true; dftd4_param param = dftd4_load_rational_damping(error, method.data(), atm); check_dftd4_error(error, "dftd4_load_rational_damping"); @@ -299,11 +297,11 @@ void Vdwd4::set_force_from_gradient(const std::vector& gradient_ha_bohr, result.has_force = true; } -void Vdwd4::set_stress_from_sigma(const std::array& sigma_ha, - VdwResult& result) const +void Vdwd4::set_stress_from_sigma(const std::array& sigma_ha, VdwResult& result) const { // Tentative mapping consistent with the current D3 convention. - // Confirm sign, transposition and volume normalization by finite-strain tests. + // Confirm sign, transposition and volume normalization by finite-strain + // tests. result.stress = ModuleBase::Matrix3(2.0 * sigma_ha[0], 2.0 * sigma_ha[1], 2.0 * sigma_ha[2], diff --git a/source/source_hamilt/module_vdw/vdwd4.h b/source/source_hamilt/module_vdw/vdwd4.h index 1c1e580e3fe..4c3ec796c5d 100644 --- a/source/source_hamilt/module_vdw/vdwd4.h +++ b/source/source_hamilt/module_vdw/vdwd4.h @@ -20,18 +20,18 @@ class Vdwd4 : public Vdw private: std::string xc_name_; std::string model_name_; - double cutoff_disp2_ = 0.0; // Bohr, two-body dispersion cutoff - double cutoff_disp3_ = 0.0; // Bohr, three-body ATM cutoff - double cutoff_cn_ = 0.0; // Bohr, coordination-number cutoff + double cutoff_disp2_ = 0.0; // Bohr, two-body dispersion cutoff + double cutoff_disp3_ = 0.0; // Bohr, three-body ATM cutoff + double cutoff_cn_ = 0.0; // Bohr, coordination-number cutoff double smooth_width_2b_ = 0.0; // Bohr, two-body cutoff smoothing width double smooth_width_3b_ = 0.0; // Bohr, three-body cutoff smoothing width - double total_charge_ = 0.0; // e, total system charge (sum zv*na - nelec) + double total_charge_ = 0.0; // e, total system charge (sum zv*na - nelec) void evaluate_impl(const VdwRequest& request, VdwResult& result) override; void set_force_from_gradient(const std::vector& gradient_ha_bohr, - const std::vector& atom_indices, - VdwResult& result) const; + const std::vector& atom_indices, + VdwResult& result) const; void set_stress_from_sigma(const std::array& sigma_ha, VdwResult& result) const; void build_structure(std::vector& numbers, From 3e898d108d17c4b38b197c87483b6404f7bb70b3 Mon Sep 17 00:00:00 2001 From: SY Wang Date: Thu, 15 Oct 2026 23:14:36 +0800 Subject: [PATCH 2/5] Add Precommit CI --- .github/workflows/precommit.yml | 20 ++++++++++++++++++++ .github/workflows/test.yml | 11 ----------- 2 files changed, 20 insertions(+), 11 deletions(-) create mode 100644 .github/workflows/precommit.yml diff --git a/.github/workflows/precommit.yml b/.github/workflows/precommit.yml new file mode 100644 index 00000000000..89cf388d635 --- /dev/null +++ b/.github/workflows/precommit.yml @@ -0,0 +1,20 @@ +name: Precommit + +on: + pull_request: + push: + branches: + - develop + +permissions: + contents: read + +jobs: + pre-commit: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v7 + - uses: actions/setup-python@v7 + with: + python-version: '3.12' + - uses: pre-commit/action@v3.0.1 diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 04bad35a5e8..2ae1826cb7a 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -135,17 +135,6 @@ jobs: -DENABLE_PEXSI=ON \ -Werror=dev -# Temporarily removed because no one maintains this now. -# And it will break the CI test workflow. - -# - uses: pre-commit/action@v3.0.1 -# with: -# extra_args: -# --from-ref ${{ github.event.pull_request.base.sha }} -# --to-ref ${{ github.event.pull_request.head.sha }} -# continue-on-error: true -# - uses: pre-commit-ci/lite-action@v1.0.3 - - name: Build run: | cmake --build build -j8 From 928d0b16a315c5ebf57e721d2398275c51c33630 Mon Sep 17 00:00:00 2001 From: SY Wang Date: Thu, 15 Oct 2026 23:07:23 +0800 Subject: [PATCH 3/5] Drop clang-tidy as it tends to cause breakage --- .clang-tidy | 43 ----------------------------------------- .pre-commit-config.yaml | 14 +++++++------- 2 files changed, 7 insertions(+), 50 deletions(-) delete mode 100644 .clang-tidy diff --git a/.clang-tidy b/.clang-tidy deleted file mode 100644 index d5587c46012..00000000000 --- a/.clang-tidy +++ /dev/null @@ -1,43 +0,0 @@ ---- -Checks: ' - modernize-deprecated-headers, - modernize-redundant-void-arg, - modernize-use-bool-literals, - modernize-use-nullptr, - readability-avoid-return-with-void-value, - readability-braces-around-statements - - # llvm-header-guard, - # bugprone-*, - # clang-analyzer-*, - # clang-diagnostic-*, - # cppcoreguidelines-avoid-non-const-global-variables, - # cppcoreguidelines-macro-usage, - # cppcoreguidelines-pro-type-member-init, - # fuchsia-multiple-inheritance, - # google-*, - # -google-runtime-references, - # misc-*, - # modernize-shrink-to-fit, - # modernize-use-starts-ends-with, - # mpi-*, - # performance-*, - # -performance-avoid-endl, - # readability-*, - # -readability-else-after-return, - # -readability-magic-numbers, - # -readability-isolate-declaration, - # -readability-braces-around-statements, - # -readability-implicit-bool-conversion, - # -google-readability-braces-around-statements, - # -misc-unused-parameters, - # -modernize-pass-by-value, - # -modernize-use-auto, - # -google-readability-casting, - # -readability-make-member-function-const, - # -cppcoreguidelines-avoid-non-const-global-variables, - # -readability-function-cognitive-complexity, - # -google-build-using-namespace, - # -bugprone-narrowing-conversions, -' -... diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 5a76d574d54..ab7a8094aac 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -6,15 +6,15 @@ repos: - id: mixed-line-ending args: [--fix=lf] exclude: '(\.bat|\.cmd)$' - - repo: https://github.com/pocc/pre-commit-hooks - rev: v1.3.5 + - repo: https://github.com/pre-commit/mirrors-clang-format + rev: v23.1.3 hooks: - id: clang-format - args: [-i, --no-diff] - # - id: clang-tidy - # args: [-p=build, --fix-errors] - # - id: cppcheck - # - id: cpplint + # - repo: https://github.com/cheshirekow/cmake-format-precommit + # rev: v0.6.13 + # hooks: + # - id: cmake-format + # additional_dependencies: [pyyaml] - repo: local hooks: - id: abacus-agent-governance From 496ff0a40f2ca5a2827a8990357bbeb08c535c2b Mon Sep 17 00:00:00 2001 From: SY Wang Date: Mon, 19 Oct 2026 01:10:43 +0800 Subject: [PATCH 4/5] Format two cpp files under tests directory --- tests/09_DeePKS/check_file.cpp | 51 +++++++++++++++++----------------- tests/integrate/check_file.cpp | 51 +++++++++++++++++----------------- 2 files changed, 50 insertions(+), 52 deletions(-) diff --git a/tests/09_DeePKS/check_file.cpp b/tests/09_DeePKS/check_file.cpp index d7e1391c9e0..049d4504da9 100644 --- a/tests/09_DeePKS/check_file.cpp +++ b/tests/09_DeePKS/check_file.cpp @@ -1,13 +1,13 @@ -#include -#include -#include +#include #include #include -#include +#include +#include +#include using namespace std; -int main(int argc, char **argv) +int main(int argc, char** argv) { string ref_file = argv[1]; string cal_file = argv[2]; @@ -20,21 +20,21 @@ int main(int argc, char **argv) ifstream ref(ref_file.c_str(), ios::in); ifstream cal(cal_file.c_str(), ios::in); - if (!ref) - { + if (!ref) + { cout << "\e[1;31m [ FAILED ] \e[0m"; - cout << "Can't find " << ref_file << " !" << endl; - return 0; - } + cout << "Can't find " << ref_file << " !" << endl; + return 0; + } - if (!cal) - { + if (!cal) + { cout << "\e[1;31m [ FAILED ] \e[0m"; - cout << "Can't find " << cal_file << " !" << endl; - return 0; - } + cout << "Can't find " << cal_file << " !" << endl; + return 0; + } - int num_ref=0, num_cal=0; + int num_ref = 0, num_cal = 0; string word_ref, word_cal; while (ref.good()) @@ -49,7 +49,7 @@ int main(int argc, char **argv) num_cal++; } - if(num_ref != num_cal) + if (num_ref != num_cal) { cout << "\e[1;31m [ FAILED ] \e[0m"; cout << "The num of data is different !" << endl; @@ -66,32 +66,31 @@ int main(int argc, char **argv) double data_ref, data_cal; - for(int index=0; index> word_ref; cal >> word_cal; - if(word_cal != word_ref) + if (word_cal != word_ref) { istringstream string_cal(word_cal); istringstream string_ref(word_ref); string_cal >> data_cal; string_ref >> data_ref; - if(data_ref==0 || data_cal==0) + if (data_ref == 0 || data_cal == 0) { cout << "\e[1;31m [ FAILED ] \e[0m"; - cout << "Data No." << index+1 << " cal=" << word_cal - << " ref=" << word_ref << endl; + cout << "Data No." << index + 1 << " cal=" << word_cal << " ref=" << word_ref << endl; } - else if(fabs(data_ref-data_cal) > threshold) + else if (fabs(data_ref - data_cal) > threshold) { cout << "\e[1;31m [ FAILED ] \e[0m"; - cout << "Data No." << index+1 << " cal=" << data_cal - << " ref=" << data_ref << " deviation=" << fabs(data_ref-data_cal) << endl; + cout << "Data No." << index + 1 << " cal=" << data_cal << " ref=" << data_ref + << " deviation=" << fabs(data_ref - data_cal) << endl; return 0; } } } - + return 1; } \ No newline at end of file diff --git a/tests/integrate/check_file.cpp b/tests/integrate/check_file.cpp index d7e1391c9e0..049d4504da9 100644 --- a/tests/integrate/check_file.cpp +++ b/tests/integrate/check_file.cpp @@ -1,13 +1,13 @@ -#include -#include -#include +#include #include #include -#include +#include +#include +#include using namespace std; -int main(int argc, char **argv) +int main(int argc, char** argv) { string ref_file = argv[1]; string cal_file = argv[2]; @@ -20,21 +20,21 @@ int main(int argc, char **argv) ifstream ref(ref_file.c_str(), ios::in); ifstream cal(cal_file.c_str(), ios::in); - if (!ref) - { + if (!ref) + { cout << "\e[1;31m [ FAILED ] \e[0m"; - cout << "Can't find " << ref_file << " !" << endl; - return 0; - } + cout << "Can't find " << ref_file << " !" << endl; + return 0; + } - if (!cal) - { + if (!cal) + { cout << "\e[1;31m [ FAILED ] \e[0m"; - cout << "Can't find " << cal_file << " !" << endl; - return 0; - } + cout << "Can't find " << cal_file << " !" << endl; + return 0; + } - int num_ref=0, num_cal=0; + int num_ref = 0, num_cal = 0; string word_ref, word_cal; while (ref.good()) @@ -49,7 +49,7 @@ int main(int argc, char **argv) num_cal++; } - if(num_ref != num_cal) + if (num_ref != num_cal) { cout << "\e[1;31m [ FAILED ] \e[0m"; cout << "The num of data is different !" << endl; @@ -66,32 +66,31 @@ int main(int argc, char **argv) double data_ref, data_cal; - for(int index=0; index> word_ref; cal >> word_cal; - if(word_cal != word_ref) + if (word_cal != word_ref) { istringstream string_cal(word_cal); istringstream string_ref(word_ref); string_cal >> data_cal; string_ref >> data_ref; - if(data_ref==0 || data_cal==0) + if (data_ref == 0 || data_cal == 0) { cout << "\e[1;31m [ FAILED ] \e[0m"; - cout << "Data No." << index+1 << " cal=" << word_cal - << " ref=" << word_ref << endl; + cout << "Data No." << index + 1 << " cal=" << word_cal << " ref=" << word_ref << endl; } - else if(fabs(data_ref-data_cal) > threshold) + else if (fabs(data_ref - data_cal) > threshold) { cout << "\e[1;31m [ FAILED ] \e[0m"; - cout << "Data No." << index+1 << " cal=" << data_cal - << " ref=" << data_ref << " deviation=" << fabs(data_ref-data_cal) << endl; + cout << "Data No." << index + 1 << " cal=" << data_cal << " ref=" << data_ref + << " deviation=" << fabs(data_ref - data_cal) << endl; return 0; } } } - + return 1; } \ No newline at end of file From 64426afecd7a45c083ed83498b94a8aad76dae82 Mon Sep 17 00:00:00 2001 From: SY Wang Date: Sat, 17 Oct 2026 00:56:33 +0800 Subject: [PATCH 5/5] add "continue-on-error" to precommit workflow for now --- .github/workflows/precommit.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/precommit.yml b/.github/workflows/precommit.yml index 89cf388d635..a4a9764dd12 100644 --- a/.github/workflows/precommit.yml +++ b/.github/workflows/precommit.yml @@ -18,3 +18,5 @@ jobs: with: python-version: '3.12' - uses: pre-commit/action@v3.0.1 + # TODO: Drop "continue-on-error" + continue-on-error: true