diff --git a/.clang-format b/.clang-format index 8970df99c29..f16ef75761c 100644 --- a/.clang-format +++ b/.clang-format @@ -1,6 +1,6 @@ --- Language: Cpp -BasedOnStyle: Microsoft +BasedOnStyle: LLVM AlwaysBreakTemplateDeclarations: Yes @@ -25,6 +25,8 @@ IncludeIsMainSourceRegex: '' IndentWidth: 4 IndentWrappedFunctionNames: true +ColumnLimit: 120 + # About Point DerivePointerAlignment: false PointerAlignment: Left 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/.github/workflows/precommit.yml b/.github/workflows/precommit.yml new file mode 100644 index 00000000000..a4a9764dd12 --- /dev/null +++ b/.github/workflows/precommit.yml @@ -0,0 +1,22 @@ +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 + # TODO: Drop "continue-on-error" + continue-on-error: true diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 87bf922499b..54634c7388a 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -99,17 +99,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 diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 8873791f2ce..ab7a8094aac 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -6,18 +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] - - id: clang-tidy - args: [-p=build, --fix-errors] - # - id: oclint - # - id: uncrustify - # - id: cppcheck - # - id: cpplint - # - id: include-what-you-use + - id: clang-format + # - 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 diff --git a/source/source_hamilt/module_vdw/test/vdw_test.cpp b/source/source_hamilt/module_vdw/test/vdw_test.cpp index f03f4fb1588..a88dec51d8b 100644 --- a/source/source_hamilt/module_vdw/test/vdw_test.cpp +++ b/source/source_hamilt/module_vdw/test/vdw_test.cpp @@ -1,104 +1,85 @@ -#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 -*/ - -pseudo::pseudo() -{ -} -pseudo::~pseudo() -{ -} -Atom::Atom() -{ -} -Atom::~Atom() -{ -} -Atom_pseudo::Atom_pseudo() -{ -} -Atom_pseudo::~Atom_pseudo() -{ -} -UnitCell::UnitCell() -{ -} -UnitCell::~UnitCell() -{ -} -Magnetism::Magnetism() -{ -} -Magnetism::~Magnetism() -{ -} - -SepPot::SepPot(){} -SepPot::~SepPot(){} + * - 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() {} +pseudo::~pseudo() {} +Atom::Atom() {} +Atom::~Atom() {} +Atom_pseudo::Atom_pseudo() {} +Atom_pseudo::~Atom_pseudo() {} +UnitCell::UnitCell() {} +UnitCell::~UnitCell() {} +Magnetism::Magnetism() {} +Magnetism::~Magnetism() {} + +SepPot::SepPot() {} +SepPot::~SepPot() {} Sep_Cell::Sep_Cell() noexcept {} Sep_Cell::~Sep_Cell() noexcept {} -struct atomtype_ -{ +struct atomtype_ { std::string atomname; std::vector> coordinate; }; -struct stru_ -{ +struct stru_ { std::vector cell; 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]; ucell.nat = 0; int nmax = 0; - for (int i = 0; i < coord.size(); i++) - { + for (int i = 0; i < coord.size(); i++) { ucell.atoms[i].label = coord[i].atomname; ucell.atoms[i].ncpp.psd = coord[i].atomname; ucell.atoms[i].na = coord[i].coordinate.size(); ucell.atoms[i].tau.resize(ucell.atoms[i].na); ucell.atoms[i].taud.resize(ucell.atoms[i].na); - for (int j = 0; j < ucell.atoms[i].na; j++) - { + for (int j = 0; j < ucell.atoms[i].na; j++) { std::vector this_atom = coord[i].coordinate[j]; ucell.atoms[i].tau[j] = ModuleBase::Vector3(this_atom[0], this_atom[1], this_atom[2]); ModuleBase::Mathzone::Cartesian_to_Direct(ucell.atoms[i].tau[j].x, @@ -118,38 +99,29 @@ void construct_ucell(stru_ &stru, UnitCell &ucell) ucell.atoms[i].taud[j].z); } ucell.nat += ucell.atoms[i].na; - if (ucell.atoms[i].na > nmax) - { + if (ucell.atoms[i].na > nmax) { nmax = ucell.atoms[i].na; } } 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}}}}}; +void construct_two_si_with_ghost(UnitCell& ucell) { + 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,22 +129,16 @@ 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); EXPECT_EQ(ghost.has_stress, reference.has_stress); - if (reference.has_force && ghost.has_force) - { + if (reference.has_force && ghost.has_force) { EXPECT_EQ(ghost.force.size(), static_cast(reference_nat + 1)); if (ghost.force.size() == static_cast(reference_nat + 1) - && reference.force.size() == static_cast(reference_nat)) - { - for (int iat = 0; iat < reference_nat; ++iat) - { + && reference.force.size() == static_cast(reference_nat)) { + for (int iat = 0; iat < reference_nat; ++iat) { EXPECT_NEAR(ghost.force[iat].x, reference.force[iat].x, 1E-12); EXPECT_NEAR(ghost.force[iat].y, reference.force[iat].y, 1E-12); EXPECT_NEAR(ghost.force[iat].z, reference.force[iat].z, 1E-12); @@ -183,8 +149,7 @@ void expect_empty_atoms_excluded(const vdw::VdwResult& reference, } } - if (reference.has_stress && ghost.has_stress) - { + if (reference.has_stress && ghost.has_stress) { EXPECT_NEAR(ghost.stress.e11, reference.stress.e11, 1E-12); EXPECT_NEAR(ghost.stress.e12, reference.stress.e12, 1E-12); EXPECT_NEAR(ghost.stress.e13, reference.stress.e13, 1E-12); @@ -197,24 +162,20 @@ void expect_empty_atoms_excluded(const vdw::VdwResult& reference, } } -class vdwd2Test: public testing::Test -{ - protected: +class vdwd2Test : public testing::Test { + 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,13 +185,10 @@ class vdwd2Test: public testing::Test input.vdw_radius_unit = "Bohr"; } - void TearDown(){ - ClearUcell(ucell); - } + void TearDown() { ClearUcell(ucell); } }; -TEST_F(vdwd2Test, D2Default) -{ +TEST_F(vdwd2Test, D2Default) { vdw::Vdwd2 vdwd2_test(ucell); vdwd2_test.parameter().initial_parameters(input); vdwd2_test.parameter().initset(ucell); @@ -238,59 +196,59 @@ 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); } -TEST_F(vdwd2Test, WrongVdwType) -{ +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) -{ +TEST_F(vdwd2Test, D2ReadFile) { input.vdw_C6_file = "c6.txt"; input.vdw_R0_file = "r0.txt"; 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) -{ +TEST_F(vdwd2Test, D2ReadFileError) { input.vdw_C6_file = "c6_wrong.txt"; input.vdw_R0_file = "r0_wrong.txt"; 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(); } -TEST_F(vdwd2Test, D2c6UniteVA6) -{ +TEST_F(vdwd2Test, D2c6UniteVA6) { input.vdw_C6_unit = "eVA6"; vdw::Vdwd2 vdwd2_test(ucell); 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) -{ +TEST_F(vdwd2Test, D2r0UnitBohr) { input.vdw_R0_unit = "Bohr"; vdw::Vdwd2 vdwd2_test(ucell); @@ -298,40 +256,40 @@ TEST_F(vdwd2Test, D2r0UnitBohr) EXPECT_EQ(vdwd2_test.parameter().R0().at("Si"), 1.716); } -TEST_F(vdwd2Test, D2WrongUnit) -{ +TEST_F(vdwd2Test, D2WrongUnit) { input.vdw_R0_unit = "B"; input.vdw_C6_unit = "eV"; 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(); } -TEST_F(vdwd2Test, D2RadiusUnitAngstrom) -{ +TEST_F(vdwd2Test, D2RadiusUnitAngstrom) { input.vdw_cutoff_radius = "56.6918"; input.vdw_radius_unit = "Angstrom"; 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) -{ +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); EXPECT_EQ(vdwd2_test.parameter().period(), input.vdw_cutoff_period); } -TEST_F(vdwd2Test, D2R0ZeroQuit) -{ +TEST_F(vdwd2Test, D2R0ZeroQuit) { vdw::Vdwd2 vdwd2_test(ucell); vdwd2_test.parameter().initial_parameters(input); // r0_zero.txt sets R0 for Si to zero through the public reader, so that @@ -343,79 +301,67 @@ TEST_F(vdwd2Test, D2R0ZeroQuit) std::string output = testing::internal::GetCapturedStdout(); } -TEST_F(vdwd2Test, D2GetEnergy) -{ +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) -{ +TEST_F(vdwd2Test, D2GetForce) { auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); EXPECT_NEAR(result.energy, -0.034526673470525196, 1E-10); 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) -{ +TEST_F(vdwd2Test, D2GetStress) { auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); EXPECT_NEAR(result.energy, -0.034526673470525196, 1E-10); 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); -} - -TEST_F(vdwd2Test, EmptyAtomsAreExcluded) -{ - const vdw::VdwResult reference = - vdw::make_vdw(ucell, input)->evaluate(vdw::VdwRequest(true, true)); + 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)); 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 -{ - protected: +class vdwd3Test : public testing::Test { + 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,13 +373,10 @@ class vdwd3Test: public testing::Test input.vdw_cn_thr = 40; } - void TearDown(){ - ClearUcell(ucell); - } + void TearDown() { ClearUcell(ucell); } }; -TEST_F(vdwd3Test, AutomaticXcInferenceFromPseudopotential) -{ +TEST_F(vdwd3Test, AutomaticXcInferenceFromPseudopotential) { input.dft_functional = "default"; ucell.atoms[0].ncpp.xc_func = "GGA_X_PBE+GGA_C_PBE"; @@ -445,8 +388,7 @@ TEST_F(vdwd3Test, AutomaticXcInferenceFromPseudopotential) EXPECT_NEAR(result.energy, -0.022043153033290883, 1E-10); } -TEST_F(vdwd3Test, FullyCustomParametersDoNotRequireKnownFunctional) -{ +TEST_F(vdwd3Test, FullyCustomParametersDoNotRequireKnownFunctional) { input.dft_functional = "not-a-known-d3-functional"; input.vdw_s6 = "1.0"; input.vdw_s8 = "0.722"; @@ -458,16 +400,14 @@ TEST_F(vdwd3Test, FullyCustomParametersDoNotRequireKnownFunctional) EXPECT_NEAR(result.energy, -0.022043153033290883, 1E-10); } -TEST_F(vdwd3Test, D30GetEnergy) -{ +TEST_F(vdwd3Test, D30GetEnergy) { 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.022043153033290883, 1E-10); } -TEST_F(vdwd3Test, D30LibxcFunctionalName) -{ +TEST_F(vdwd3Test, D30LibxcFunctionalName) { input.dft_functional = "XC_GGA_X_PBE+XC_GGA_C_PBE"; auto vdw_solver = vdw::make_vdw(ucell, input); @@ -475,58 +415,52 @@ TEST_F(vdwd3Test, D30LibxcFunctionalName) EXPECT_NEAR(result.energy, -0.022043153033290883, 1E-10); } -TEST_F(vdwd3Test, D30GetForce) -{ +TEST_F(vdwd3Test, D30GetForce) { auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); EXPECT_NEAR(result.energy, -0.022043153033290883, 1E-10); 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) -{ +TEST_F(vdwd3Test, D30GetStress) { auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); EXPECT_NEAR(result.energy, -0.022043153033290883, 1E-10); 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); -} - -TEST_F(vdwd3Test, EmptyAtomsAreExcluded) -{ - const vdw::VdwResult reference = - vdw::make_vdw(ucell, input)->evaluate(vdw::VdwRequest(true, true)); + 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)); 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); } -TEST_F(vdwd3Test, D3bjGetEnergy) -{ +TEST_F(vdwd3Test, D3bjGetEnergy) { input.vdw_method = "d3_bj"; auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false)); @@ -534,8 +468,7 @@ TEST_F(vdwd3Test, D3bjGetEnergy) EXPECT_NEAR(ene, -0.047425367813039881, 1E-10); } -TEST_F(vdwd3Test, D3bjGetForce) -{ +TEST_F(vdwd3Test, D3bjGetForce) { input.vdw_method = "d3_bj"; auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); @@ -543,16 +476,15 @@ 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) -{ +TEST_F(vdwd3Test, D3bjGetStress) { input.vdw_method = "d3_bj"; auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); @@ -560,30 +492,28 @@ 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 -{ - protected: +class vdwd3abcTest : public testing::Test { + 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,68 +525,58 @@ class vdwd3abcTest: public testing::Test input.vdw_cn_thr = 40; } - void TearDown(){ - ClearUcell(ucell); - } + void TearDown() { ClearUcell(ucell); } }; - -TEST_F(vdwd3abcTest, InconsistentPseudopotentialXcIsRejected) -{ +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) -{ +TEST_F(vdwd3abcTest, D30GetEnergy) { 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.015864127638792386, 1E-10); } -TEST_F(vdwd3abcTest, D30GetForce) -{ +TEST_F(vdwd3abcTest, D30GetForce) { auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); EXPECT_NEAR(result.energy, -0.015864127638792386, 1E-10); 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) -{ +TEST_F(vdwd3abcTest, D30GetStress) { auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); EXPECT_NEAR(result.energy, -0.015864127638792386, 1E-10); 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); -} - -TEST_F(vdwd3abcTest, D3bjGetEnergy) -{ + 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) { input.vdw_method = "d3_bj"; auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false)); @@ -664,8 +584,7 @@ TEST_F(vdwd3abcTest, D3bjGetEnergy) EXPECT_NEAR(ene, -0.030643056581152218, 1E-10); } -TEST_F(vdwd3abcTest, D3bjGetForce) -{ +TEST_F(vdwd3abcTest, D3bjGetForce) { input.vdw_method = "d3_bj"; auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); @@ -673,16 +592,15 @@ 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) -{ +TEST_F(vdwd3abcTest, D3bjGetStress) { input.vdw_method = "d3_bj"; auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); @@ -690,25 +608,22 @@ 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); -} - -TEST_F(vdwd3abcTest, EmptyAtomsAreExcludedWithThreeBodyTerm) -{ + 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 +635,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 +643,16 @@ TEST_F(vdwd3abcTest, EmptyAtomsAreExcludedWithThreeBodyTerm) #ifdef __DFTD4 -class vdwd4Test: public testing::Test -{ - protected: +class vdwd4Test : public testing::Test { + 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,21 +666,17 @@ class vdwd4Test: public testing::Test input.vdw_cn_thr = 30; } - void TearDown(){ - ClearUcell(ucell); - } + void TearDown() { ClearUcell(ucell); } }; -TEST_F(vdwd4Test, D4GetEnergy) -{ +TEST_F(vdwd4Test, D4GetEnergy) { 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.049988366372568399, 1E-10); } -TEST_F(vdwd4Test, D4LibxcFunctionalName) -{ +TEST_F(vdwd4Test, D4LibxcFunctionalName) { input.vdw_d4_xc = "default"; input.dft_functional = "XC_GGA_X_PBE+XC_GGA_C_PBE"; @@ -779,8 +685,7 @@ TEST_F(vdwd4Test, D4LibxcFunctionalName) EXPECT_NEAR(result.energy, -0.049988366372568399, 1E-10); } -TEST_F(vdwd4Test, D4GetEnergyForChargedSystem) -{ +TEST_F(vdwd4Test, D4GetEnergyForChargedSystem) { input.nelec = 7.0; auto vdw_solver = vdw::make_vdw(ucell, input); @@ -789,8 +694,7 @@ TEST_F(vdwd4Test, D4GetEnergyForChargedSystem) EXPECT_NEAR(ene, -0.043594505741178208, 1E-10); } -TEST_F(vdwd4Test, D4GetForce) -{ +TEST_F(vdwd4Test, D4GetForce) { auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); EXPECT_NEAR(result.energy, -0.049988366372568399, 1E-10); @@ -805,8 +709,7 @@ TEST_F(vdwd4Test, D4GetForce) EXPECT_NEAR(force[1].z, 0.0, 1e-12); } -TEST_F(vdwd4Test, D4GetStress) -{ +TEST_F(vdwd4Test, D4GetStress) { auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); EXPECT_NEAR(result.energy, -0.049988366372568399, 1E-10); @@ -824,22 +727,18 @@ TEST_F(vdwd4Test, D4GetStress) EXPECT_NEAR(stress.e33, 0.00016713814230248975, 1e-12); } -TEST_F(vdwd4Test, EmptyAtomsAreExcluded) -{ - const vdw::VdwResult reference = - vdw::make_vdw(ucell, input)->evaluate(vdw::VdwRequest(true, true)); +TEST_F(vdwd4Test, EmptyAtomsAreExcluded) { + 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); } -TEST_F(vdwd4Test, D4SGetEnergy) -{ +TEST_F(vdwd4Test, D4SGetEnergy) { input.vdw_d4_model = "d4s"; auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false)); @@ -847,8 +746,7 @@ TEST_F(vdwd4Test, D4SGetEnergy) EXPECT_NEAR(ene, -0.056385156117091439, 1E-10); } -TEST_F(vdwd4Test, D4SGetForce) -{ +TEST_F(vdwd4Test, D4SGetForce) { input.vdw_d4_model = "d4s"; auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); @@ -864,8 +762,7 @@ TEST_F(vdwd4Test, D4SGetForce) EXPECT_NEAR(force[1].z, 0.0, 1e-12); } -TEST_F(vdwd4Test, D4SGetStress) -{ +TEST_F(vdwd4Test, D4SGetStress) { input.vdw_d4_model = "d4s"; auto vdw_solver = vdw::make_vdw(ucell, input); const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); @@ -884,16 +781,13 @@ TEST_F(vdwd4Test, D4SGetStress) EXPECT_NEAR(stress.e33, 0.0001579183447576555, 1e-12); } -TEST_F(vdwd4Test, EmptyAtomsAreExcludedForD4S) -{ +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,8 +795,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); int result = RUN_ALL_TESTS(); 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..616d5d06fb5 100644 --- a/source/source_hamilt/module_vdw/test/vdwd3_evaluator_test.cpp +++ b/source/source_hamilt/module_vdw/test/vdwd3_evaluator_test.cpp @@ -1,16 +1,15 @@ -#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 -namespace -{ +namespace { using vdw::d3::Cutoffs; using vdw::d3::Damping; @@ -19,42 +18,34 @@ using vdw::d3::Result; using vdw::d3::Structure; using vdw::d3::Vec3; -double& coordinate(Vec3& vector, int component) -{ - if (component == 0) - { +double& coordinate(Vec3& vector, int component) { + if (component == 0) { return vector.x; } - if (component == 1) - { + if (component == 1) { return vector.y; } return vector.z; } -double coordinate(const Vec3& vector, int component) -{ - if (component == 0) - { +double coordinate(const Vec3& vector, int component) { + if (component == 0) { return vector.x; } - if (component == 1) - { + if (component == 1) { return vector.y; } return vector.z; } -double energy(const Structure& structure, const Parameters& parameters, const Cutoffs& cutoffs) -{ +double energy(const Structure& structure, const Parameters& parameters, const Cutoffs& cutoffs) { Result result; std::string error; EXPECT_TRUE(vdw::d3::evaluate(structure, parameters, cutoffs, false, result, error)) << error; return result.energy; } -Structure mindless01() -{ +Structure mindless01() { Structure structure; structure.atomic_numbers = {11, 1, 8, 1, 9, 1, 1, 8, 7, 1, 1, 17, 5, 5, 7, 13}; structure.positions = { @@ -78,8 +69,7 @@ Structure mindless01() return structure; } -Structure mindless17() -{ +Structure mindless17() { Structure structure; structure.atomic_numbers = {12, 5, 1, 6, 14, 5, 9, 9, 1, 1, 9, 1, 1, 8, 1, 1}; structure.positions = { @@ -103,8 +93,7 @@ Structure mindless17() return structure; } -Structure mindless09() -{ +Structure mindless09() { Structure structure; structure.atomic_numbers = {1, 1, 1, 1, 3, 1, 6, 5, 1, 1, 14, 1, 17, 9, 1, 5}; structure.positions = { @@ -128,11 +117,9 @@ Structure mindless09() return structure; } -Structure actinides() -{ +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}, @@ -155,12 +142,10 @@ Structure actinides() return structure; } -Structure x23_acetic() -{ +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,29 +181,22 @@ 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; } } // namespace -TEST(VdwXcName, NormalizesAbacusAndLibxcSpellings) -{ +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) -{ +TEST(D3Parameters, PbeAndLibxcAliases) { Parameters zero; std::string canonical; ASSERT_TRUE(vdw::d3::lookup_parameters("pbe", Damping::Zero, zero, canonical)); @@ -229,10 +207,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,93 +217,62 @@ 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) -{ +TEST(D3Data, ReferencePackingAndSymmetry) { EXPECT_EQ(vdw::d3::data::max_element, 103); EXPECT_EQ(vdw::d3::data::max_reference, 7); int maximum_reference_count = 0; - for (int atomic_number_i = 1; atomic_number_i <= vdw::d3::data::max_element; ++atomic_number_i) - { + for (int atomic_number_i = 1; atomic_number_i <= vdw::d3::data::max_element; ++atomic_number_i) { const int references_i = vdw::d3::data::reference_count(atomic_number_i); ASSERT_GE(references_i, 1); ASSERT_LE(references_i, vdw::d3::data::max_reference); maximum_reference_count = std::max(maximum_reference_count, references_i); EXPECT_TRUE(std::isfinite(vdw::d3::data::covalent_radius(atomic_number_i))); EXPECT_TRUE(std::isfinite(vdw::d3::data::r4r2(atomic_number_i))); - for (int reference_i = 0; reference_i < references_i; ++reference_i) - { + for (int reference_i = 0; reference_i < references_i; ++reference_i) { EXPECT_TRUE(std::isfinite(vdw::d3::data::reference_cn(atomic_number_i, reference_i))); } - for (int atomic_number_j = 1; atomic_number_j <= atomic_number_i; ++atomic_number_j) - { + for (int atomic_number_j = 1; atomic_number_j <= atomic_number_i; ++atomic_number_j) { const int references_j = vdw::d3::data::reference_count(atomic_number_j); EXPECT_DOUBLE_EQ(vdw::d3::data::vdw_radius(atomic_number_i, atomic_number_j), vdw::d3::data::vdw_radius(atomic_number_j, atomic_number_i)); - for (int reference_i = 0; reference_i < references_i; ++reference_i) - { - 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); + for (int reference_i = 0; reference_i < references_i; ++reference_i) { + 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); 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)); } } } @@ -337,32 +281,27 @@ TEST(D3Data, ReferencePackingAndSymmetry) EXPECT_EQ(vdw::d3::data::reference_count(103), 7); } -TEST(D3Evaluator, MatchesSdftd3Mindless01) -{ +TEST(D3Evaluator, MatchesSdftd3Mindless01) { Parameters parameters; std::string canonical; ASSERT_TRUE(vdw::d3::lookup_parameters("pbe", Damping::Rational, parameters, canonical)); 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) -{ +TEST(D3Evaluator, MatchesSdftd3Actinides) { Parameters parameters; std::string canonical; ASSERT_TRUE(vdw::d3::lookup_parameters("pbe", Damping::Rational, parameters, canonical)); 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) -{ +TEST(D3Evaluator, MatchesSdftd3PeriodicX23Acetic) { Parameters parameters; std::string canonical; ASSERT_TRUE(vdw::d3::lookup_parameters("pbe", Damping::Rational, parameters, canonical)); @@ -374,46 +313,36 @@ 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) -{ +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) -{ +TEST(D3Evaluator, MatchesSdftd3ZeroDampingMindless09) { Parameters parameters; std::string canonical; ASSERT_TRUE(vdw::d3::lookup_parameters("rpbe", Damping::Zero, parameters, canonical)); 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). +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). 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,51 +361,43 @@ 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) - { - const double magnitude = std::sqrt(gradient.x * gradient.x - + gradient.y * gradient.y - + gradient.z * gradient.z); + for (const Vec3& gradient: analytic.gradient) { + 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); } const double step = 1.0e-5; - for (std::size_t atom = 0; atom < structure.positions.size(); ++atom) - { - for (int component = 0; component < 3; ++component) - { + for (std::size_t atom = 0; atom < structure.positions.size(); ++atom) { + for (int component = 0; component < 3; ++component) { Structure plus = structure; 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); } } Structure plus = structure; Structure minus = structure; - for (std::size_t atom = 0; atom < structure.positions.size(); ++atom) - { + for (std::size_t atom = 0; atom < structure.positions.size(); ++atom) { const double displacement = step * structure.positions[atom].x; plus.positions[atom].x += displacement; minus.positions[atom].x -= displacement; } - for (int vector = 0; vector < 3; ++vector) - { + for (int vector = 0; vector < 3; ++vector) { const double displacement = step * structure.lattice[vector].x; 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); } -TEST(D3Evaluator, SmoothCutoffGradientAndVirial) -{ +TEST(D3Evaluator, SmoothCutoffGradientAndVirial) { Structure structure; structure.atomic_numbers = {6, 8, 7}; structure.positions = {{0.0, 0.0, 0.0}, {5.5, 0.4, 0.2}, {2.7, 4.8, -0.3}}; @@ -499,50 +420,40 @@ TEST(D3Evaluator, SmoothCutoffGradientAndVirial) const double step = 1.0e-5; double maximum_gradient_error = 0.0; - for (std::size_t atom = 0; atom < structure.positions.size(); ++atom) - { - for (int component = 0; component < 3; ++component) - { + for (std::size_t atom = 0; atom < structure.positions.size(); ++atom) { + for (int component = 0; component < 3; ++component) { Structure plus = structure; 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); double maximum_virial_error = 0.0; - for (int strain_row = 0; strain_row < 3; ++strain_row) - { - for (int strain_column = 0; strain_column < 3; ++strain_column) - { + for (int strain_row = 0; strain_row < 3; ++strain_row) { + for (int strain_column = 0; strain_column < 3; ++strain_column) { Structure plus = structure; Structure minus = structure; - for (std::size_t atom = 0; atom < structure.positions.size(); ++atom) - { + for (std::size_t atom = 0; atom < structure.positions.size(); ++atom) { const double displacement = step * coordinate(structure.positions[atom], strain_column); 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); } -TEST(D3Evaluator, PeriodicSmoothCutoffVirial) -{ +TEST(D3Evaluator, PeriodicSmoothCutoffVirial) { Structure structure; structure.atomic_numbers = {6, 8, 7}; structure.positions = {{1.1, 1.4, 1.7}, {5.0, 1.8, 2.1}, {2.6, 5.4, 4.9}}; @@ -567,56 +478,45 @@ TEST(D3Evaluator, PeriodicSmoothCutoffVirial) const double step = 1.0e-5; double maximum_gradient_error = 0.0; - for (std::size_t atom = 0; atom < structure.positions.size(); ++atom) - { - for (int component = 0; component < 3; ++component) - { + for (std::size_t atom = 0; atom < structure.positions.size(); ++atom) { + for (int component = 0; component < 3; ++component) { Structure plus = structure; 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); double maximum_virial_error = 0.0; - for (int strain_row = 0; strain_row < 3; ++strain_row) - { - for (int strain_column = 0; strain_column < 3; ++strain_column) - { + for (int strain_row = 0; strain_row < 3; ++strain_row) { + for (int strain_column = 0; strain_column < 3; ++strain_column) { Structure plus = structure; Structure minus = structure; - for (std::size_t atom = 0; atom < structure.positions.size(); ++atom) - { + for (std::size_t atom = 0; atom < structure.positions.size(); ++atom) { const double displacement = step * coordinate(structure.positions[atom], strain_column); coordinate(plus.positions[atom], strain_row) += displacement; coordinate(minus.positions[atom], strain_row) -= displacement; } - for (int vector = 0; vector < 3; ++vector) - { + for (int vector = 0; vector < 3; ++vector) { const double displacement = step * coordinate(structure.lattice[vector], strain_column); 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); } -TEST(D3Evaluator, SharpCutoffCompatibilityAndWidthValidation) -{ +TEST(D3Evaluator, SharpCutoffCompatibilityAndWidthValidation) { Structure structure; structure.atomic_numbers = {6, 8}; structure.positions = {{0.0, 0.0, 0.0}, {5.0, 0.0, 0.0}}; @@ -636,8 +536,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..f94d09b8b50 100644 --- a/source/source_hamilt/module_vdw/vdw.cpp +++ b/source/source_hamilt/module_vdw/vdw.cpp @@ -1,92 +1,75 @@ -#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) -{ - if (xc_input != "default") - { +std::string parse_xcname(const std::string& xc_input, const std::vector& xc_psp) { + if (xc_input != "default") { return xc_input; } - if (xc_psp.size() <= 0) - { - ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::parse_xcname", - "XC name automatic inference failed: no pseudopotential files are found"); + if (xc_psp.size() <= 0) { + 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"); + + 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"); } 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 -{ +namespace vdw { -std::unique_ptr make_vdw(const UnitCell &ucell, - const Input_para &input, - std::ofstream* plog) -{ - // NOTE: the following lines are incorrect! +std::unique_ptr make_vdw(const UnitCell& ucell, const Input_para& input, std::ofstream* plog) { + // 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") - { + if (input.vdw_method == "d2") { std::unique_ptr vdw_ptr = make_unique_compat(ucell); vdw_ptr->parameter().initial_parameters(input, plog); vdw_ptr->parameter().initset(ucell); return vdw_ptr; - } - else if (input.vdw_method == "d3_0" || input.vdw_method == "d3_bj") - { + } else if (input.vdw_method == "d3_0" || input.vdw_method == "d3_bj") { std::vector xc_psp(ucell.ntype); - for (int it = 0; it < ucell.ntype; it++) - { + for (int it = 0; it < ucell.ntype; it++) { xc_psp[it] = ucell.atoms[it].ncpp.xc_func; } - return make_unique_compat(ucell, - parse_xcname(input.dft_functional, xc_psp), - input, - plog); - } - else if (input.vdw_method == "d4") - { + return make_unique_compat(ucell, parse_xcname(input.dft_functional, xc_psp), input, plog); + } else if (input.vdw_method == "d4") { #ifdef __DFTD4 std::vector xc_psp(ucell.ntype); - for (int it = 0; it < ucell.ntype; it++) - { + for (int it = 0; it < ucell.ntype; it++) { xc_psp[it] = ucell.atoms[it].ncpp.xc_func; } std::string xc_name = input.vdw_d4_xc; - if (xc_name == "default") - { + if (xc_name == "default") { xc_name = parse_xcname(input.dft_functional, xc_psp); } @@ -96,11 +79,9 @@ std::unique_ptr make_vdw(const UnitCell &ucell, "DFT-D4 support was not enabled at build time. " "Rebuild ABACUS with -DENABLE_DFTD4=ON."); #endif - } - else if (input.vdw_method != "none") - { - ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::make_vdw", - "Unrecognized Van der Waals correction method: " + input.vdw_method); + } else if (input.vdw_method != "none") { + 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..2cad1142253 100644 --- a/source/source_hamilt/module_vdw/vdw.h +++ b/source/source_hamilt/module_vdw/vdw.h @@ -1,37 +1,30 @@ #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 -{ +namespace vdw { template -std::unique_ptr make_unique_compat(Args&&... args) -{ +std::unique_ptr make_unique_compat(Args&&... args) { return std::unique_ptr(new T(std::forward(args)...)); } -struct VdwRequest -{ +struct VdwRequest { VdwRequest(const bool force_in, const bool stress_in) : force(force_in), stress(stress_in) {} bool force; bool stress; }; -struct VdwResult -{ - VdwResult() : energy(0.0), has_force(false), has_stress(false) - { - stress.Zero(); - } +struct VdwResult { + VdwResult() : energy(0.0), has_force(false), has_stress(false) { stress.Zero(); } double energy; std::vector> force; @@ -40,15 +33,13 @@ struct VdwResult bool has_stress; }; -class Vdw -{ +class Vdw { public: Vdw(const UnitCell& unit_in) : ucell_(unit_in) {} virtual ~Vdw() = default; - VdwResult evaluate(const VdwRequest& request) - { + VdwResult evaluate(const VdwRequest& request) { VdwResult result; evaluate_impl(request, result); return result; @@ -68,9 +59,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..d6efefa2525 100644 --- a/source/source_hamilt/module_vdw/vdw_parameters.h +++ b/source/source_hamilt/module_vdw/vdw_parameters.h @@ -5,17 +5,15 @@ #include -namespace vdw -{ +namespace vdw { -class VdwParameters -{ +class VdwParameters { public: 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..1c378c3016c 100644 --- a/source/source_hamilt/module_vdw/vdw_xcname.cpp +++ b/source/source_hamilt/module_vdw/vdw_xcname.cpp @@ -4,36 +4,30 @@ #include #include -namespace vdw -{ -namespace -{ +namespace vdw { +namespace { -std::string trim_ascii(const std::string& value) -{ +std::string trim_ascii(const std::string& value) { const auto first = std::find_if_not(value.begin(), value.end(), [](unsigned char character) { return std::isspace(character) != 0; }); - if (first == value.end()) - { + if (first == value.end()) { return std::string(); } 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); } -std::string normalize_component(std::string value) -{ +std::string normalize_component(std::string value) { value = trim_ascii(value); std::transform(value.begin(), value.end(), value.begin(), [](unsigned char character) { return static_cast(std::tolower(character)); }); - if (value.compare(0, 3, "xc_") == 0) - { + if (value.compare(0, 3, "xc_") == 0) { value.erase(0, 3); } return value; @@ -41,19 +35,16 @@ std::string normalize_component(std::string value) } // namespace -std::string normalize_xc_name(const std::string& input) -{ +std::string normalize_xc_name(const std::string& input) { const std::size_t plus = input.find('+'); const std::size_t colon = input.find(':'); const std::size_t separator = plus != std::string::npos ? plus : colon; - if (separator == std::string::npos) - { + if (separator == std::string::npos) { 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/vdw_xcname.h b/source/source_hamilt/module_vdw/vdw_xcname.h index 865b77a5224..f2d7ea01923 100644 --- a/source/source_hamilt/module_vdw/vdw_xcname.h +++ b/source/source_hamilt/module_vdw/vdw_xcname.h @@ -3,8 +3,7 @@ #include -namespace vdw -{ +namespace vdw { /** * @brief Normalize an XC-functional spelling at the ABACUS/libXC boundary. diff --git a/source/source_hamilt/module_vdw/vdwd2.cpp b/source/source_hamilt/module_vdw/vdwd2.cpp index 9489ea62a75..d169c40312b 100644 --- a/source/source_hamilt/module_vdw/vdwd2.cpp +++ b/source/source_hamilt/module_vdw/vdwd2.cpp @@ -5,27 +5,24 @@ //========================================================== #include "vdwd2.h" + #include "source_base/timer.h" #include #include -namespace vdw -{ +namespace vdw { -void Vdwd2::evaluate_impl(const VdwRequest& request, VdwResult& result) -{ +void Vdwd2::evaluate_impl(const VdwRequest& request, VdwResult& result) { ModuleBase::TITLE("Vdwd2", "evaluate"); ModuleBase::timer::start("Vdwd2", "evaluate"); para_.initset(ucell_); - if (request.force) - { + if (request.force) { result.force.resize(ucell_.nat); } - if (request.stress) - { + if (request.stress) { result.stress.Zero(); } @@ -43,22 +40,18 @@ void Vdwd2::evaluate_impl(const VdwRequest& request, VdwResult& result) result.energy -= C6_product / pow(r_sqr, 3) / tmp_damp_recip / 2.0; - if (!need_derivatives) - { + if (!need_derivatives) { return; } 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) - { + if (request.force) { result.force[ucell_.itia2iat(it1, ia1)] += tmp_factor * (tau1 - tau2); } - if (request.stress) - { + if (request.stress) { const ModuleBase::Vector3 dr = tau2 - tau1; result.stress += tmp_factor / 2.0 * ModuleBase::Matrix3(dr.x * dr.x, @@ -77,16 +70,14 @@ void Vdwd2::evaluate_impl(const VdwRequest& request, VdwResult& result) result.energy *= para_.scaling(); - if (request.force) - { + if (request.force) { std::for_each(result.force.begin(), result.force.end(), [&](ModuleBase::Vector3& force) { force *= para_.scaling() / ucell_.lat0; }); result.has_force = true; } - if (request.stress) - { + if (request.stress) { result.stress *= para_.scaling() / ucell_.omega; result.has_stress = true; } diff --git a/source/source_hamilt/module_vdw/vdwd2.h b/source/source_hamilt/module_vdw/vdwd2.h index 43d93615db5..6432d56c4b7 100644 --- a/source/source_hamilt/module_vdw/vdwd2.h +++ b/source/source_hamilt/module_vdw/vdwd2.h @@ -10,11 +10,9 @@ #include "vdw.h" #include "vdwd2_parameters.h" -namespace vdw -{ +namespace vdw { -class Vdwd2 : public Vdw -{ +class Vdwd2 : public Vdw { public: Vdwd2(const UnitCell& unit_in) : Vdw(unit_in) {} @@ -30,22 +28,17 @@ class Vdwd2 : public Vdw void evaluate_impl(const VdwRequest& request, VdwResult& result) override; template - void index_loops(F&& f) - { + void index_loops(F&& f) { int xidx = para_.period().x / 2; int yidx = para_.period().y / 2; int zidx = para_.period().z / 2; - for (int it1 = 0; it1 < ucell_.ntype; ++it1) - { - if (ucell_.atoms[it1].flag_empty_element) - { + for (int it1 = 0; it1 < ucell_.ntype; ++it1) { + if (ucell_.atoms[it1].flag_empty_element) { continue; } - for (int it2 = 0; it2 < ucell_.ntype; ++it2) - { - if (ucell_.atoms[it2].flag_empty_element) - { + for (int it2 = 0; it2 < ucell_.ntype; ++it2) { + if (ucell_.atoms[it2].flag_empty_element) { continue; } const double C6_product @@ -54,25 +47,18 @@ class Vdwd2 : public Vdw const double R0_sum = (para_.R0().at(ucell_.atoms[it1].ncpp.psd) + para_.R0().at(ucell_.atoms[it2].ncpp.psd)) / ucell_.lat0; - if (!R0_sum) - { + if (!R0_sum) { ModuleBase::WARNING_QUIT("Input", "R0_sum can not be 0"); } - for (int ia1 = 0; ia1 != ucell_.atoms[it1].na; ++ia1) - { + for (int ia1 = 0; ia1 != ucell_.atoms[it1].na; ++ia1) { const ModuleBase::Vector3 tau1 = ucell_.atoms[it1].tau[ia1]; - for (int ia2 = 0; ia2 != ucell_.atoms[it2].na; ++ia2) - { + for (int ia2 = 0; ia2 != ucell_.atoms[it2].na; ++ia2) { ModuleBase::Vector3 ilat_loop; - for (ilat_loop.x = -xidx; ilat_loop.x <= xidx; ++ilat_loop.x) - { - for (ilat_loop.y = -yidx; ilat_loop.y <= yidx; ++ilat_loop.y) - { - for (ilat_loop.z = -zidx; ilat_loop.z <= zidx; ++ilat_loop.z) - { + for (ilat_loop.x = -xidx; ilat_loop.x <= xidx; ++ilat_loop.x) { + for (ilat_loop.y = -yidx; ilat_loop.y <= yidx; ++ilat_loop.y) { + 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)) - { + && (ia1 == ia2)) { continue; } const ModuleBase::Vector3 tau2 diff --git a/source/source_hamilt/module_vdw/vdwd2_parameters.cpp b/source/source_hamilt/module_vdw/vdwd2_parameters.cpp index 3ac48e18199..c90402c2f0b 100644 --- a/source/source_hamilt/module_vdw/vdwd2_parameters.cpp +++ b/source/source_hamilt/module_vdw/vdwd2_parameters.cpp @@ -6,45 +6,35 @@ #include "vdwd2_parameters.h" -namespace vdw -{ +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; C6_input(input.vdw_C6_file, input.vdw_C6_unit); R0_input(input.vdw_R0_file, input.vdw_R0_unit); model_ = input.vdw_cutoff_type; - if (input.vdw_cutoff_type == "radius") - { + if (input.vdw_cutoff_type == "radius") { radius_ = std::stod(input.vdw_cutoff_radius); - if (!(input.vdw_radius_unit == "Bohr")) - { + if (!(input.vdw_radius_unit == "Bohr")) { radius_ /= ModuleBase::BOHR_TO_A; } - } - else if (input.vdw_cutoff_type == "period") - { + } else if (input.vdw_cutoff_type == "period") { period_ = input.vdw_cutoff_period; } } -void Vdwd2Parameters::initset(const UnitCell &ucell) -{ - if (model_ == "radius") - { +void Vdwd2Parameters::initset(const UnitCell& ucell) { + if (model_ == "radius") { period_.x = 2 * ceil(radius_ / ucell.lat0 / sqrt(ucell.a1.norm2())) + 1; period_.y = 2 * ceil(radius_ / ucell.lat0 / sqrt(ucell.a2.norm2())) + 1; period_.z = 2 * ceil(radius_ / ucell.lat0 / sqrt(ucell.a3.norm2())) + 1; } } -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") - { + if (file != "default") { std::ifstream ifs(file); if (!ifs) ModuleBase::WARNING_QUIT("Vdwd2::C6_input", @@ -55,8 +45,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) / ModuleBase::Ry_to_eV; @@ -67,11 +56,9 @@ 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") - { + if (file != "default") { std::ifstream ifs(file.c_str()); if (!ifs) ModuleBase::WARNING_QUIT("Vdwd2::R0_input", @@ -82,8 +69,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; else if (unit == "Bohr") diff --git a/source/source_hamilt/module_vdw/vdwd2_parameters.h b/source/source_hamilt/module_vdw/vdwd2_parameters.h index a22fe350320..ca231c188ed 100644 --- a/source/source_hamilt/module_vdw/vdwd2_parameters.h +++ b/source/source_hamilt/module_vdw/vdwd2_parameters.h @@ -7,20 +7,18 @@ #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" -namespace vdw -{ +#include + +namespace vdw { -class Vdwd2Parameters : public VdwParameters -{ +class Vdwd2Parameters : public VdwParameters { public: - Vdwd2Parameters() : VdwParameters() - { + Vdwd2Parameters() : VdwParameters() { C6_ = C6_default_; R0_ = R0_default_; damping_ = 0.0; @@ -30,19 +28,20 @@ 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_; } diff --git a/source/source_hamilt/module_vdw/vdwd3.cpp b/source/source_hamilt/module_vdw/vdwd3.cpp index 923ea619300..9f00da9410c 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 @@ -13,125 +13,82 @@ #include #include -namespace vdw -{ -namespace -{ +namespace vdw { +namespace { -int atomic_number_from_symbol(const std::string& symbol) -{ - for (std::size_t index = 0; index < ModuleBase::element_name.size(); ++index) - { - if (symbol == ModuleBase::element_name[index]) - { +int atomic_number_from_symbol(const std::string& symbol) { + for (std::size_t index = 0; index < ModuleBase::element_name.size(); ++index) { + if (symbol == ModuleBase::element_name[index]) { return static_cast(index) + 1; } } ModuleBase::WARNING_QUIT("Vdwd3::atomic_number_from_symbol", "Unknown element symbol: " + symbol); } -double length_to_bohr(double value, const std::string& unit) -{ - if (unit == "Bohr") - { +double length_to_bohr(double value, const std::string& unit) { + if (unit == "Bohr") { return value; } - if (unit == "A") - { + if (unit == "A") { return value / ModuleBase::BOHR_TO_A; } ModuleBase::WARNING_QUIT("Vdwd3::length_to_bohr", "Unsupported length unit: " + unit); } -double cutoff_to_bohr(const std::string& value, const std::string& unit) -{ +double cutoff_to_bohr(const std::string& value, const std::string& unit) { return length_to_bohr(std::stod(value), unit); } -d3::Vec3 to_d3_vector(const ModuleBase::Vector3& value) -{ - return d3::Vec3(value.x, value.y, value.z); -} +d3::Vec3 to_d3_vector(const ModuleBase::Vector3& value) { return d3::Vec3(value.x, value.y, value.z); } -const char* damping_name(d3::Damping damping) -{ - return damping == d3::Damping::Rational ? "rational (BJ)" : "zero"; -} +const char* damping_name(d3::Damping damping) { return damping == d3::Damping::Rational ? "rational (BJ)" : "zero"; } } // namespace -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") - { +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") { ModuleBase::WARNING_QUIT("Vdwd3::Vdwd3", "DFT-D3 requires vdw_cutoff_type=radius; " "the legacy period-count cutoff is no longer supported."); } - if (input.vdw_method == "d3_bj") - { + if (input.vdw_method == "d3_bj") { parameters_.damping = d3::Damping::Rational; - } - else if (input.vdw_method == "d3_0") - { + } else if (input.vdw_method == "d3_0") { parameters_.damping = d3::Damping::Zero; - } - else - { + } else { 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) - { + 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") - { + if (input.vdw_s6 != "default") { parameters_.s6 = std::stod(input.vdw_s6); } - if (input.vdw_s8 != "default") - { + if (input.vdw_s8 != "default") { parameters_.s8 = std::stod(input.vdw_s8); } - if (input.vdw_a1 != "default") - { - if (parameters_.damping == d3::Damping::Rational) - { + if (input.vdw_a1 != "default") { + if (parameters_.damping == d3::Damping::Rational) { parameters_.a1 = std::stod(input.vdw_a1); - } - else - { + } else { parameters_.rs6 = std::stod(input.vdw_a1); } } - if (input.vdw_a2 != "default") - { - if (parameters_.damping == d3::Damping::Rational) - { + if (input.vdw_a2 != "default") { + if (parameters_.damping == d3::Damping::Rational) { parameters_.a2 = std::stod(input.vdw_a2); - } - else - { + } else { parameters_.rs8 = std::stod(input.vdw_a2); } } @@ -143,14 +100,12 @@ Vdwd3::Vdwd3(const UnitCell& unit_in, cutoffs_.width2 = input.vdw_cutoff_width2; cutoffs_.width3 = input.vdw_cutoff_width3; - if (cutoffs_.width2 < 0.0 || cutoffs_.width2 > cutoffs_.disp2) - { + if (cutoffs_.width2 < 0.0 || cutoffs_.width2 > cutoffs_.disp2) { ModuleBase::WARNING_QUIT("Vdwd3::Vdwd3", "vdw_cutoff_width2 must satisfy " "0 <= width <= two-body cutoff"); } - if (cutoffs_.width3 < 0.0 || cutoffs_.width3 > cutoffs_.disp3) - { + if (cutoffs_.width3 < 0.0 || cutoffs_.width3 > cutoffs_.disp3) { ModuleBase::WARNING_QUIT("Vdwd3::Vdwd3", "vdw_cutoff_width3 must satisfy " "0 <= width <= three-body cutoff"); @@ -159,35 +114,25 @@ Vdwd3::Vdwd3(const UnitCell& unit_in, write_parameters(plog); } -void Vdwd3::write_parameters(std::ofstream* plog) const -{ - if (plog == nullptr) - { +void Vdwd3::write_parameters(std::ofstream* plog) const { + if (plog == nullptr) { return; } *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; - if (parameters_.damping == d3::Damping::Rational) - { + << 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; - } - else - { + } else { *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 -{ +d3::Structure Vdwd3::build_structure(std::vector& atom_indices) const { d3::Structure structure; atom_indices.clear(); structure.atomic_numbers.reserve(ucell_.nat); @@ -195,17 +140,14 @@ d3::Structure Vdwd3::build_structure(std::vector& atom_indices) const atom_indices.reserve(ucell_.nat); int iat = 0; - for (int it = 0; it < ucell_.ntype; ++it) - { - if (ucell_.atoms[it].flag_empty_element) - { + for (int it = 0; it < ucell_.ntype; ++it) { + if (ucell_.atoms[it].flag_empty_element) { iat += ucell_.atoms[it].na; continue; } const int atomic_number = atomic_number_from_symbol(ucell_.atoms[it].ncpp.psd); - for (int ia = 0; ia < ucell_.atoms[it].na; ++ia) - { + for (int ia = 0; ia < ucell_.atoms[it].na; ++ia) { structure.atomic_numbers.push_back(atomic_number); structure.positions.push_back(to_d3_vector(ucell_.atoms[it].tau[ia] * ucell_.lat0)); atom_indices.push_back(iat); @@ -220,8 +162,7 @@ d3::Structure Vdwd3::build_structure(std::vector& atom_indices) const return structure; } -void Vdwd3::evaluate_impl(const VdwRequest& request, VdwResult& result) -{ +void Vdwd3::evaluate_impl(const VdwRequest& request, VdwResult& result) { ModuleBase::TITLE("Vdwd3", "evaluate"); ModuleBase::timer::start("Vdwd3", "evaluate"); @@ -230,19 +171,16 @@ void Vdwd3::evaluate_impl(const VdwRequest& request, VdwResult& result) const bool derivatives = request.force || request.stress; std::vector atom_indices; const d3::Structure structure = build_structure(atom_indices); - if (!d3::evaluate(structure, parameters_, cutoffs_, derivatives, d3_result, error)) - { + if (!d3::evaluate(structure, parameters_, cutoffs_, derivatives, d3_result, error)) { ModuleBase::WARNING_QUIT("Vdwd3::evaluate", error); } // The native core follows s-dftd3 and returns Hartree-based quantities. result.energy = 2.0 * d3_result.energy; - if (request.force) - { + if (request.force) { result.force.assign(ucell_.nat, ModuleBase::Vector3(0.0, 0.0, 0.0)); - for (std::size_t index = 0; index < atom_indices.size(); ++index) - { + for (std::size_t index = 0; index < atom_indices.size(); ++index) { const int iat = atom_indices[index]; result.force[iat] = ModuleBase::Vector3(-2.0 * d3_result.gradient[index].x, -2.0 * d3_result.gradient[index].y, @@ -251,8 +189,7 @@ void Vdwd3::evaluate_impl(const VdwRequest& request, VdwResult& result) result.has_force = true; } - if (request.stress) - { + if (request.stress) { const d3::Matrix3& sigma = d3_result.virial; result.stress = ModuleBase::Matrix3(2.0 * sigma.value[0][0], 2.0 * sigma.value[0][1], diff --git a/source/source_hamilt/module_vdw/vdwd3.h b/source/source_hamilt/module_vdw/vdwd3.h index e8a1e50a9e7..b4363b65383 100644 --- a/source/source_hamilt/module_vdw/vdwd3.h +++ b/source/source_hamilt/module_vdw/vdwd3.h @@ -1,22 +1,17 @@ #ifndef VDWD3_H #define VDWD3_H -#include "vdwd3_types.h" #include "vdw.h" +#include "vdwd3_types.h" #include #include -namespace vdw -{ +namespace vdw { -class Vdwd3 : public 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..5415eafbd23 100644 --- a/source/source_hamilt/module_vdw/vdwd3_data.cpp +++ b/source/source_hamilt/module_vdw/vdwd3_data.cpp @@ -3,44 +3,31 @@ #include #include -namespace vdw -{ -namespace d3 -{ -namespace data -{ -namespace -{ +namespace vdw { +namespace d3 { +namespace data { +namespace { #include "data/d3_reference.inc" -inline std::size_t pair_index(int high, int low) -{ - return static_cast(high * (high - 1) / 2 + low - 1); -} +inline std::size_t pair_index(int high, int low) { return static_cast(high * (high - 1) / 2 + low - 1); } -inline void check_atomic_number(int atomic_number) -{ - assert(atomic_number > 0 && atomic_number <= max_element); -} +inline void check_atomic_number(int atomic_number) { assert(atomic_number > 0 && atomic_number <= max_element); } } // namespace -int reference_count(int atomic_number) -{ +int reference_count(int atomic_number) { check_atomic_number(atomic_number); return static_cast(kReferenceCount[atomic_number]); } -double reference_cn(int atomic_number, int reference) -{ +double reference_cn(int atomic_number, int reference) { check_atomic_number(atomic_number); assert(reference >= 0 && reference < reference_count(atomic_number)); return kReferenceCn[atomic_number * max_reference + reference]; } -double reference_c6(int atomic_number_i, int reference_i, int atomic_number_j, int reference_j) -{ +double reference_c6(int atomic_number_i, int reference_i, int atomic_number_j, int reference_j) { check_atomic_number(atomic_number_i); check_atomic_number(atomic_number_j); assert(reference_i >= 0 && reference_i < reference_count(atomic_number_i)); @@ -50,8 +37,7 @@ double reference_c6(int atomic_number_i, int reference_i, int atomic_number_j, i int low = atomic_number_j; int high_reference = reference_i; int low_reference = reference_j; - if (high < low) - { + if (high < low) { const int atomic_number = high; high = low; low = atomic_number; @@ -61,25 +47,21 @@ 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]; } -double covalent_radius(int atomic_number) -{ +double covalent_radius(int atomic_number) { check_atomic_number(atomic_number); return kCovalentRadius[atomic_number]; } -double r4r2(int atomic_number) -{ +double r4r2(int atomic_number) { check_atomic_number(atomic_number); return kR4R2[atomic_number]; } -double vdw_radius(int atomic_number_i, int atomic_number_j) -{ +double vdw_radius(int atomic_number_i, int atomic_number_j) { check_atomic_number(atomic_number_i); check_atomic_number(atomic_number_j); const int high = atomic_number_i > atomic_number_j ? atomic_number_i : atomic_number_j; diff --git a/source/source_hamilt/module_vdw/vdwd3_data.h b/source/source_hamilt/module_vdw/vdwd3_data.h index ba6de0a8d7c..73a8212b312 100644 --- a/source/source_hamilt/module_vdw/vdwd3_data.h +++ b/source/source_hamilt/module_vdw/vdwd3_data.h @@ -1,12 +1,9 @@ #ifndef ABACUS_D3_DATA_H #define ABACUS_D3_DATA_H -namespace vdw -{ -namespace d3 -{ -namespace data -{ +namespace vdw { +namespace d3 { +namespace data { constexpr int max_element = 103; constexpr int max_reference = 7; diff --git a/source/source_hamilt/module_vdw/vdwd3_evaluator.cpp b/source/source_hamilt/module_vdw/vdwd3_evaluator.cpp index f130971b37b..145097b6193 100644 --- a/source/source_hamilt/module_vdw/vdwd3_evaluator.cpp +++ b/source/source_hamilt/module_vdw/vdwd3_evaluator.cpp @@ -11,88 +11,61 @@ #include #include -namespace vdw -{ -namespace d3 -{ -namespace -{ +namespace vdw { +namespace d3 { +namespace { constexpr double kCoordinationSteepness = 16.0; constexpr double kReferenceWeight = 4.0; constexpr double kAtmRadiusScale = 4.0 / 3.0; -double dot(const Vec3& lhs, const Vec3& rhs) -{ - return lhs.x * rhs.x + lhs.y * rhs.y + lhs.z * rhs.z; -} +double dot(const Vec3& lhs, const Vec3& rhs) { return lhs.x * rhs.x + lhs.y * rhs.y + lhs.z * rhs.z; } -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); +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); } -double norm(const Vec3& value) -{ - return std::sqrt(dot(value, value)); -} +double norm(const Vec3& value) { return std::sqrt(dot(value, value)); } -void add_outer(Matrix3& matrix, const Vec3& lhs, const Vec3& rhs, double scale = 1.0) -{ +void add_outer(Matrix3& matrix, const Vec3& lhs, const Vec3& rhs, double scale = 1.0) { const double left[3] = {lhs.x, lhs.y, lhs.z}; const double right[3] = {rhs.x, rhs.y, rhs.z}; - for (int row = 0; row < 3; ++row) - { - for (int column = 0; column < 3; ++column) - { + for (int row = 0; row < 3; ++row) { + for (int column = 0; column < 3; ++column) { matrix.value[row][column] += scale * left[row] * right[column]; } } } -bool validate(const Structure& structure, const Cutoffs& cutoffs, std::string& error) -{ - if (structure.atomic_numbers.size() != structure.positions.size()) - { +bool validate(const Structure& structure, const Cutoffs& cutoffs, std::string& error) { + if (structure.atomic_numbers.size() != structure.positions.size()) { error = "atomic-number and position arrays have different sizes"; return false; } - for (int atomic_number : structure.atomic_numbers) - { - if (atomic_number <= 0 || atomic_number > data::max_element) - { + for (int atomic_number: structure.atomic_numbers) { + if (atomic_number <= 0 || atomic_number > data::max_element) { error = "DFT-D3 supports atomic numbers 1 through 103"; 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; } return true; } -std::vector lattice_points(const Structure& structure, double cutoff) -{ - if (!structure.periodic[0] && !structure.periodic[1] && !structure.periodic[2]) - { +std::vector lattice_points(const Structure& structure, double cutoff) { + if (!structure.periodic[0] && !structure.periodic[1] && !structure.periodic[2]) { return std::vector(1, Vec3()); } @@ -100,25 +73,18 @@ std::vector lattice_points(const Structure& structure, double cutoff) cross(structure.lattice[2], structure.lattice[0]), cross(structure.lattice[0], structure.lattice[1])}}; std::array repeat; - for (int direction = 0; direction < 3; ++direction) - { + 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))); - 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.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 + structure.lattice[2] * iz); } } @@ -126,76 +92,57 @@ std::vector lattice_points(const Structure& structure, double cutoff) return translations; } -void smooth_cutoff(double distance, double cutoff, double width, double& value, double& derivative) -{ +void smooth_cutoff(double distance, double cutoff, double width, double& value, double& derivative) { // Match s-dftd3 exactly: width == cutoff selects the legacy sharp path. - if (width <= 0.0 || width >= cutoff) - { + if (width <= 0.0 || width >= cutoff) { value = 1.0; derivative = 0.0; return; } const double inner = cutoff - width; - if (distance <= inner) - { + if (distance <= inner) { value = 1.0; derivative = 0.0; - } - else if (distance >= cutoff) - { + } else if (distance >= cutoff) { value = 0.0; derivative = 0.0; - } - else - { + } else { const double x = (cutoff - distance) / width; value = x * x * x * (10.0 + x * (-15.0 + 6.0 * x)); derivative = -30.0 * x * x * (1.0 - x) * (1.0 - x) / width; } } -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); +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); 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); +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 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, - const std::vector& translations, - double cutoff) -{ +std::vector + coordination_numbers(const Structure& structure, const std::vector& translations, double cutoff) { const std::size_t atoms = structure.positions.size(); const double cutoff2 = cutoff * cutoff; std::vector coordination(atoms, 0.0); - for (std::size_t i = 0; i < atoms; ++i) - { - for (std::size_t j = 0; j <= i; ++j) - { - for (const Vec3& translation : translations) - { + for (std::size_t i = 0; i < atoms; ++i) { + for (std::size_t j = 0; j <= i; ++j) { + for (const Vec3& translation: translations) { const Vec3 vector = structure.positions[i] - (structure.positions[j] + translation); const double distance2 = dot(vector, vector); - if (distance2 > cutoff2 || distance2 < 1.0e-12) - { + if (distance2 > cutoff2 || distance2 < 1.0e-12) { continue; } const double count = coordination_count(structure.atomic_numbers[i], structure.atomic_numbers[j], std::sqrt(distance2)); coordination[i] += count; - if (i != j) - { + if (i != j) { coordination[j] += count; } } @@ -208,55 +155,48 @@ void reference_weights(const Structure& structure, const std::vector& coordination, bool derivatives, std::vector& weights, - std::vector& weight_derivatives) -{ + std::vector& weight_derivatives) { const std::size_t atoms = structure.positions.size(); weights.assign(atoms * data::max_reference, 0.0); - if (derivatives) - { + if (derivatives) { weight_derivatives.assign(atoms * data::max_reference, 0.0); } - for (std::size_t atom = 0; atom < atoms; ++atom) - { + for (std::size_t atom = 0; atom < atoms; ++atom) { const int atomic_number = structure.atomic_numbers[atom]; const int references = data::reference_count(atomic_number); double raw[data::max_reference]; double max_exponent = -std::numeric_limits::infinity(); - for (int reference = 0; reference < references; ++reference) - { + for (int reference = 0; reference < references; ++reference) { const double delta = data::reference_cn(atomic_number, reference) - coordination[atom]; raw[reference] = -kReferenceWeight * delta * delta; 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. + for (int reference = 0; reference < references; ++reference) { + // Shifting by the maximum exponent prevents all weights from + // underflowing. raw[reference] = std::exp(raw[reference] - max_exponent); weight_sum += raw[reference]; - if (derivatives) - { + if (derivatives) { weighted_reference_cn += raw[reference] * data::reference_cn(atomic_number, reference); } } const double inverse_weight_sum = 1.0 / weight_sum; weighted_reference_cn *= inverse_weight_sum; - for (int reference = 0; reference < references; ++reference) - { + for (int reference = 0; reference < references; ++reference) { const std::size_t index = atom * data::max_reference + reference; const double weight = raw[reference] * inverse_weight_sum; weights[index] = weight; - if (derivatives) - { + 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); } } } @@ -267,49 +207,40 @@ void atomic_c6(const Structure& structure, const std::vector& weight_derivatives, bool derivatives, std::vector& c6, - std::vector& dc6dcn) -{ + std::vector& dc6dcn) { const std::size_t atoms = structure.positions.size(); c6.assign(atoms * atoms, 0.0); - if (derivatives) - { + if (derivatives) { dc6dcn.assign(atoms * atoms, 0.0); } - for (std::size_t i = 0; i < atoms; ++i) - { + for (std::size_t i = 0; i < atoms; ++i) { const int references_i = data::reference_count(structure.atomic_numbers[i]); - for (std::size_t j = 0; j <= i; ++j) - { + for (std::size_t j = 0; j <= i; ++j) { const int references_j = data::reference_count(structure.atomic_numbers[j]); double coefficient = 0.0; double derivative_i = 0.0; double derivative_j = 0.0; - for (int reference_i = 0; reference_i < references_i; ++reference_i) - { - for (int reference_j = 0; reference_j < references_j; ++reference_j) - { + for (int reference_i = 0; reference_i < references_i; ++reference_i) { + 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; + 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; } } } c6[i * atoms + j] = coefficient; c6[j * atoms + i] = coefficient; - if (derivatives) - { + if (derivatives) { dc6dcn[i * atoms + j] = derivative_i; dc6dcn[j * atoms + i] = derivative_j; } @@ -326,17 +257,14 @@ void pairwise_dispersion(const Structure& structure, bool derivatives, std::vector& atomic_energy, std::vector& dEdcn, - Result& result) -{ + Result& result) { const std::size_t atoms = structure.positions.size(); const double cutoff2 = cutoffs.disp2 * cutoffs.disp2; const double epsilon = std::numeric_limits::epsilon(); - for (std::size_t i = 0; i < atoms; ++i) - { + for (std::size_t i = 0; i < atoms; ++i) { const int atomic_number_i = structure.atomic_numbers[i]; - for (std::size_t j = 0; j <= i; ++j) - { + for (std::size_t j = 0; j <= i; ++j) { const int atomic_number_j = structure.atomic_numbers[j]; const double rrij = 3.0 * data::r4r2(atomic_number_i) * data::r4r2(atomic_number_j); const double coefficient = c6[j * atoms + i]; @@ -344,52 +272,37 @@ void pairwise_dispersion(const Structure& structure, double r0 = 0.0; double r0_sixth = 0.0; double r0_eighth = 0.0; - if (parameters.damping == Damping::Rational) - { + if (parameters.damping == Damping::Rational) { r0 = parameters.a1 * std::sqrt(rrij) + parameters.a2; const double r0_squared = r0 * r0; r0_sixth = r0_squared * r0_squared * r0_squared; r0_eighth = r0_sixth * r0_squared; - } - else - { + } else { 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); - if (distance2 > cutoff2 || distance2 < epsilon) - { + if (distance2 > cutoff2 || distance2 < epsilon) { continue; } 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; - if (parameters.damping == Damping::Rational) - { + if (parameters.damping == Damping::Rational) { const double fourth = distance2 * distance2; 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; - } - else - { + 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 { const double sixth = distance2 * distance2 * distance2; const double eighth = sixth * distance2; const double exponent6 = parameters.alp; @@ -398,33 +311,25 @@ 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; const double energy = -coefficient * dispersion * 0.5; atomic_energy[i] += energy; - if (i != j) - { + if (i != j) { atomic_energy[j] += energy; } - if (derivatives) - { - const double radial = switching * bare_radial - + switching_derivative * bare_energy / distance; + if (derivatives) { + 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); - if (i != j) - { + if (i != j) { dEdcn[j] -= dc6dcn[j * atoms + i] * dispersion; result.gradient[i] += gradient; result.gradient[j] -= gradient; @@ -436,10 +341,8 @@ void pairwise_dispersion(const Structure& structure, } } -double triple_scale(std::size_t i, std::size_t j, std::size_t k) -{ - if (i == j) - { +double triple_scale(std::size_t i, std::size_t j, std::size_t k) { + if (i == j) { return i == k ? 1.0 / 6.0 : 0.5; } return i != k && j != k ? 1.0 : 0.5; @@ -454,10 +357,8 @@ void atm_dispersion(const Structure& structure, bool derivatives, std::vector& atomic_energy, std::vector& dEdcn, - Result& result) -{ - if (std::abs(parameters.s9) < std::numeric_limits::epsilon()) - { + Result& result) { + if (std::abs(parameters.s9) < std::numeric_limits::epsilon()) { return; } @@ -467,20 +368,16 @@ void atm_dispersion(const Structure& structure, const double alpha = parameters.alp + 2.0; const double alpha_third = alpha / 3.0; - for (std::size_t i = 0; i < atoms; ++i) - { + for (std::size_t i = 0; i < atoms; ++i) { const int atomic_number_i = structure.atomic_numbers[i]; - for (std::size_t j = 0; j <= i; ++j) - { + for (std::size_t j = 0; j <= i; ++j) { 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); - if (r2ij > cutoff2 || r2ij < epsilon) - { + if (r2ij > cutoff2 || r2ij < epsilon) { continue; } const double rij = std::sqrt(r2ij); @@ -488,25 +385,20 @@ void atm_dispersion(const Structure& structure, double dswij = 0.0; smooth_cutoff(rij, cutoffs.disp3, cutoffs.width3, swij, dswij); - for (std::size_t k = 0; k <= j; ++k) - { + for (std::size_t k = 0; k <= j; ++k) { const int atomic_number_k = structure.atomic_numbers[k]; 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); - if (r2ik > cutoff2 || r2ik < epsilon) - { + if (r2ik > cutoff2 || r2ik < epsilon) { continue; } const double rik = std::sqrt(r2ik); @@ -516,8 +408,7 @@ void atm_dispersion(const Structure& structure, const Vec3 vjk = vik - vij; const double r2jk = dot(vjk, vjk); - if (r2jk > cutoff2 || r2jk < epsilon) - { + if (r2jk > cutoff2 || r2jk < epsilon) { continue; } const double rjk = std::sqrt(r2jk); @@ -532,10 +423,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; @@ -543,50 +433,43 @@ void atm_dispersion(const Structure& structure, atomic_energy[j] -= energy / 3.0; atomic_energy[k] -= energy / 3.0; - if (!derivatives) - { + if (!derivatives) { 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 +479,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); } } } @@ -616,28 +493,22 @@ void add_coordination_derivatives(const Structure& structure, const std::vector& translations, double cutoff, const std::vector& dEdcn, - Result& result) -{ + Result& result) { const std::size_t atoms = structure.positions.size(); const double cutoff2 = cutoff * cutoff; - for (std::size_t i = 0; i < atoms; ++i) - { - for (std::size_t j = 0; j <= i; ++j) - { - for (const Vec3& translation : translations) - { + for (std::size_t i = 0; i < atoms; ++i) { + for (std::size_t j = 0; j <= i; ++j) { + for (const Vec3& translation: translations) { const Vec3 vector = structure.positions[i] - (structure.positions[j] + translation); const double distance2 = dot(vector, vector); - if (distance2 > cutoff2 || distance2 < 1.0e-12) - { + if (distance2 > cutoff2 || distance2 < 1.0e-12) { continue; } 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; @@ -655,28 +526,23 @@ bool evaluate(const Structure& structure, const Cutoffs& cutoffs, bool derivatives, Result& result, - std::string& error) -{ + std::string& error) { error.clear(); result = Result(); - if (!validate(structure, cutoffs, error)) - { + if (!validate(structure, cutoffs, error)) { return false; } const std::size_t atoms = structure.positions.size(); - if (derivatives) - { + if (derivatives) { result.gradient.assign(atoms, Vec3()); } - if (atoms == 0) - { + if (atoms == 0) { return true; } 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; @@ -712,8 +578,7 @@ bool evaluate(const Structure& structure, dEdcn, result); - if (derivatives) - { + if (derivatives) { add_coordination_derivatives(structure, cn_translations, cutoffs.cn, dEdcn, result); } result.energy = std::accumulate(atomic_energy.begin(), atomic_energy.end(), 0.0); diff --git a/source/source_hamilt/module_vdw/vdwd3_evaluator.h b/source/source_hamilt/module_vdw/vdwd3_evaluator.h index bf91ca83d24..bbcc9adddac 100644 --- a/source/source_hamilt/module_vdw/vdwd3_evaluator.h +++ b/source/source_hamilt/module_vdw/vdwd3_evaluator.h @@ -5,10 +5,8 @@ #include -namespace vdw -{ -namespace d3 -{ +namespace vdw { +namespace d3 { bool evaluate(const Structure& structure, const Parameters& parameters, diff --git a/source/source_hamilt/module_vdw/vdwd3_parameters.cpp b/source/source_hamilt/module_vdw/vdwd3_parameters.cpp index 5f6a4bffbf1..086e29e9358 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 @@ -7,15 +8,11 @@ #include #include -namespace vdw -{ -namespace d3 -{ -namespace -{ +namespace vdw { +namespace d3 { +namespace { -struct DampingParameterRecord -{ +struct DampingParameterRecord { const char* method_id; Damping damping; double s6; @@ -28,8 +25,7 @@ struct DampingParameterRecord double alp; }; -struct MethodAlias -{ +struct MethodAlias { const char* alias; const char* method_id; }; @@ -37,55 +33,46 @@ struct MethodAlias #include "data/d3_damping_parameters.inc" #include "data/d3_method_aliases.inc" -std::string lower_ascii(std::string value) -{ +std::string lower_ascii(std::string value) { std::transform(value.begin(), value.end(), value.begin(), [](unsigned char character) { return static_cast(std::tolower(character)); }); return value; } -std::string normalized(std::string value) -{ +std::string normalized(std::string value) { value = lower_ascii(value); std::string result; result.reserve(value.size()); - for (unsigned char character : value) - { - if (std::isalnum(character) || character == '/') - { + for (unsigned char character: value) { + if (std::isalnum(character) || character == '/') { result.push_back(static_cast(character)); } } return result; } -std::string libxc_component(const std::string& component, const std::string& marker) -{ +std::string libxc_component(const std::string& component, const std::string& marker) { const std::string lower = lower_ascii(component); const std::size_t position = lower.find(marker); - if (position == std::string::npos) - { + if (position == std::string::npos) { return normalized(component); } return normalized(component.substr(position + marker.size())); } -std::string canonicalize_libxc_pair(const std::string& input) -{ +std::string canonicalize_libxc_pair(const std::string& input) { const std::size_t separator = input.find(':'); std::string first = input.substr(0, separator); std::string second = input.substr(separator + 1); const std::string first_lower = lower_ascii(first); - if (first_lower.find("_c_") != std::string::npos) - { + if (first_lower.find("_c_") != std::string::npos) { std::swap(first, second); } const std::string exchange = libxc_component(first, "_x_"); const std::string correlation = libxc_component(second, "_c_"); - if (exchange == correlation) - { + if (exchange == correlation) { return exchange; } @@ -93,11 +80,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); @@ -106,31 +98,25 @@ std::string canonicalize_libxc_pair(const std::string& input) } // namespace -std::string canonicalize_method_name(const std::string& input) -{ +std::string canonicalize_method_name(const std::string& input) { const std::string xc_name = normalize_xc_name(input); std::string canonical; - if (xc_name.find(':') != std::string::npos) - { + if (xc_name.find(':') != std::string::npos) { const std::string pair = canonicalize_libxc_pair(xc_name); - if (!pair.empty()) - { + if (!pair.empty()) { canonical = pair; } } - if (canonical.empty()) - { + 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 // pseudopotential metadata. Keep that established spelling at the adapter // boundary while using s-dftd3's canonical method identifier internally. - if (canonical == "hse") - { + if (canonical == "hse") { canonical = "hse06"; } return canonical; @@ -139,8 +125,7 @@ std::string canonicalize_method_name(const std::string& input) bool lookup_parameters(const std::string& method, Damping damping, Parameters& parameters, - std::string& canonical_method) -{ + std::string& canonical_method) { canonical_method = canonicalize_method_name(method); // Preserve the two historical LibXC spellings whose D3 parameter set is @@ -149,33 +134,25 @@ bool lookup_parameters(const std::string& method, // lookup key for the damping variant where the legacy ABACUS mapping was // defined. std::string lookup_method = canonical_method; - if (damping == Damping::Rational && lookup_method == "wb97xv") - { + if (damping == Damping::Rational && lookup_method == "wb97xv") { lookup_method = "wb97x"; - } - else if (damping == Damping::Zero && lookup_method == "wb97xd3") - { + } else if (damping == Damping::Zero && lookup_method == "wb97xd3") { lookup_method = "wb97x"; } const char* method_id = nullptr; - for (const MethodAlias& alias : kMethodAliases) - { - if (normalized(alias.alias) == lookup_method) - { + for (const MethodAlias& alias: kMethodAliases) { + if (normalized(alias.alias) == lookup_method) { method_id = alias.method_id; break; } } - if (method_id == nullptr) - { + if (method_id == nullptr) { return false; } - for (const DampingParameterRecord& record : kDampingParameters) - { - if (record.damping == damping && std::string(record.method_id) == method_id) - { + for (const DampingParameterRecord& record: kDampingParameters) { + if (record.damping == damping && std::string(record.method_id) == method_id) { parameters.damping = damping; parameters.s6 = record.s6; parameters.s8 = record.s8; diff --git a/source/source_hamilt/module_vdw/vdwd3_parameters.h b/source/source_hamilt/module_vdw/vdwd3_parameters.h index 5080180b6cb..f04777ce885 100644 --- a/source/source_hamilt/module_vdw/vdwd3_parameters.h +++ b/source/source_hamilt/module_vdw/vdwd3_parameters.h @@ -5,10 +5,8 @@ #include -namespace vdw -{ -namespace d3 -{ +namespace vdw { +namespace d3 { std::string canonicalize_method_name(const std::string& input); diff --git a/source/source_hamilt/module_vdw/vdwd3_types.h b/source/source_hamilt/module_vdw/vdwd3_types.h index 803ec8ca69d..611371734b8 100644 --- a/source/source_hamilt/module_vdw/vdwd3_types.h +++ b/source/source_hamilt/module_vdw/vdwd3_types.h @@ -4,13 +4,10 @@ #include #include -namespace vdw -{ -namespace d3 -{ +namespace vdw { +namespace d3 { -struct Vec3 -{ +struct Vec3 { double x = 0.0; double y = 0.0; double z = 0.0; @@ -18,16 +15,14 @@ struct Vec3 Vec3() = default; Vec3(double x_in, double y_in, double z_in) : x(x_in), y(y_in), z(z_in) {} - Vec3& operator+=(const Vec3& rhs) - { + Vec3& operator+=(const Vec3& rhs) { x += rhs.x; y += rhs.y; z += rhs.z; return *this; } - Vec3& operator-=(const Vec3& rhs) - { + Vec3& operator-=(const Vec3& rhs) { x -= rhs.x; y -= rhs.y; z -= rhs.z; @@ -35,59 +30,40 @@ struct Vec3 } }; -inline Vec3 operator+(Vec3 lhs, const Vec3& rhs) -{ +inline Vec3 operator+(Vec3 lhs, const Vec3& rhs) { lhs += rhs; return lhs; } -inline Vec3 operator-(Vec3 lhs, const Vec3& rhs) -{ +inline Vec3 operator-(Vec3 lhs, const Vec3& rhs) { lhs -= rhs; return lhs; } -inline Vec3 operator-(const Vec3& value) -{ - return Vec3(-value.x, -value.y, -value.z); -} +inline Vec3 operator-(const Vec3& value) { return Vec3(-value.x, -value.y, -value.z); } -inline Vec3 operator*(const Vec3& value, double factor) -{ +inline Vec3 operator*(const Vec3& value, double factor) { return Vec3(value.x * factor, value.y * factor, value.z * factor); } -inline Vec3 operator*(double factor, const Vec3& value) -{ - return value * factor; -} +inline Vec3 operator*(double factor, const Vec3& value) { return value * factor; } -inline Vec3 operator/(const Vec3& value, double divisor) -{ - return value * (1.0 / divisor); -} +inline Vec3 operator/(const Vec3& value, double divisor) { return value * (1.0 / divisor); } -struct Matrix3 -{ +struct Matrix3 { double value[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; }; -struct Structure -{ +struct Structure { std::vector atomic_numbers; std::vector positions; std::array lattice; std::array periodic = {{false, false, false}}; }; -enum class Damping -{ - Zero, - Rational -}; +enum class Damping { Zero, Rational }; -struct Parameters -{ +struct Parameters { Damping damping = Damping::Rational; double s6 = 1.0; double s8 = 1.0; @@ -99,8 +75,7 @@ struct Parameters double alp = 14.0; }; -struct Cutoffs -{ +struct Cutoffs { double disp2 = 60.0; double disp3 = 40.0; double cn = 40.0; @@ -108,8 +83,7 @@ struct Cutoffs double width3 = 0.0; }; -struct Result -{ +struct Result { double energy = 0.0; std::vector gradient; Matrix3 virial; diff --git a/source/source_hamilt/module_vdw/vdwd4.cpp b/source/source_hamilt/module_vdw/vdwd4.cpp index ad4a93e9925..f95d9c002f6 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" @@ -22,24 +23,19 @@ #include #endif -namespace vdw -{ +namespace vdw { -namespace -{ +namespace { -std::string to_lower(std::string value) -{ +std::string to_lower(std::string value) { std::transform(value.begin(), value.end(), value.begin(), [](unsigned char c) { return static_cast(std::tolower(c)); }); return value; } -void check_dftd4_error(dftd4_error error, const std::string& where) -{ - if (dftd4_check_error(error)) - { +void check_dftd4_error(dftd4_error error, const std::string& where) { + if (dftd4_check_error(error)) { char buffer[1024]; std::memset(buffer, 0, sizeof(buffer)); int buffersize = static_cast(sizeof(buffer)); @@ -49,12 +45,9 @@ void check_dftd4_error(dftd4_error error, const std::string& where) } } -int atomic_number_from_symbol(const std::string& symbol) -{ - for (int i = 0; i < static_cast(ModuleBase::element_name.size()); ++i) - { - if (symbol == ModuleBase::element_name[i]) - { +int atomic_number_from_symbol(const std::string& symbol) { + for (int i = 0; i < static_cast(ModuleBase::element_name.size()); ++i) { + if (symbol == ModuleBase::element_name[i]) { return i + 1; // DFT-D4 expects the true atomic number. } } @@ -63,14 +56,11 @@ int atomic_number_from_symbol(const std::string& symbol) return 0; } -double length_to_bohr(double value, const std::string& unit) -{ - if (unit == "A") - { +double length_to_bohr(double value, const std::string& unit) { + if (unit == "A") { return value / ModuleBase::BOHR_TO_A; } - if (unit == "Bohr") - { + if (unit == "Bohr") { return value; } @@ -78,38 +68,30 @@ double length_to_bohr(double value, const std::string& unit) return value; } -double cutoff_to_bohr(const std::string& value, const std::string& unit) -{ +double cutoff_to_bohr(const std::string& value, const std::string& unit) { return length_to_bohr(std::stod(value), unit); } } // namespace Vdwd4::Vdwd4(const UnitCell& unit_in, const std::string& xc_name, const Input_para& input) - : Vdw(unit_in), xc_name_(normalize_xc_name(xc_name)), model_name_(to_lower(input.vdw_d4_model)) -{ + : Vdw(unit_in), xc_name_(normalize_xc_name(xc_name)), model_name_(to_lower(input.vdw_d4_model)) { cutoff_disp2_ = cutoff_to_bohr(input.vdw_cutoff_radius, input.vdw_radius_unit); cutoff_disp3_ = std::min(40.0, cutoff_disp2_); cutoff_cn_ = length_to_bohr(input.vdw_cn_thr, input.vdw_cn_thr_unit); smooth_width_2b_ = input.vdw_cutoff_width2; smooth_width_3b_ = input.vdw_cutoff_width3; - 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"); + 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"); } - 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"); + 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"); } double valence_charge = 0.0; - for (int it = 0; it < ucell_.ntype; ++it) - { - if (ucell_.atoms[it].flag_empty_element) - { + for (int it = 0; it < ucell_.ntype; ++it) { + if (ucell_.atoms[it].flag_empty_element) { continue; } valence_charge += ucell_.atoms[it].ncpp.zv * ucell_.atoms[it].na; @@ -121,8 +103,7 @@ void Vdwd4::build_structure(std::vector& numbers, std::vector& positions, std::vector& lattice, std::array& periodic, - std::vector& atom_indices) const -{ + std::vector& atom_indices) const { numbers.clear(); positions.clear(); lattice.clear(); @@ -134,18 +115,15 @@ void Vdwd4::build_structure(std::vector& numbers, atom_indices.reserve(ucell_.nat); int iat = 0; - for (int it = 0; it < ucell_.ntype; ++it) - { - if (ucell_.atoms[it].flag_empty_element) - { + for (int it = 0; it < ucell_.ntype; ++it) { + if (ucell_.atoms[it].flag_empty_element) { iat += ucell_.atoms[it].na; continue; } const int atomic_number = atomic_number_from_symbol(ucell_.atoms[it].ncpp.psd); - for (int ia = 0; ia < ucell_.atoms[it].na; ++ia) - { + for (int ia = 0; ia < ucell_.atoms[it].na; ++ia) { const ModuleBase::Vector3 position = ucell_.atoms[it].tau[ia] * ucell_.lat0; numbers.push_back(atomic_number); @@ -180,8 +158,7 @@ void Vdwd4::build_structure(std::vector& numbers, void Vdwd4::compute(double& energy_ha, std::vector* gradient_ha_bohr, std::array* sigma_ha, - std::vector& atom_indices) -{ + std::vector& atom_indices) { #ifdef __DFTD4 std::vector numbers; std::vector positions; @@ -190,11 +167,8 @@ 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())) - { - ModuleBase::WARNING_QUIT("Vdwd4::compute", - "gradient_ha_bohr is too small for the filtered atom set."); + 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."); } // These vectors own all arrays passed to DFT-D4. Their data() pointers @@ -211,18 +185,13 @@ void Vdwd4::compute(double& energy_ha, check_dftd4_error(error, "dftd4_new_structure"); dftd4_model model = nullptr; - if (model_name_ == "d4") - { + if (model_name_ == "d4") { model = dftd4_new_d4_model(error, mol); check_dftd4_error(error, "dftd4_new_d4_model"); - } - else if (model_name_ == "d4s") - { + } else if (model_name_ == "d4s") { model = dftd4_new_d4s_model(error, mol); check_dftd4_error(error, "dftd4_new_d4s_model"); - } - else - { + } else { ModuleBase::WARNING_QUIT("Vdwd4::compute", "Unsupported DFT-D4 model: " + model_name_); } @@ -248,7 +217,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"); @@ -262,12 +232,10 @@ void Vdwd4::compute(double& energy_ha, #ifdef __MPI Parallel_Reduce::reduce_all(energy_ha); - if (gradient != nullptr) - { + if (gradient != nullptr) { Parallel_Reduce::reduce_all(gradient, 3 * ucell_.nat); } - if (sigma != nullptr) - { + if (sigma != nullptr) { Parallel_Reduce::reduce_all(sigma, 9); } #endif @@ -283,12 +251,10 @@ void Vdwd4::compute(double& energy_ha, void Vdwd4::set_force_from_gradient(const std::vector& gradient_ha_bohr, const std::vector& atom_indices, - VdwResult& result) const -{ + VdwResult& result) const { result.force.assign(ucell_.nat, ModuleBase::Vector3(0.0, 0.0, 0.0)); - for (std::size_t index = 0; index < atom_indices.size(); ++index) - { + for (std::size_t index = 0; index < atom_indices.size(); ++index) { const int iat = atom_indices[index]; // DFT-D4 returns dE/dR in Ha/Bohr; ABACUS forces are -dE/dR in Ry/Bohr. result.force[iat].x = -2.0 * gradient_ha_bohr[3 * index + 0]; @@ -299,11 +265,10 @@ 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], @@ -317,15 +282,13 @@ void Vdwd4::set_stress_from_sigma(const std::array& sigma_ha, result.has_stress = true; } -void Vdwd4::evaluate_impl(const VdwRequest& request, VdwResult& result) -{ +void Vdwd4::evaluate_impl(const VdwRequest& request, VdwResult& result) { ModuleBase::TITLE("Vdwd4", "evaluate"); ModuleBase::timer::start("Vdwd4", "evaluate"); double energy_ha = 0.0; std::vector atom_indices; - if (request.force || request.stress) - { + if (request.force || request.stress) { std::vector gradient(3 * ucell_.nat, 0.0); std::array sigma; sigma.fill(0.0); @@ -334,17 +297,13 @@ void Vdwd4::evaluate_impl(const VdwRequest& request, VdwResult& result) // Keep all requested quantities from this single call. compute(energy_ha, &gradient, &sigma, atom_indices); - if (request.force) - { + if (request.force) { set_force_from_gradient(gradient, atom_indices, result); } - if (request.stress) - { + if (request.stress) { set_stress_from_sigma(sigma, result); } - } - else - { + } else { compute(energy_ha, nullptr, nullptr, atom_indices); } diff --git a/source/source_hamilt/module_vdw/vdwd4.h b/source/source_hamilt/module_vdw/vdwd4.h index 1c1e580e3fe..a861a5a6e65 100644 --- a/source/source_hamilt/module_vdw/vdwd4.h +++ b/source/source_hamilt/module_vdw/vdwd4.h @@ -7,11 +7,9 @@ #include #include -namespace vdw -{ +namespace vdw { -class Vdwd4 : public Vdw -{ +class Vdwd4 : public Vdw { public: Vdwd4(const UnitCell& unit_in, const std::string& xc_name, const Input_para& input); @@ -20,18 +18,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, diff --git a/tests/09_DeePKS/check_file.cpp b/tests/09_DeePKS/check_file.cpp index d7e1391c9e0..6022200d265 100644 --- a/tests/09_DeePKS/check_file.cpp +++ b/tests/09_DeePKS/check_file.cpp @@ -1,14 +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,37 +19,32 @@ 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()) - { + while (ref.good()) { ref >> word_ref; num_ref++; } - while (cal.good()) - { + while (cal.good()) { cal >> word_cal; 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; return 0; @@ -66,32 +60,26 @@ 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; - } - else if(fabs(data_ref-data_cal) > threshold) - { + cout << "Data No." << index + 1 << " cal=" << word_cal << " ref=" << word_ref << endl; + } 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..6022200d265 100644 --- a/tests/integrate/check_file.cpp +++ b/tests/integrate/check_file.cpp @@ -1,14 +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,37 +19,32 @@ 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()) - { + while (ref.good()) { ref >> word_ref; num_ref++; } - while (cal.good()) - { + while (cal.good()) { cal >> word_cal; 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; return 0; @@ -66,32 +60,26 @@ 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; - } - else if(fabs(data_ref-data_cal) > threshold) - { + cout << "Data No." << index + 1 << " cal=" << word_cal << " ref=" << word_ref << endl; + } 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 +}