diff --git a/CMakeLists.txt b/CMakeLists.txt index 8c6fd2078..d60c8966c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -171,7 +171,7 @@ if(Git_FOUND AND EXISTS "${PROJECT_SOURCE_DIR}/.git") # Update submodules as needed option(GIT_SUBMODULE "Check submodules during build" ON) - if(GIT_SUBMODULE) + if(GIT_SUBMODULE AND NOT EXISTS "${PROJECT_SOURCE_DIR}/external/mstd/include/") message(STATUS "Submodule update") execute_process(COMMAND ${GIT_EXECUTABLE} submodule update --init --recursive WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} diff --git a/benchmarks/perf/perfBenchSetup.hpp b/benchmarks/perf/perfBenchSetup.hpp index ac59e8bb8..82f42efbb 100644 --- a/benchmarks/perf/perfBenchSetup.hpp +++ b/benchmarks/perf/perfBenchSetup.hpp @@ -38,6 +38,7 @@ #include "matrix.hpp" #include "molecule.hpp" #include "simulationBox.hpp" +#include "strongTypes.hpp" namespace potential { @@ -109,7 +110,7 @@ namespace benchSetup atom->setShiftForce({0.0, 0.0, 0.0}); atom->setMass(12.0); atom->setAtomType(i % 2); - atom->setInternalGlobalVDWType(i % 2); + atom->setInternalGlobalVDWType(VdwType{i % 2}); atom->setPartialCharge((i % 2 == 0) ? 0.4 : -0.4); molecule.addAtom(atom); @@ -131,7 +132,12 @@ namespace benchSetup ) ); - auto pair = potential::LennardJonesPair(0UL, 1UL, 12.0, 2.0, 3.0); + auto pair = potential::LennardJonesPair( + ExtVdwType(0), + ExtVdwType(1), + 12.0, + LJParams{.c6 = 2.0, .c12 = 3.0} + ); potential.setNonCoulombPairsMatrix(0, 1, pair); potential.setNonCoulombPairsMatrix(1, 0, pair); diff --git a/benchmarks/perf/perfInterWater.cpp b/benchmarks/perf/perfInterWater.cpp index 9e384e8e4..c58597fca 100644 --- a/benchmarks/perf/perfInterWater.cpp +++ b/benchmarks/perf/perfInterWater.cpp @@ -29,6 +29,8 @@ #include #include +#include "strongTypes.hpp" + #ifdef PQ_WITH_CALLGRIND #include #else @@ -90,7 +92,7 @@ int main() atom->setAtomicNumber(atomicNumber); atom->setPosition(pos); atom->setAtomType(0); - atom->setInternalGlobalVDWType(0); + atom->setInternalGlobalVDWType(VdwType{0}); atom->setPartialCharge(charge); atom->setForceToZero(); return atom; @@ -130,14 +132,20 @@ int main() } InterWaterState state; - state._oxygenCharge = -0.82; - state._hydrogenCharge = 0.41; - state._nonCoulombPairOO = - std::make_unique(CUTOFF, -2.0, 4.0); - state._nonCoulombPairOH = - std::make_unique(CUTOFF, -0.5, 1.5); - state._nonCoulombPairHH = - std::make_unique(CUTOFF, -0.2, 0.8); + state._oxygenCharge = -0.82; + state._hydrogenCharge = 0.41; + state._nonCoulombPairOO = std::make_unique( + CUTOFF, + LJParams{.c6 = 2.0, .c12 = 4.0} + ); + state._nonCoulombPairOH = std::make_unique( + CUTOFF, + LJParams{.c6 = 0.5, .c12 = 1.5} + ); + state._nonCoulombPairHH = std::make_unique( + CUTOFF, + LJParams{.c6 = 0.2, .c12 = 0.8} + ); InterWater interWater( std::move(state), diff --git a/benchmarks/perf/perfNonCoulombPairs.cpp b/benchmarks/perf/perfNonCoulombPairs.cpp index 0686ce303..3f860c601 100644 --- a/benchmarks/perf/perfNonCoulombPairs.cpp +++ b/benchmarks/perf/perfNonCoulombPairs.cpp @@ -42,9 +42,19 @@ static constexpr std::uint64_t ITERATIONS = 20000; int main() { - auto lj = potential::LennardJonesPair(9.0, 2.0, 3.0); - auto buck = potential::BuckinghamPair(9.0, 1.0, 0.3, 2.0); - auto morse = potential::MorsePair(9.0, 1.0, 2.0, 1.5); + auto lj = potential::LennardJonesPair(9.0, LJParams{.c6 = 2.0, .c12 = 3.0}); + auto buck = potential::BuckinghamPair( + 9.0, + BuckinghamParams{.scaling = 1.0, .dRho = 0.3, .c6 = 2.0} + ); + auto morse = potential::MorsePair( + 9.0, + MorseParams{ + .dissociationEnergy = 1.0, + .wellWidth = 2.0, + .equilibriumDistance = 1.5 + } + ); CALLGRIND_ZERO_STATS; diff --git a/benchmarks/src/benchmarkSetup.hpp b/benchmarks/src/benchmarkSetup.hpp index 398eb8110..c471a6ba7 100644 --- a/benchmarks/src/benchmarkSetup.hpp +++ b/benchmarks/src/benchmarkSetup.hpp @@ -28,6 +28,7 @@ #include "atom.hpp" #include "molecule.hpp" #include "simulationBox.hpp" +#include "strongTypes.hpp" #include "vector3d.hpp" namespace benchmarkSetup @@ -68,7 +69,7 @@ namespace benchmarkSetup atom->setForce({0.1, -0.2, 0.05}); atom->setMass(12.0); atom->setAtomType(0); - atom->setInternalGlobalVDWType(0); + atom->setInternalGlobalVDWType(VdwType{0}); atom->setPartialCharge(atomIndex++ % 2 == 0 ? 0.4 : -0.4); atom->setShiftForce({0.0, 0.0, 0.0}); diff --git a/benchmarks/src/potential/benchmarkForceCalculation.cpp b/benchmarks/src/potential/benchmarkForceCalculation.cpp index fb902f406..fab3bf611 100644 --- a/benchmarks/src/potential/benchmarkForceCalculation.cpp +++ b/benchmarks/src/potential/benchmarkForceCalculation.cpp @@ -50,8 +50,7 @@ namespace const auto pair = std::make_shared( benchmarkSetup::cutOff, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ); nonCoulomb->setGuffNonCoulPair({1, 1, 0, 0}, pair); diff --git a/benchmarks/src/potential/benchmarkPairPotentials.cpp b/benchmarks/src/potential/benchmarkPairPotentials.cpp index be9cd60ef..705b8b841 100644 --- a/benchmarks/src/potential/benchmarkPairPotentials.cpp +++ b/benchmarks/src/potential/benchmarkPairPotentials.cpp @@ -67,19 +67,32 @@ namespace void BM_LennardJones(benchmark::State& state) { - potential::LennardJonesPair potential(9.0, 2.0, 3.0); + potential::LennardJonesPair potential( + 9.0, + LJParams{.c6 = 2.0, .c12 = 3.0} + ); runNonCoulombBenchmark(state, potential); } void BM_Buckingham(benchmark::State& state) { - potential::BuckinghamPair potential(9.0, 1.0, 0.3, 2.0); + potential::BuckinghamPair potential( + 9.0, + BuckinghamParams{.scaling = 1.0, .dRho = 0.3, .c6 = 2.0} + ); runNonCoulombBenchmark(state, potential); } void BM_Morse(benchmark::State& state) { - potential::MorsePair potential(9.0, 1.0, 2.0, 1.5); + potential::MorsePair potential( + 9.0, + MorseParams{ + .dissociationEnergy = 1.0, + .wellWidth = 2.0, + .equilibriumDistance = 1.5 + } + ); runNonCoulombBenchmark(state, potential); } diff --git a/changes/developer/enhancement.strong-types.md b/changes/developer/enhancement.strong-types.md index 720b9acbf..f55ed6f01 100644 --- a/changes/developer/enhancement.strong-types.md +++ b/changes/developer/enhancement.strong-types.md @@ -2,3 +2,8 @@ - add strong type for `BondId` - add strong type for `AngleId` - add strong type for `DihedralId` +- add strong types for `ExtVdwType` and `VdwType` +- add strong type for morse params aka `MorseParams` +- add strong type for lennard jones params aka `LJParams` +- add strong type for bukcingham params aka `BuckinghamParams` +- make guff pari used `std::array` instead of `std::vector` diff --git a/external/mstd b/external/mstd index 2e2393c74..a197a764d 160000 --- a/external/mstd +++ b/external/mstd @@ -1 +1 @@ -Subproject commit 2e2393c74a5d33babe9242ec2e68430ce7c6ebb4 +Subproject commit a197a764da7b7f98d5683c226554566ffae7ec8d diff --git a/include/config/defaults.hpp b/include/config/defaults.hpp index 43db9003b..5ded2f22d 100644 --- a/include/config/defaults.hpp +++ b/include/config/defaults.hpp @@ -131,6 +131,8 @@ namespace defaults static constexpr double ABS_ENERGY_CONV_DEFAULT = 1.0e-6; static constexpr double MAX_FORCE_CONV_DEFAULT = 1.0e-6; static constexpr double RMS_FORCE_CONV_DEFAULT = 1.0e-6; + + static constexpr auto NUM_GUFF_COEFFICIENTS = 22; // clang-format on } // namespace defaults diff --git a/include/config/guffCoefficients.hpp b/include/config/guffCoefficients.hpp index 4585b4548..e93cd051e 100644 --- a/include/config/guffCoefficients.hpp +++ b/include/config/guffCoefficients.hpp @@ -24,39 +24,51 @@ #define _GUFF_COEFFICIENTS_HPP_ +#include "strongTypes.hpp" + namespace constants { // clang-format off // values calculated from the original SPC publication static constexpr double SPC_LJ_C6_OO = -625.455653639347; // A^6 kcal mol^-1 static constexpr double SPC_LJ_C12_OO = 629358.472583307; // A^12 kcal mol^-1 + static constexpr LJParams SPC_LJ_PARAMS_OO{.c6 = SPC_LJ_C6_OO, .c12 = SPC_LJ_C12_OO}; static constexpr double SPC_E_LJ_C6_OO = SPC_LJ_C6_OO; // A^6 kcal mol^-1 static constexpr double SPC_E_LJ_C12_OO = SPC_LJ_C12_OO; // A^12 kcal mol^-1 + static constexpr LJParams SPC_E_LJ_PARAMS_OO{.c6 = SPC_E_LJ_C6_OO, .c12 = SPC_E_LJ_C12_OO}; static constexpr double SPC_FW_LJ_C6_OO = -625.5024676571352; // A^6 kcal mol^-1 static constexpr double SPC_FW_LJ_C12_OO = 629326.9774051674; // A^12 kcal mol^-1 + static constexpr LJParams SPC_FW_LJ_PARAMS_OO{.c6 = SPC_FW_LJ_C6_OO, .c12 = SPC_FW_LJ_C12_OO}; static constexpr double QSPC_FW_LJ_C6_OO = -625.5020652114152; // A^6 kcal mol^-1 static constexpr double QSPC_FW_LJ_C12_OO = 629326.5724987736; // A^12 kcal mol^-1 + static constexpr LJParams QSPC_FW_LJ_PARAMS_OO{.c6 = QSPC_FW_LJ_C6_OO, .c12 = QSPC_FW_LJ_C12_OO}; static constexpr double SPC_DC_LJ_C6_OO = -779.8414665600154; // A^6 kcal mol^-1 static constexpr double SPC_DC_LJ_C12_OO = 773048.5510230307; // A^12 kcal mol^-1 + static constexpr LJParams SPC_DC_LJ_PARAMS_OO{.c6 = SPC_DC_LJ_C6_OO, .c12 = SPC_DC_LJ_C12_OO}; static constexpr double H2O_DC_LJ_C6_OO = -590.6923729027751; // A^6 kcal mol^-1 static constexpr double H2O_DC_LJ_C12_OO = 615459.8371975797; // A^12 kcal mol^-1 + static constexpr LJParams H2O_DC_LJ_PARAMS_OO{.c6 = H2O_DC_LJ_C6_OO, .c12 = H2O_DC_LJ_C12_OO}; static constexpr double TIP3P_LJ_C6_OO = -595.067688427684; // A^6 kcal mol^-1 static constexpr double TIP3P_LJ_C12_OO = 582015.099443679; // A^12 kcal mol^-1 + static constexpr LJParams TIP3P_LJ_PARAMS_OO{.c6 = TIP3P_LJ_C6_OO, .c12 = TIP3P_LJ_C12_OO}; static constexpr double OPC3_LJ_C6_OO = -668.637501352773; // A^6 kcal mol^-1 static constexpr double OPC3_LJ_C12_OO = 683996.561589584; // A^12 kcal mol^-1 + static constexpr LJParams OPC3_LJ_PARAMS_OO{.c6 = OPC3_LJ_C6_OO, .c12 = OPC3_LJ_C12_OO}; static constexpr double SPC_MTR_LJ_C6_OO = -613.527724665392; // A^6 kcal mol^-1 static constexpr double SPC_MTR_LJ_C12_OO = 629302.103250478; // A^12 kcal mol^-1 + static constexpr LJParams SPC_MTR_LJ_PARAMS_OO{.c6 = SPC_MTR_LJ_C6_OO, .c12 = SPC_MTR_LJ_C12_OO}; static constexpr double TIP3P_MTR_LJ_C6_OO = -454.1108986615679; // A^6 kcal mol^-1 static constexpr double TIP3P_MTR_LJ_C12_OO = 513862.3326959847; // A^12 kcal mol^-1 + static constexpr LJParams TIP3P_MTR_LJ_PARAMS_OO{.c6 = TIP3P_MTR_LJ_C6_OO, .c12 = TIP3P_MTR_LJ_C12_OO}; // clang-format on } // namespace constants diff --git a/include/forceField/forcefield.tpp.hpp b/include/forceField/forcefield.tpp.hpp index 946ee54c7..a0b69e8bb 100644 --- a/include/forceField/forcefield.tpp.hpp +++ b/include/forceField/forcefield.tpp.hpp @@ -88,10 +88,10 @@ namespace forceField const auto vdwType1 = molecule1->getInternalGlobalVDWType(atomIndex1); const auto vdwType2 = molecule2->getInternalGlobalVDWType(atomIndex2); - const auto indices = - {molType1, molType2, atomType1, atomType2, vdwType1, vdwType2}; + const auto indices = {molType1, molType2, atomType1, atomType2}; - const auto nonCoulombPair = nonCoulPot.getNonCoulPair(indices); + const auto nonCoulombPair = + nonCoulPot.getNonCoulPair(indices, {vdwType1, vdwType2}); if (distance < nonCoulombPair->getRadialCutOff()) { diff --git a/include/input/guffDatReader.hpp b/include/input/guffDatReader.hpp index 9ec0890f6..708c47066 100644 --- a/include/input/guffDatReader.hpp +++ b/include/input/guffDatReader.hpp @@ -106,12 +106,13 @@ namespace input::guffdat const double rncCutOff ); void addGuffPair( - const size_t molType1, - const size_t molType2, - const size_t atomType1, - const size_t atomType2, - const std::vector &coefficients, - const double rncCutOff + const size_t molType1, + const size_t molType2, + const size_t atomType1, + const size_t atomType2, + const std::array + &coefficients, + const double rncCutOff ); /******************** diff --git a/include/molsys/atom.hpp b/include/molsys/atom.hpp index 63258b990..803066c3d 100644 --- a/include/molsys/atom.hpp +++ b/include/molsys/atom.hpp @@ -49,8 +49,8 @@ namespace molsys std::string _name; std::string _atomTypeName; - size_t _externalGlobalVDWType; - size_t _internalGlobalVDWType; + ExtVdwType _externalGlobalVDWType; + VdwType _internalGlobalVDWType; size_t _externalAtomType; size_t _atomType; @@ -125,8 +125,8 @@ namespace molsys [[nodiscard]] size_t getExternalAtomType() const; [[nodiscard]] size_t getAtomType() const; - [[nodiscard]] size_t getExternalGlobalVDWType() const; - [[nodiscard]] size_t getInternalGlobalVDWType() const; + [[nodiscard]] ExtVdwType getExternalGlobalVDWType() const; + [[nodiscard]] VdwType getInternalGlobalVDWType() const; [[nodiscard]] AtomNumber getAtomicNumber() const @@ -166,8 +166,8 @@ namespace molsys void setAtomType(const size_t atomType); void setExternalAtomType(const size_t externalAtomType); - void setExternalGlobalVDWType(const size_t externalGlobalVDWType); - void setInternalGlobalVDWType(const size_t internalGlobalVDWType); + void setExternalGlobalVDWType(const ExtVdwType externalGlobalVDWType); + void setInternalGlobalVDWType(const VdwType internalGlobalVDWType); void setPosition(const linearAlgebra::Vec3D &position); void setVelocity(const linearAlgebra::Vec3D &velocity); diff --git a/include/molsys/molecule.hpp b/include/molsys/molecule.hpp index 12b7a8ac6..859a98ec2 100644 --- a/include/molsys/molecule.hpp +++ b/include/molsys/molecule.hpp @@ -32,6 +32,7 @@ #include // for vector #include "atom.hpp" // for Atom +#include "strongTypes.hpp" namespace molsys { @@ -109,8 +110,8 @@ namespace molsys void scale(const linearAlgebra::tensor3D &, const Box &); void scaleVelocity(const linearAlgebra::tensor3D &, const Box &); - [[nodiscard]] size_t getNumberOfAtomTypes(); - [[nodiscard]] std::vector getExternalGlobalVDWTypes() const; + [[nodiscard]] size_t getNumberOfAtomTypes(); + [[nodiscard]] std::vector getExternalGlobalVDWTypes() const; [[nodiscard]] std::vector getAtomMasses() const; [[nodiscard]] std::vector getPartialCharges() const; @@ -184,12 +185,13 @@ namespace molsys const size_t index ) const; - [[nodiscard]] AtomNumber getAtomicNumber(const size_t index) const; - [[nodiscard]] double getAtomMass(const size_t index) const; - [[nodiscard]] double getPartialCharge(const size_t index) const; - [[nodiscard]] size_t getAtomType(const size_t index) const; - [[nodiscard]] size_t getInternalGlobalVDWType(const size_t index) const; + [[nodiscard]] AtomNumber getAtomicNumber(const size_t index) const; + [[nodiscard]] double getAtomMass(const size_t index) const; + [[nodiscard]] double getPartialCharge(const size_t index) const; + [[nodiscard]] size_t getAtomType(const size_t index) const; [[nodiscard]] std::string getAtomName(const size_t index) const; + [[nodiscard]] + VdwType getInternalGlobalVDWType(const size_t index) const; /*************************** * standard getter methods * diff --git a/include/molsys/moleculeType.hpp b/include/molsys/moleculeType.hpp index 2754cebb0..f32d61e5e 100644 --- a/include/molsys/moleculeType.hpp +++ b/include/molsys/moleculeType.hpp @@ -29,6 +29,8 @@ #include // for string_view #include // for vector +#include "strongTypes.hpp" + namespace molsys { /** @@ -49,7 +51,7 @@ namespace molsys std::vector _atomNames; std::vector _atomTypes; std::vector _externalAtomTypes; - std::vector _externalGlobalVDWTypes; + std::vector _externalGlobalVDWTypes; std::vector _partialCharges; std::map _externalToInternalAtomTypes; @@ -68,7 +70,7 @@ namespace molsys void addAtomName(const std::string &atomName); void addExternalAtomType(const size_t externalAtomType); void addPartialCharge(const double partialCharge); - void addExternalGlobalVDWType(const size_t externalGlobalVDWType); + void addExternalGlobalVDWType(const ExtVdwType externalGlobalVDWType); void addExternalToInternalAtomTypeElement(const size_t, const size_t); void addAtomType(const size_t atomType); @@ -104,7 +106,7 @@ namespace molsys [[nodiscard]] std::vector getAtomNames() const; [[nodiscard]] std::vector &getExternalAtomTypes(); - [[nodiscard]] std::vector &getExternalGlobalVDWTypes(); + [[nodiscard]] std::vector &getExternalGlobalVDWTypes(); [[nodiscard]] std::vector &getPartialCharges(); [[nodiscard]] std::map getExternalToInternalAtomTypes( diff --git a/include/molsys/simulationBox.hpp b/include/molsys/simulationBox.hpp index b4d78a8de..571e5879e 100644 --- a/include/molsys/simulationBox.hpp +++ b/include/molsys/simulationBox.hpp @@ -24,11 +24,11 @@ #define _SIMULATION_BOX_HPP_ -#include // for map #include #include // for optional #include #include // for string +#include #include // for vector #include "atom.hpp" // for Atom @@ -37,6 +37,7 @@ #include "moleculeType.hpp" // for MoleculeType #include "orthorhombicBox.hpp" // for OrthorhombicBox #include "simulationBoxView.hpp" // for SimulationBoxView +#include "strongTypes.hpp" /** * @namespace molsys @@ -86,8 +87,9 @@ namespace molsys std::vector _molecules; std::vector _moleculeTypes; - std::vector _externalGlobalVdwTypes; - std::map _externalToInternalGlobalVDWTypes; + std::vector _externalGlobalVdwTypes; + std::unordered_map + _externalToInternalGlobalVDWTypes; public: void copy(const SimulationBox&); @@ -206,9 +208,12 @@ namespace molsys [[nodiscard]] const std::vector& getMolecules() const; [[nodiscard]] std::vector& getMoleculeTypes(); - [[nodiscard]] std::vector& getExternalGlobalVdwTypes(); - [[nodiscard]] std::map& getExternalToInternalGlobalVDWTypes( - ); + [[nodiscard]] + std::vector& getExternalGlobalVdwTypes(); + [[nodiscard]] + std::unordered_map< + ExtVdwType, + VdwType>& getExternalToInternalGlobalVDWTypes(); [[nodiscard]] Box& getBox(); [[nodiscard]] Box& getBox() const; diff --git a/include/potential/nonCoulomb/buckinghamPair.hpp b/include/potential/nonCoulomb/buckinghamPair.hpp index 85057fb73..44acb1573 100644 --- a/include/potential/nonCoulomb/buckinghamPair.hpp +++ b/include/potential/nonCoulomb/buckinghamPair.hpp @@ -24,10 +24,12 @@ #define _BUCKINGHAM_PAIR_HPP_ -#include // size_t #include // pair #include "nonCoulombPair.hpp" +#include "strongTypes.hpp" + +struct TestBuckinghamPairUtils; // forward declaration namespace potential { @@ -40,46 +42,46 @@ namespace potential class BuckinghamPair : public NonCoulombPair { private: - double _a; - double _dRho; - double _c6; + BuckinghamParams _params; public: explicit BuckinghamPair( - const size_t vanDerWaalsType1, - const size_t vanDerWaalsType2, - const double cutOff, - const double a, - const double dRho, - const double c6 + const ExtVdwType vanDerWaalsType1, + const ExtVdwType vanDerWaalsType2, + const double cutOff, + const BuckinghamParams& params ); explicit BuckinghamPair( - const double cutOff, - const double a, - const double dRho, - const double c6 + const double cutOff, + const BuckinghamParams& params ); explicit BuckinghamPair( - const double cutOff, - const double energyCutoff, - const double forceCutoff, - const double a, - const double dRho, - const double c6 + const double cutOff, + const double energyCutoff, + const double forceCutoff, + const BuckinghamParams& params ); - [[nodiscard]] bool operator==(const BuckinghamPair &other) const; + // TODO: we need to explicitly delete it to not implicitly create it + // with the wrong types!!! Needs cleanup + explicit BuckinghamPair( + const size_t, + const size_t, + const double, + const BuckinghamParams& params + ) = delete; + + [[nodiscard]] bool operator==(const BuckinghamPair& other) const; - [[nodiscard]] std::pair calculate(const double distance + [[nodiscard]] std::pair calculate( + const double distance ) const override; - [[nodiscard]] double getA() const; - [[nodiscard]] double getDRho() const; - [[nodiscard]] double getC6() const; + friend struct ::TestBuckinghamPairUtils; }; } // namespace potential -#endif // _BUCKINGHAM_PAIR_HPP_ \ No newline at end of file +#endif // _BUCKINGHAM_PAIR_HPP_ diff --git a/include/potential/nonCoulomb/forceFieldNonCoulomb.hpp b/include/potential/nonCoulomb/forceFieldNonCoulomb.hpp index 02efe72a4..3c801ab4c 100644 --- a/include/potential/nonCoulomb/forceFieldNonCoulomb.hpp +++ b/include/potential/nonCoulomb/forceFieldNonCoulomb.hpp @@ -25,10 +25,11 @@ #define _FORCE_FIELD_NON_COULOMB_HPP_ #include // for size_t -#include // for map #include +#include #include "nonCoulombPotential.hpp" +#include "strongTypes.hpp" class TestNonCoulombPotentialFF; // forward declaration @@ -47,9 +48,6 @@ namespace potential struct matrix; std::unique_ptr _nonCoulPairsMatPtr; - static constexpr auto _globalVdwType1Index = 4; - static constexpr auto _globalVdwType2Index = 5; - public: ForceFieldNonCoulomb(); ~ForceFieldNonCoulomb() override; @@ -60,7 +58,9 @@ namespace potential ForceFieldNonCoulomb &operator=(ForceFieldNonCoulomb &&) noexcept; void setupNonCoulombicCutoffs(); - void determineInternalGlobalVdwTypes(const std::map &); + void determineInternalGlobalVdwTypes( + const std::unordered_map & + ); void fillDiagOfNonCoulPairsMatrix( std::vector> & ); @@ -68,7 +68,7 @@ namespace potential void sortNonCoulombicsPairs( std::vector> &diagonalElements ); - void setOffDiagonalElement(const size_t, const size_t); + void setOffDiagonalElement(VdwType atomType1, VdwType atomType2); [[nodiscard]] std::vector> findNonCoulPairByInternalTypes( - const size_t, - const size_t + const VdwType intType1, + const VdwType intType2 ) const; void addNonCoulombicPair(const std::shared_ptr &pair); @@ -89,14 +89,10 @@ namespace potential [[nodiscard]] std::shared_ptr getNonCoulPair( - const std::vector &indices + const std::vector &indices, + const std::pair &vdwTypes ) override; - [[nodiscard]] - size_t getGlobalVdwType1(const std::vector &) const; - [[nodiscard]] - size_t getGlobalVdwType2(const std::vector &) const; - [[nodiscard]] std::vector> &getNonCoulombPairsVector( ); diff --git a/include/potential/nonCoulomb/guffNonCoulomb.hpp b/include/potential/nonCoulomb/guffNonCoulomb.hpp index afc9596cd..9756e8c58 100644 --- a/include/potential/nonCoulomb/guffNonCoulomb.hpp +++ b/include/potential/nonCoulomb/guffNonCoulomb.hpp @@ -65,7 +65,8 @@ namespace potential [[nodiscard]] std::shared_ptr getNonCoulPair( - const std::vector &indices + const std::vector &indices, + const std::pair &vdwTypes ) override; [[nodiscard]] std::vector // pair -#include // vector +#include "defaults.hpp" #include "nonCoulombPair.hpp" +struct TestGuffPairUtils; // forward declaration + namespace potential { /** @@ -45,18 +47,28 @@ namespace potential class GuffPair : public NonCoulombPair { private: - std::vector _coefficients; + std::array _coefficients{}; public: - explicit GuffPair(const double, const std::vector &); - explicit GuffPair(const double, const double, const double, const std::vector &); + explicit GuffPair( + const double, + const std::array& + ); + explicit GuffPair( + const double, + const double, + const double, + const std::array& + ); - [[nodiscard]] std::pair calculate(const double distance + [[nodiscard]] + std::pair calculate( + const double distance ) const override; - [[nodiscard]] std::vector getCoefficients() const; + friend struct ::TestGuffPairUtils; }; } // namespace potential -#endif // _GUFF_PAIR_HPP_ \ No newline at end of file +#endif // _GUFF_PAIR_HPP_ diff --git a/include/potential/nonCoulomb/lennardJonesPair.hpp b/include/potential/nonCoulomb/lennardJonesPair.hpp index cbd7306b8..764f841e5 100644 --- a/include/potential/nonCoulomb/lennardJonesPair.hpp +++ b/include/potential/nonCoulomb/lennardJonesPair.hpp @@ -24,13 +24,16 @@ #define _LENNARD_JONES_PAIR_HPP_ -#include // size_t #include // pair #include "nonCoulombPair.hpp" +#include "strongTypes.hpp" + +struct TestLJPairUtils; namespace potential { + /** * @class LennardJonesPair * @@ -41,41 +44,44 @@ namespace potential class LennardJonesPair : public NonCoulombPair { private: - double _c6; - double _c12; + LJParams _params; public: explicit LennardJonesPair( - const size_t vanDerWaalsType1, - const size_t vanDerWaalsType2, - const double cutOff, - const double c6, - const double c12 + const ExtVdwType vanDerWaalsType1, + const ExtVdwType vanDerWaalsType2, + const double cutOff, + const LJParams ¶ms ); + explicit LennardJonesPair(const double cutOff, const LJParams ¶ms); + explicit LennardJonesPair( - const double cutOff, - const double c6, - const double c12 + const double cutOff, + const double energyCutoff, + const double forceCutoff, + const LJParams ¶ms ); + // TODO: we need to explicitly delete it to not implictly create it with + // the wrong types!!! Needs cleanup explicit LennardJonesPair( - const double cutOff, - const double energyCutoff, - const double forceCutoff, - const double c6, - const double c12 - ); + const size_t, + const size_t, + const double, + const LJParams & + ) = delete; - [[nodiscard]] bool operator==(const LennardJonesPair &other) const; + [[nodiscard]] + bool operator==(const LennardJonesPair &other) const; - [[nodiscard]] std::pair calculate(const double distance + [[nodiscard]] std::pair calculate( + const double distance ) const override; - [[nodiscard]] double getC6() const; - [[nodiscard]] double getC12() const; + friend struct ::TestLJPairUtils; }; } // namespace potential -#endif // _LENNARD_JONES_PAIR_HPP_ \ No newline at end of file +#endif // _LENNARD_JONES_PAIR_HPP_ diff --git a/include/potential/nonCoulomb/morsePair.hpp b/include/potential/nonCoulomb/morsePair.hpp index 3e3cbaf66..c78230b5a 100644 --- a/include/potential/nonCoulomb/morsePair.hpp +++ b/include/potential/nonCoulomb/morsePair.hpp @@ -24,10 +24,12 @@ #define _MORSE_PAIR_HPP_ -#include // size_t #include // pair #include "nonCoulombPair.hpp" +#include "strongTypes.hpp" + +struct TestMorsePairUtils; // forward declaration namespace potential { @@ -40,46 +42,43 @@ namespace potential class MorsePair : public NonCoulombPair { private: - double _dissociationEnergy; - double _wellWidth; - double _equilibriumDistance; + MorseParams _params; public: explicit MorsePair( - const size_t vanDerWaalsType1, - const size_t vanDerWaalsType2, - const double cutOff, - const double dissociationEnergy, - const double wellWidth, - const double equilibriumDistance + const ExtVdwType vanDerWaalsType1, + const ExtVdwType vanDerWaalsType2, + const double cutOff, + const MorseParams ¶ms ); + explicit MorsePair(const double cutOff, const MorseParams ¶ms); + explicit MorsePair( - const double cutOff, - const double dissociationEnergy, - const double wellWidth, - const double equilibriumDistance + const double cutOff, + const double energyCutoff, + const double forceCutoff, + const MorseParams ¶ms ); + // TODO: we need to explicitly delete it to not implicitly create it + // with the wrong types!!! Needs cleanup explicit MorsePair( - const double cutOff, - const double energyCutoff, - const double forceCutoff, - const double dissociationEnergy, - const double wellWidth, - const double equilibriumDistance - ); + const size_t, + const size_t, + const double, + const MorseParams & + ) = delete; [[nodiscard]] bool operator==(const MorsePair &other) const; - [[nodiscard]] std::pair calculate(const double distance + [[nodiscard]] std::pair calculate( + const double distance ) const override; - [[nodiscard]] double getDissociationEnergy() const; - [[nodiscard]] double getWellWidth() const; - [[nodiscard]] double getEquilibriumDistance() const; + friend struct ::TestMorsePairUtils; }; } // namespace potential -#endif // _MORSE_PAIR_HPP_ \ No newline at end of file +#endif // _MORSE_PAIR_HPP_ diff --git a/include/potential/nonCoulomb/nonCoulombPair.hpp b/include/potential/nonCoulomb/nonCoulombPair.hpp index c30f52546..a660766f5 100644 --- a/include/potential/nonCoulomb/nonCoulombPair.hpp +++ b/include/potential/nonCoulomb/nonCoulombPair.hpp @@ -24,9 +24,10 @@ #define _NON_COULOMB_PAIR_HPP_ -#include // for size_t #include // for pair +#include "strongTypes.hpp" + namespace potential { /** @@ -42,17 +43,21 @@ namespace potential class NonCoulombPair { protected: - size_t _vanDerWaalsType1 = 0; - size_t _vanDerWaalsType2 = 0; - size_t _internalType1 = 0; - size_t _internalType2 = 0; + ExtVdwType _vanDerWaalsType1{0}; + ExtVdwType _vanDerWaalsType2{0}; + VdwType _internalType1{0}; + VdwType _internalType2{0}; double _radialCutOff; double _energyCutOff = 0.0; double _forceCutOff = 0.0; public: - explicit NonCoulombPair(const size_t, const size_t, const double); + explicit NonCoulombPair( + const ExtVdwType, + const ExtVdwType, + const double + ); explicit NonCoulombPair(const double); explicit NonCoulombPair(const double, const double, const double); @@ -68,8 +73,8 @@ namespace potential * standard setters * ********************/ - void setInternalType1(const size_t internalType1); - void setInternalType2(const size_t internalType2); + void setInternalType1(const VdwType internalType1); + void setInternalType2(const VdwType internalType2); void setRadialCutOff(const double radialCutoff); void setEnergyCutOff(const double energyCutoff); void setForceCutOff(const double forceCutoff); @@ -78,10 +83,10 @@ namespace potential * standard getters * ********************/ - [[nodiscard]] size_t getVanDerWaalsType1() const; - [[nodiscard]] size_t getVanDerWaalsType2() const; - [[nodiscard]] size_t getInternalType1() const; - [[nodiscard]] size_t getInternalType2() const; + [[nodiscard]] ExtVdwType getVanDerWaalsType1() const; + [[nodiscard]] ExtVdwType getVanDerWaalsType2() const; + [[nodiscard]] VdwType getInternalType1() const; + [[nodiscard]] VdwType getInternalType2() const; [[nodiscard]] double getRadialCutOff() const { return _radialCutOff; } [[nodiscard]] double getEnergyCutOff() const; [[nodiscard]] double getForceCutOff() const; diff --git a/include/potential/nonCoulomb/nonCoulombPotential.hpp b/include/potential/nonCoulomb/nonCoulombPotential.hpp index a6672939d..851e50943 100644 --- a/include/potential/nonCoulomb/nonCoulombPotential.hpp +++ b/include/potential/nonCoulomb/nonCoulombPotential.hpp @@ -29,6 +29,7 @@ #include // for vector #include "nonCoulombPair.hpp" +#include "strongTypes.hpp" namespace potential { @@ -62,7 +63,8 @@ namespace potential [[nodiscard]] virtual std::shared_ptr getNonCoulPair( - const std::vector &indices + const std::vector &indices, + const std::pair &vdwTypes ) = 0; [[nodiscard]] MixingRule getMixingRule() const; diff --git a/include/potential/potential.tpp.hpp b/include/potential/potential.tpp.hpp index 28449151a..4e589bb0d 100644 --- a/include/potential/potential.tpp.hpp +++ b/include/potential/potential.tpp.hpp @@ -93,14 +93,8 @@ namespace potential const auto moltype_i = mol1.getMoltype(); const auto moltype_j = mol2.getMoltype(); - const auto combinedIdx = { - moltype_i, - moltype_j, - atomType_i, - atomType_j, - globalVdwType_i, - globalVdwType_j - }; + const auto combinedIdx = + {moltype_i, moltype_j, atomType_i, atomType_j}; const auto charge_i = getPartialCharge(atom1); const auto charge_j = getPartialCharge(atom2); @@ -111,8 +105,10 @@ namespace potential _coulombPotential->calculate(distance, coulombPreFactor); coulombEnergy = e; - const auto nonCoulPair = - _nonCoulombPot->getNonCoulPair(combinedIdx); + const auto nonCoulPair = _nonCoulombPot->getNonCoulPair( + combinedIdx, + {globalVdwType_i, globalVdwType_j} + ); const auto rncCutOff = nonCoulPair->getRadialCutOff(); if (distance < rncCutOff) @@ -259,14 +255,8 @@ namespace potential const auto moltype_i = mol1.getMoltype(); const auto moltype_j = mol2.getMoltype(); - const auto combinedIdx = { - moltype_i, - moltype_j, - atomType_i, - atomType_j, - globalVdwType_i, - globalVdwType_j - }; + const auto combinedIdx = + {moltype_i, moltype_j, atomType_i, atomType_j}; const auto charge_i = getPartialCharge(atom1); const auto charge_j = getPartialCharge(atom2); @@ -277,8 +267,10 @@ namespace potential _coulombPotential->calculate(distance, coulombPreFactor); coulombEnergy = e; - const auto nonCoulPair = - _nonCoulombPot->getNonCoulPair(combinedIdx); + const auto nonCoulPair = _nonCoulombPot->getNonCoulPair( + combinedIdx, + {globalVdwType_i, globalVdwType_j} + ); const auto rncCutOff = nonCoulPair->getRadialCutOff(); if (distance < rncCutOff) diff --git a/include/utilities/mathUtilities.hpp b/include/utilities/mathUtilities.hpp index 617a70210..219dc0327 100644 --- a/include/utilities/mathUtilities.hpp +++ b/include/utilities/mathUtilities.hpp @@ -48,7 +48,8 @@ namespace utilities * @return false */ template - [[nodiscard]] bool compare(const T &a, const T &b, const T &tolerance) + [[nodiscard]] + bool compare(const T &a, const T &b, const T &tolerance) { return std::abs(a - b) < tolerance; } @@ -69,12 +70,14 @@ namespace utilities * @return false */ template - [[nodiscard]] bool compare(const T &a, const T &b) + [[nodiscard]] + bool compare(const T &a, const T &b) { return std::fabs(a - b) < std::numeric_limits::epsilon(); } - [[nodiscard]] bool compare( + [[nodiscard]] + bool compare( const linearAlgebra::Vector3D &a, const linearAlgebra::Vector3D &b ); diff --git a/include/utilities/strongTypes.hpp b/include/utilities/strongTypes.hpp index 0f8c0a3c9..713b6540d 100644 --- a/include/utilities/strongTypes.hpp +++ b/include/utilities/strongTypes.hpp @@ -65,4 +65,70 @@ struct DihedralIdTag }; using DihedralId = StrongSizeT; +struct ExtVdwTypeTag +{ + static std::string toString(const size_t &value) + { + return std::format("ExtVdwType({})", value); + } +}; +using ExtVdwType = StrongSizeT; + +struct VdwTypeTag +{ + static std::string toString(const size_t &value) + { + return std::format("VdwType({})", value); + } +}; +using VdwType = StrongSizeT; + +/** + * @struct LJParams + * + * @brief Struct to hold Lennard-Jones parameters c6 and c12 + * + */ +struct LJParams +{ + double c6{0.0}; + double c12{0.0}; + + [[nodiscard]] + bool operator==(const LJParams &other) const; +}; + +/** + * @struct MorseParams + * + * @brief Struct to hold Morse parameters dissociationEnergy, wellWidth and + * equilibriumDistance + * + */ +struct MorseParams +{ + double dissociationEnergy; + double wellWidth; + double equilibriumDistance; + + [[nodiscard]] + bool operator==(const MorseParams &other) const; +}; + +/** + * @struct BuckinghamParams + * + * @brief Struct to hold Buckingham parameters a, dRho and c6 + * + */ +struct BuckinghamParams +{ + double scaling; + double dRho; + double c6; + + [[nodiscard]] + bool operator==(const BuckinghamParams &other) const; +}; + #endif // _STRONG_TYPES_HPP_ diff --git a/include/waterModels/inter/interWaterParamters.hpp b/include/waterModels/inter/interWaterParamters.hpp index 826ac4655..c543ef178 100644 --- a/include/waterModels/inter/interWaterParamters.hpp +++ b/include/waterModels/inter/interWaterParamters.hpp @@ -77,13 +77,12 @@ namespace waterModel inline static const auto _nonCoulombPairOO = potential::LennardJonesPair( defaults::COULOMB_CUT_OFF_DEFAULT, - constants::SPC_LJ_C6_OO, - constants::SPC_LJ_C12_OO + constants::SPC_LJ_PARAMS_OO ); inline static const auto _nonCoulombPairOH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); inline static const auto _nonCoulombPairHH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); }; struct SPCEInterParam @@ -94,13 +93,12 @@ namespace waterModel inline static const auto _nonCoulombPairOO = potential::LennardJonesPair( defaults::COULOMB_CUT_OFF_DEFAULT, - constants::SPC_E_LJ_C6_OO, - constants::SPC_E_LJ_C12_OO + constants::SPC_E_LJ_PARAMS_OO ); inline static const auto _nonCoulombPairOH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); inline static const auto _nonCoulombPairHH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); }; struct SPCFwInterParam @@ -111,13 +109,12 @@ namespace waterModel inline static const auto _nonCoulombPairOO = potential::LennardJonesPair( defaults::COULOMB_CUT_OFF_DEFAULT, - constants::SPC_FW_LJ_C6_OO, - constants::SPC_FW_LJ_C12_OO + constants::SPC_FW_LJ_PARAMS_OO ); inline static const auto _nonCoulombPairOH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); inline static const auto _nonCoulombPairHH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); }; struct qSPCFwInterParam @@ -128,13 +125,12 @@ namespace waterModel inline static const auto _nonCoulombPairOO = potential::LennardJonesPair( defaults::COULOMB_CUT_OFF_DEFAULT, - constants::QSPC_FW_LJ_C6_OO, - constants::QSPC_FW_LJ_C12_OO + constants::QSPC_FW_LJ_PARAMS_OO ); inline static const auto _nonCoulombPairOH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); inline static const auto _nonCoulombPairHH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); }; struct SPCDCInterParam @@ -145,13 +141,12 @@ namespace waterModel inline static const auto _nonCoulombPairOO = potential::LennardJonesPair( defaults::COULOMB_CUT_OFF_DEFAULT, - constants::SPC_DC_LJ_C6_OO, - constants::SPC_DC_LJ_C12_OO + constants::SPC_DC_LJ_PARAMS_OO ); inline static const auto _nonCoulombPairOH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); inline static const auto _nonCoulombPairHH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); }; struct H2ODCInterParam @@ -162,13 +157,12 @@ namespace waterModel inline static const auto _nonCoulombPairOO = potential::LennardJonesPair( defaults::COULOMB_CUT_OFF_DEFAULT, - constants::H2O_DC_LJ_C6_OO, - constants::H2O_DC_LJ_C12_OO + constants::H2O_DC_LJ_PARAMS_OO ); inline static const auto _nonCoulombPairOH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); inline static const auto _nonCoulombPairHH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); }; struct TIP3PInterParam @@ -179,13 +173,12 @@ namespace waterModel inline static const auto _nonCoulombPairOO = potential::LennardJonesPair( defaults::COULOMB_CUT_OFF_DEFAULT, - constants::TIP3P_LJ_C6_OO, - constants::TIP3P_LJ_C12_OO + constants::TIP3P_LJ_PARAMS_OO ); inline static const auto _nonCoulombPairOH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); inline static const auto _nonCoulombPairHH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); }; struct OPC3InterParam @@ -196,13 +189,12 @@ namespace waterModel inline static const auto _nonCoulombPairOO = potential::LennardJonesPair( defaults::COULOMB_CUT_OFF_DEFAULT, - constants::OPC3_LJ_C6_OO, - constants::OPC3_LJ_C12_OO + constants::OPC3_LJ_PARAMS_OO ); inline static const auto _nonCoulombPairOH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); inline static const auto _nonCoulombPairHH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); }; struct SPCmTRInterParam @@ -213,13 +205,12 @@ namespace waterModel inline static const auto _nonCoulombPairOO = potential::LennardJonesPair( defaults::COULOMB_CUT_OFF_DEFAULT, - constants::SPC_MTR_LJ_C6_OO, - constants::SPC_MTR_LJ_C12_OO + constants::SPC_MTR_LJ_PARAMS_OO ); inline static const auto _nonCoulombPairOH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); inline static const auto _nonCoulombPairHH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); }; struct TIP3PmTRInterParam @@ -230,13 +221,12 @@ namespace waterModel inline static const auto _nonCoulombPairOO = potential::LennardJonesPair( defaults::COULOMB_CUT_OFF_DEFAULT, - constants::TIP3P_MTR_LJ_C6_OO, - constants::TIP3P_MTR_LJ_C12_OO + constants::TIP3P_MTR_LJ_PARAMS_OO ); inline static const auto _nonCoulombPairOH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); inline static const auto _nonCoulombPairHH = - potential::LennardJonesPair(0.01, 0.0, 0.0); + potential::LennardJonesPair(0.01, LJParams{}); }; } // namespace waterModel diff --git a/src/input/guffDatReader.cpp b/src/input/guffDatReader.cpp index 296353e06..c8780b9f3 100644 --- a/src/input/guffDatReader.cpp +++ b/src/input/guffDatReader.cpp @@ -375,12 +375,16 @@ void GuffDatReader::addNonCoulombPair( } case GUFF: { + std::array coeffs; + for (size_t i = 0; i < coefficients.size(); ++i) + coeffs.at(i) = coefficients[i]; + addGuffPair( molType1, molType2, atomType1, atomType2, - coefficients, + coeffs, rncCutOff ); break; @@ -423,31 +427,19 @@ void GuffDatReader::addLennardJonesPair( _engine.getPotential()->getNonCoulombPotential() ); - const auto LJPair = - LennardJonesPair(rncCutOff, coefficients[0], coefficients[2]); + const auto params = LJParams{.c6 = coefficients[0], .c12 = coefficients[2]}; + const auto LJPair = LennardJonesPair(rncCutOff, params); const auto [eCutOff, fCutOff] = LJPair.calculate(rncCutOff); guffNonCoulomb.setGuffNonCoulPair( {molType1, molType2, atomType1, atomType2}, - std::make_shared( - rncCutOff, - eCutOff, - fCutOff, - coefficients[0], - coefficients[2] - ) + std::make_shared(rncCutOff, eCutOff, fCutOff, params) ); guffNonCoulomb.setGuffNonCoulPair( {molType2, molType1, atomType2, atomType1}, - std::make_shared( - rncCutOff, - eCutOff, - fCutOff, - coefficients[0], - coefficients[2] - ) + std::make_shared(rncCutOff, eCutOff, fCutOff, params) ); } @@ -476,36 +468,22 @@ void GuffDatReader::addBuckinghamPair( _engine.getPotential()->getNonCoulombPotential() ); - const auto buckPair = BuckinghamPair( - rncCutOff, - coefficients[0], - coefficients[1], - coefficients[2] - ); + const auto params = BuckinghamParams{ + .scaling = coefficients[0], + .dRho = coefficients[1], + .c6 = coefficients[2] + }; + const auto buckPair = BuckinghamPair(rncCutOff, params); const auto [eCutOff, fCutOff] = buckPair.calculate(rncCutOff); guffNonCoulomb.setGuffNonCoulPair( {molType1, molType2, atomType1, atomType2}, - std::make_shared( - rncCutOff, - eCutOff, - fCutOff, - coefficients[0], - coefficients[1], - coefficients[2] - ) + std::make_shared(rncCutOff, eCutOff, fCutOff, params) ); guffNonCoulomb.setGuffNonCoulPair( {molType2, molType1, atomType2, atomType1}, - std::make_shared( - rncCutOff, - eCutOff, - fCutOff, - coefficients[0], - coefficients[1], - coefficients[2] - ) + std::make_shared(rncCutOff, eCutOff, fCutOff, params) ); } @@ -535,21 +513,18 @@ void GuffDatReader::addMorsePair( _engine.getPotential()->getNonCoulombPotential() ); - // clang-format off - const auto morsePair = MorsePair(rncCutOff, coeffs[0], coeffs[1], coeffs[2]); + const auto params = MorseParams{ + .dissociationEnergy = coeffs[0], + .wellWidth = coeffs[1], + .equilibriumDistance = coeffs[2] + }; + + const auto morsePair = MorsePair(rncCutOff, params); const auto [eCutOff, fCutOff] = morsePair.calculate(rncCutOff); - // clang-format on guffNonCoulomb.setGuffNonCoulPair( {molType1, molType2, atomType1, atomType2}, - std::make_shared( - rncCutOff, - eCutOff, - fCutOff, - coeffs[0], - coeffs[1], - coeffs[2] - ) + std::make_shared(rncCutOff, eCutOff, fCutOff, params) ); guffNonCoulomb.setGuffNonCoulPair( @@ -559,14 +534,7 @@ void GuffDatReader::addMorsePair( atomType2, atomType1, }, - std::make_shared( - rncCutOff, - eCutOff, - fCutOff, - coeffs[0], - coeffs[1], - coeffs[2] - ) + std::make_shared(rncCutOff, eCutOff, fCutOff, params) ); } @@ -581,12 +549,12 @@ void GuffDatReader::addMorsePair( * @param rncCutOff */ void GuffDatReader::addGuffPair( - const size_t molType1, - const size_t molType2, - const size_t atomType1, - const size_t atomType2, - const std::vector &coefficients, - const double rncCutOff + const size_t molType1, + const size_t molType2, + const size_t atomType1, + const size_t atomType2, + const std::array &coefficients, + const double rncCutOff ) { auto &guffNonCoulomb = dynamic_cast( diff --git a/src/input/moldescriptorReader.cpp b/src/input/moldescriptorReader.cpp index fd34c459c..73d89c76f 100644 --- a/src/input/moldescriptorReader.cpp +++ b/src/input/moldescriptorReader.cpp @@ -248,7 +248,8 @@ void MoldescriptorReader::processMolecule( ); } - molecule.addExternalGlobalVDWType(stoul(lineElements[3])); + const auto vdwType = ExtVdwType{stoul(lineElements[3])}; + molecule.addExternalGlobalVDWType(vdwType); } } diff --git a/src/input/parameterFileReader/nonCoulombicsSection.cpp b/src/input/parameterFileReader/nonCoulombicsSection.cpp index 005067143..5206fcb3f 100644 --- a/src/input/parameterFileReader/nonCoulombicsSection.cpp +++ b/src/input/parameterFileReader/nonCoulombicsSection.cpp @@ -34,6 +34,7 @@ #include "morsePair.hpp" // for MorsePair #include "potentialSettings.hpp" // for PotentialSettings #include "stringUtilities.hpp" // for toLowerCopy +#include "strongTypes.hpp" using namespace input::parameterFile; using namespace exc; @@ -172,10 +173,10 @@ void NonCoulombicsSection::processLJ( ); } - const size_t atomType1 = stoul(lineElements[0]); - const size_t atomType2 = stoul(lineElements[1]); - const auto c6 = stod(lineElements[2]); - const auto c12 = stod(lineElements[3]); + const auto atomType1 = ExtVdwType{stoul(lineElements[0])}; + const auto atomType2 = ExtVdwType{stoul(lineElements[1])}; + const auto c6 = stod(lineElements[2]); + const auto c12 = stod(lineElements[3]); auto cutOff = 5 == lineElements.size() ? stod(lineElements[4]) : -1.0; // NOLINTEND(cppcoreguidelines-avoid-magic-numbers) @@ -187,9 +188,10 @@ void NonCoulombicsSection::processLJ( auto &pot = engine.getPotential()->getNonCoulombPotential(); auto &potential = dynamic_cast(pot); + const auto params = LJParams{.c6 = c6, .c12 = c12}; + potential.addNonCoulombicPair( - std::make_shared< - LennardJonesPair>(atomType1, atomType2, cutOff, c6, c12) + std::make_shared(atomType1, atomType2, cutOff, params) ); } @@ -222,7 +224,7 @@ void NonCoulombicsSection::processBuckingham( { throw ParameterFileException( std::format( - "Wrong number of arguments in parameter file in Lennard Jones " + "Wrong number of arguments in parameter file in Buckingham " "nonCoulombics section at line {} - number of " "elements has to be 5 or 6!", _lineNumber @@ -230,11 +232,11 @@ void NonCoulombicsSection::processBuckingham( ); } - const size_t atomType1 = stoul(lineElements[0]); - const size_t atomType2 = stoul(lineElements[1]); - const auto a = stod(lineElements[2]); - const auto dRho = stod(lineElements[3]); - const auto c6 = stod(lineElements[4]); + const auto atomType1 = ExtVdwType{stoul(lineElements[0])}; + const auto atomType2 = ExtVdwType{stoul(lineElements[1])}; + const auto a = stod(lineElements[2]); + const auto dRho = stod(lineElements[3]); + const auto c6 = stod(lineElements[4]); auto cutOff = 6 == lineElements.size() ? stod(lineElements[5]) : -1.0; // NOLINTEND(cppcoreguidelines-avoid-magic-numbers) @@ -243,12 +245,12 @@ void NonCoulombicsSection::processBuckingham( cutOff = cutOff < 0.0 ? coulombCutOff : cutOff; - auto &pot = engine.getPotential()->getNonCoulombPotential(); - auto &potential = dynamic_cast(pot); + auto &pot = engine.getPotential()->getNonCoulombPotential(); + auto &potential = dynamic_cast(pot); + const auto params = BuckinghamParams{.scaling = a, .dRho = dRho, .c6 = c6}; potential.addNonCoulombicPair( - std::make_shared< - BuckinghamPair>(atomType1, atomType2, cutOff, a, dRho, c6) + std::make_shared(atomType1, atomType2, cutOff, params) ); } @@ -281,7 +283,7 @@ void NonCoulombicsSection::processMorse( { throw ParameterFileException( std::format( - "Wrong number of arguments in parameter file in Lennard Jones " + "Wrong number of arguments in parameter file in Morse " "nonCoulombics section at line {} - number of " "elements has to be 5 or 6!", _lineNumber @@ -289,11 +291,11 @@ void NonCoulombicsSection::processMorse( ); } - const size_t atomType1 = stoul(lineElements[0]); - const size_t atomType2 = stoul(lineElements[1]); - const auto dissociationEnergy = stod(lineElements[2]); - const auto wellWidth = stod(lineElements[3]); - const auto equilibriumDistance = stod(lineElements[4]); + const auto atomType1 = ExtVdwType{stoul(lineElements[0])}; + const auto atomType2 = ExtVdwType{stoul(lineElements[1])}; + const auto dissociationEnergy = stod(lineElements[2]); + const auto wellWidth = stod(lineElements[3]); + const auto equilibriumDistance = stod(lineElements[4]); auto cutOff = 6 == lineElements.size() ? stod(lineElements[5]) : -1.0; // NOLINTEND(cppcoreguidelines-avoid-magic-numbers) @@ -305,14 +307,13 @@ void NonCoulombicsSection::processMorse( auto &pot = engine.getPotential()->getNonCoulombPotential(); auto &potential = dynamic_cast(pot); + const auto params = MorseParams{ + .dissociationEnergy = dissociationEnergy, + .wellWidth = wellWidth, + .equilibriumDistance = equilibriumDistance + }; + potential.addNonCoulombicPair( - std::make_shared( - atomType1, - atomType2, - cutOff, - dissociationEnergy, - wellWidth, - equilibriumDistance - ) + std::make_shared(atomType1, atomType2, cutOff, params) ); } diff --git a/src/intraNonBonded/intraNonBondedMap.cpp b/src/intraNonBonded/intraNonBondedMap.cpp index e3cabd07d..40d20cff5 100644 --- a/src/intraNonBonded/intraNonBondedMap.cpp +++ b/src/intraNonBonded/intraNonBondedMap.cpp @@ -166,16 +166,12 @@ std::pair IntraNonBondedMap::calculateSingleInteraction( const auto moltype = _molecule->getMoltype(); - const auto combinedIdx = { - moltype, - moltype, - atomType1, - atomType2, - globalVdwType1, - globalVdwType2 - }; - - const auto nonCoulombicPair = nonCoulPot->getNonCoulPair(combinedIdx); + const auto combinedIdx = {moltype, moltype, atomType1, atomType2}; + + const auto nonCoulombicPair = nonCoulPot->getNonCoulPair( + combinedIdx, + {globalVdwType1, globalVdwType2} + ); if (distance < nonCoulombicPair->getRadialCutOff()) { diff --git a/src/molsys/atom.cpp b/src/molsys/atom.cpp index 3e9d6b319..ca0675fd3 100644 --- a/src/molsys/atom.cpp +++ b/src/molsys/atom.cpp @@ -262,16 +262,22 @@ size_t Atom::getAtomType() const { return _atomType; } /** * @brief return the external global VDW type * - * @return size_t + * @return ExtVdwType */ -size_t Atom::getExternalGlobalVDWType() const { return _externalGlobalVDWType; } +ExtVdwType Atom::getExternalGlobalVDWType() const +{ + return _externalGlobalVDWType; +} /** * @brief return the internal global VDW type * - * @return size_t + * @return VdwType */ -size_t Atom::getInternalGlobalVDWType() const { return _internalGlobalVDWType; } +VdwType Atom::getInternalGlobalVDWType() const +{ + return _internalGlobalVDWType; +} /** * @brief return the mass of the atom @@ -424,7 +430,7 @@ void Atom::setExternalAtomType(const size_t externalAtomType) * * @param externalGlobalVDWType */ -void Atom::setExternalGlobalVDWType(const size_t externalGlobalVDWType) +void Atom::setExternalGlobalVDWType(const ExtVdwType externalGlobalVDWType) { _externalGlobalVDWType = externalGlobalVDWType; } @@ -434,7 +440,7 @@ void Atom::setExternalGlobalVDWType(const size_t externalGlobalVDWType) * * @param internalGlobalVDWType */ -void Atom::setInternalGlobalVDWType(const size_t internalGlobalVDWType) +void Atom::setInternalGlobalVDWType(const VdwType internalGlobalVDWType) { _internalGlobalVDWType = internalGlobalVDWType; } diff --git a/src/molsys/molecule.cpp b/src/molsys/molecule.cpp index 9757f494b..f4a375506 100644 --- a/src/molsys/molecule.cpp +++ b/src/molsys/molecule.cpp @@ -193,14 +193,14 @@ void Molecule::scaleVelocity(const tensor3D &scalingTensor, const Box &box) /** * @brief returns the external global vdw types of the atoms in the molecule * - * @return std::vector + * @return std::vector */ -std::vector Molecule::getExternalGlobalVDWTypes() const +std::vector Molecule::getExternalGlobalVDWTypes() const { - std::vector externalGlobalVDWTypes(getNumberOfAtoms()); + std::vector externalGlobalVDWTypes; - for (size_t i = 0; i < getNumberOfAtoms(); ++i) - externalGlobalVDWTypes[i] = _atoms[i]->getExternalGlobalVDWType(); + for (const auto &atom : _atoms) + externalGlobalVDWTypes.push_back(atom->getExternalGlobalVDWType()); return externalGlobalVDWTypes; } @@ -521,9 +521,9 @@ size_t Molecule::getAtomType(const size_t index) const * @brief returns the internal global vdw type of the atom by index * * @param index - * @return size_t + * @return VdwType */ -size_t Molecule::getInternalGlobalVDWType(const size_t index) const +VdwType Molecule::getInternalGlobalVDWType(const size_t index) const { return _atoms[index]->getInternalGlobalVDWType(); } diff --git a/src/molsys/moleculeType.cpp b/src/molsys/moleculeType.cpp index 1fdb2fc33..0c9dc070d 100644 --- a/src/molsys/moleculeType.cpp +++ b/src/molsys/moleculeType.cpp @@ -91,7 +91,9 @@ void MoleculeType::addPartialCharge(const double partialCharge) * * @param externalGlobalVDWType */ -void MoleculeType::addExternalGlobalVDWType(const size_t externalGlobalVDWType) +void MoleculeType::addExternalGlobalVDWType( + const ExtVdwType externalGlobalVDWType +) { _externalGlobalVDWTypes.push_back(externalGlobalVDWType); } @@ -293,9 +295,9 @@ std::vector &MoleculeType::getExternalAtomTypes() /** * @brief get the external global VDW types of the molecule * - * @return std::vector& + * @return std::vector& */ -std::vector &MoleculeType::getExternalGlobalVDWTypes() +std::vector &MoleculeType::getExternalGlobalVDWTypes() { return _externalGlobalVDWTypes; } diff --git a/src/molsys/simulationBox_standardMethods.cpp b/src/molsys/simulationBox_standardMethods.cpp index 000b19859..edf951d48 100644 --- a/src/molsys/simulationBox_standardMethods.cpp +++ b/src/molsys/simulationBox_standardMethods.cpp @@ -295,9 +295,9 @@ namespace molsys /** * @brief get the external global VDW types * - * @return std::vector& + * @return std::vector& */ - std::vector &SimulationBox::getExternalGlobalVdwTypes() + std::vector &SimulationBox::getExternalGlobalVdwTypes() { return _externalGlobalVdwTypes; } @@ -305,9 +305,9 @@ namespace molsys /** * @brief get the external to internal global VDW types map * - * @return std::unordered_map& + * @return std::unordered_map& */ - std::map &SimulationBox:: + std::unordered_map &SimulationBox:: getExternalToInternalGlobalVDWTypes() { return _externalToInternalGlobalVDWTypes; diff --git a/src/potential/nonCoulomb/buckinghamPair.cpp b/src/potential/nonCoulomb/buckinghamPair.cpp index a0ae77b76..b51c04a50 100644 --- a/src/potential/nonCoulomb/buckinghamPair.cpp +++ b/src/potential/nonCoulomb/buckinghamPair.cpp @@ -24,10 +24,7 @@ #include // for exp -#include "mathUtilities.hpp" // for compare - using namespace potential; -using namespace utilities; /** * @brief Construct a new Buckingham Pair:: Buckingham Pair object @@ -40,17 +37,13 @@ using namespace utilities; * @param c6 */ BuckinghamPair::BuckinghamPair( - const size_t vanDerWaalsType1, - const size_t vanDerWaalsType2, - const double cutOff, - const double a, - const double dRho, - const double c6 + const ExtVdwType vanDerWaalsType1, + const ExtVdwType vanDerWaalsType2, + const double cutOff, + const BuckinghamParams& params ) : NonCoulombPair(vanDerWaalsType1, vanDerWaalsType2, cutOff), - _a(a), - _dRho(dRho), - _c6(c6) + _params(params) { } @@ -63,12 +56,10 @@ BuckinghamPair::BuckinghamPair( * @param c6 */ BuckinghamPair::BuckinghamPair( - const double cutOff, - const double a, - const double dRho, - const double c6 + const double cutOff, + const BuckinghamParams& params ) - : NonCoulombPair(cutOff), _a(a), _dRho(dRho), _c6(c6) + : NonCoulombPair(cutOff), _params(params) { } @@ -83,17 +74,12 @@ BuckinghamPair::BuckinghamPair( * @param c6 */ BuckinghamPair::BuckinghamPair( - const double cutOff, - const double energyCutoff, - const double forceCutoff, - const double a, - const double dRho, - const double c6 + const double cutOff, + const double energyCutoff, + const double forceCutoff, + const BuckinghamParams& params ) - : NonCoulombPair(cutOff, energyCutoff, forceCutoff), - _a(a), - _dRho(dRho), - _c6(c6) + : NonCoulombPair(cutOff, energyCutoff, forceCutoff), _params(params) { } @@ -104,14 +90,9 @@ BuckinghamPair::BuckinghamPair( * @return true * @return false */ -bool BuckinghamPair::operator==(const BuckinghamPair &other) const +bool BuckinghamPair::operator==(const BuckinghamPair& other) const { - auto isEqual = NonCoulombPair::operator==(other); - isEqual = isEqual && compare(_a, other._a); - isEqual = isEqual && compare(_dRho, other._dRho); - isEqual = isEqual && compare(_c6, other._c6); - - return isEqual; + return NonCoulombPair::operator==(other) && _params == other._params; } /** @@ -126,42 +107,15 @@ std::pair BuckinghamPair::calculate(const double distance) const { const auto distanceThird = distance * distance * distance; const auto distanceSixth = distanceThird * distanceThird; - const auto expTerm = _a * ::exp(_dRho * distance); + const auto expTerm = _params.scaling * ::exp(_params.dRho * distance); - auto energy = expTerm + _c6 / distanceSixth - _energyCutOff; + auto energy = expTerm + _params.c6 / distanceSixth - _energyCutOff; energy -= _forceCutOff * (_radialCutOff - distance); - auto force = -_dRho * expTerm; + auto force = -_params.dRho * expTerm; // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers) - force += 6.0 * _c6 / (distanceSixth * distance) - _forceCutOff; + force += 6.0 * _params.c6 / (distanceSixth * distance) - _forceCutOff; return {energy, force}; } - -/*************************** - * * - * standard getter methods * - * * - ***************************/ - -/** - * @brief get the A parameter - * - * @return double - */ -double BuckinghamPair::getA() const { return _a; } - -/** - * @brief get the dRho parameter - * - * @return double - */ -double BuckinghamPair::getDRho() const { return _dRho; } - -/** - * @brief get the C6 parameter - * - * @return double - */ -double BuckinghamPair::getC6() const { return _c6; } diff --git a/src/potential/nonCoulomb/forceFieldNonCoulomb.cpp b/src/potential/nonCoulomb/forceFieldNonCoulomb.cpp index 678c4143b..74a6e15cd 100644 --- a/src/potential/nonCoulomb/forceFieldNonCoulomb.cpp +++ b/src/potential/nonCoulomb/forceFieldNonCoulomb.cpp @@ -24,7 +24,6 @@ #include // for copy, max #include // for std::format -#include // for map #include // for __find_if_fn, find_if #include // for string_view @@ -33,6 +32,7 @@ #include "lennardJonesPair.hpp" // IWYU pragma: keep -- for template instantiation #include "matrix.hpp" #include "nonCoulombPair.hpp" // for NonCoulombPair +#include "strongTypes.hpp" using namespace potential; using namespace exc; @@ -115,7 +115,7 @@ void ForceFieldNonCoulomb::setupNonCoulombicCutoffs() * */ void ForceFieldNonCoulomb::determineInternalGlobalVdwTypes( - const std::map &extToIntGlobalVDWTypes + const std::unordered_map &extToIntGlobalVDWTypes ) { auto setIntVDWType = [&extToIntGlobalVDWTypes](auto &nonCoulombPair) @@ -172,8 +172,8 @@ void ForceFieldNonCoulomb::sortNonCoulombicsPairs( "Non-coulombic pairs with global van der Waals types {} and {} " "in " "the parameter file are defined twice", - (*iter)->getVanDerWaalsType1(), - (*iter)->getVanDerWaalsType2() + (*iter)->getVanDerWaalsType1().toString(), + (*iter)->getVanDerWaalsType2().toString() ) ); } @@ -220,8 +220,8 @@ void ForceFieldNonCoulomb::fillDiagOfNonCoulPairsMatrix( * index combinations and it has different parameters */ void ForceFieldNonCoulomb::setOffDiagonalElement( - const size_t atomType1, - const size_t atomType2 + const VdwType atomType1, + const VdwType atomType2 ) { auto nonCoulPair1 = findNonCoulPairByInternalTypes(atomType1, atomType2); @@ -250,26 +250,32 @@ void ForceFieldNonCoulomb::setOffDiagonalElement( "{} " "and {}, {} in the parameter file have " "different parameters", - vdwType1, - vdwType2, - vdwType2, - vdwType1 + vdwType1.toString(), + vdwType2.toString(), + vdwType2.toString(), + vdwType1.toString() ) ); } - _nonCoulPairsMatPtr->matrix(atomType1, atomType2) = *nonCoulPair1; - _nonCoulPairsMatPtr->matrix(atomType2, atomType1) = *nonCoulPair1; + _nonCoulPairsMatPtr->matrix(atomType1.get(), atomType2.get()) = + *nonCoulPair1; + _nonCoulPairsMatPtr->matrix(atomType2.get(), atomType1.get()) = + *nonCoulPair1; } else if (nonCoulPair1 != std::nullopt) { - _nonCoulPairsMatPtr->matrix(atomType1, atomType2) = *nonCoulPair1; - _nonCoulPairsMatPtr->matrix(atomType2, atomType1) = *nonCoulPair1; + _nonCoulPairsMatPtr->matrix(atomType1.get(), atomType2.get()) = + *nonCoulPair1; + _nonCoulPairsMatPtr->matrix(atomType2.get(), atomType1.get()) = + *nonCoulPair1; } else { - _nonCoulPairsMatPtr->matrix(atomType1, atomType2) = *nonCoulPair2; - _nonCoulPairsMatPtr->matrix(atomType2, atomType1) = *nonCoulPair2; + _nonCoulPairsMatPtr->matrix(atomType1.get(), atomType2.get()) = + *nonCoulPair2; + _nonCoulPairsMatPtr->matrix(atomType2.get(), atomType1.get()) = + *nonCoulPair2; } } @@ -282,7 +288,8 @@ void ForceFieldNonCoulomb::fillOffDiagOfNonCoulPairsMatrix() const auto &[rows, cols] = _nonCoulPairsMatPtr->matrix.shape(); for (size_t i = 0; i < rows; ++i) - for (size_t j = i + 1; j < cols; ++j) setOffDiagonalElement(i, j); + for (size_t j = i + 1; j < cols; ++j) + setOffDiagonalElement(VdwType{i}, VdwType{j}); } /** @@ -323,8 +330,8 @@ std::vector> ForceFieldNonCoulomb:: */ std::optional> ForceFieldNonCoulomb:: findNonCoulPairByInternalTypes( - const size_t intType1, - const size_t intType2 + const VdwType intType1, + const VdwType intType2 ) const { auto findByIntAtomTypes = [intType1, intType2](const auto &nonCoulPair) @@ -358,8 +365,8 @@ std::optional> ForceFieldNonCoulomb:: "Non coulombic pair with global van der waals types {} and " "{} " "is defined twice in the parameter file.", - vdwType1, - vdwType2 + vdwType1.toString(), + vdwType2.toString() ) ); } @@ -398,41 +405,16 @@ void ForceFieldNonCoulomb::addNonCoulombicPair( * @return std::shared_ptr */ std::shared_ptr ForceFieldNonCoulomb::getNonCoulPair( - const std::vector &indices + const std::vector & /*indices*/, + const std::pair &vdwTypes ) { - const auto idx1 = getGlobalVdwType1(indices); - const auto idx2 = getGlobalVdwType2(indices); + const auto idx1 = vdwTypes.first.get(); + const auto idx2 = vdwTypes.second.get(); return _nonCoulPairsMatPtr->matrix(idx1, idx2); } -/** - * @brief get the global van der Waals type 1 - * - * @param indices - * @return size_t - */ -size_t ForceFieldNonCoulomb::getGlobalVdwType1( - const std::vector &indices -) const -{ - return indices[_globalVdwType1Index]; -} - -/** - * @brief get the global van der Waals type 2 - * - * @param indices - * @return size_t - */ -size_t ForceFieldNonCoulomb::getGlobalVdwType2( - const std::vector &indices -) const -{ - return indices[_globalVdwType2Index]; -} - /** * @brief Get the Non Coulomb Pairs Vector object * diff --git a/src/potential/nonCoulomb/guffNonCoulomb.cpp b/src/potential/nonCoulomb/guffNonCoulomb.cpp index 3b4e26ecd..07ad656d6 100644 --- a/src/potential/nonCoulomb/guffNonCoulomb.cpp +++ b/src/potential/nonCoulomb/guffNonCoulomb.cpp @@ -122,7 +122,8 @@ void GuffNonCoulomb::setGuffNonCoulPair( * @return std::shared_ptr */ std::shared_ptr GuffNonCoulomb::getNonCoulPair( - const std::vector &indices + const std::vector &indices, + const std::pair & /*vdwTypes*/ ) { const auto m1 = getMolType1(indices) - 1; diff --git a/src/potential/nonCoulomb/guffPair.cpp b/src/potential/nonCoulomb/guffPair.cpp index 7dd71e9fd..fe5aa1609 100644 --- a/src/potential/nonCoulomb/guffPair.cpp +++ b/src/potential/nonCoulomb/guffPair.cpp @@ -32,7 +32,10 @@ using namespace potential; * @param cutOff * @param coefficients */ -GuffPair::GuffPair(const double cutOff, const std::vector &coefficients) +GuffPair::GuffPair( + const double cutOff, + const std::array &coefficients +) : NonCoulombPair(cutOff), _coefficients(coefficients) { } @@ -46,10 +49,10 @@ GuffPair::GuffPair(const double cutOff, const std::vector &coefficients) * @param coefficients */ GuffPair::GuffPair( - const double cutOff, - const double energyCutoff, - const double forceCutoff, - const std::vector &coefficients + const double cutOff, + const double energyCutoff, + const double forceCutoff, + const std::array &coefficients ) : NonCoulombPair(cutOff, energyCutoff, forceCutoff), _coefficients(coefficients) @@ -72,40 +75,40 @@ std::pair GuffPair::calculate(const double distance) const double energy = 0.0; double force = 0.0; - if (const double c1 = _coefficients[0]; c1 != 0.0) + if (const double c1 = _coefficients.at(0); c1 != 0.0) { - const double n2 = _coefficients[1]; - const double distance_n2 = ::pow(distance, n2); - energy += c1 / distance_n2; - force += n2 * c1 / (distance_n2 * distance); + const double n2 = _coefficients.at(1); + const double distance_n2 = ::pow(distance, n2); + energy += c1 / distance_n2; + force += n2 * c1 / (distance_n2 * distance); } - if (const double c3 = _coefficients[2]; c3 != 0.0) + if (const double c3 = _coefficients.at(2); c3 != 0.0) { - const double n4 = _coefficients[3]; - const double distance_n4 = ::pow(distance, n4); - energy += c3 / distance_n4; - force += n4 * c3 / (distance_n4 * distance); + const double n4 = _coefficients.at(3); + const double distance_n4 = ::pow(distance, n4); + energy += c3 / distance_n4; + force += n4 * c3 / (distance_n4 * distance); } - if (const double c5 = _coefficients[4]; c5 != 0.0) + if (const double c5 = _coefficients.at(4); c5 != 0.0) { - const double n6 = _coefficients[5]; - const double distance_n6 = ::pow(distance, n6); - energy += c5 / distance_n6; - force += n6 * c5 / (distance_n6 * distance); + const double n6 = _coefficients.at(5); + const double distance_n6 = ::pow(distance, n6); + energy += c5 / distance_n6; + force += n6 * c5 / (distance_n6 * distance); } - if (const double c7 = _coefficients[6]; c7 != 0.0) + if (const double c7 = _coefficients.at(6); c7 != 0.0) { - const double n8 = _coefficients[7]; - const double distance_n8 = ::pow(distance, n8); - energy += c7 / distance_n8; - force += n8 * c7 / (distance_n8 * distance); + const double n8 = _coefficients.at(7); + const double distance_n8 = ::pow(distance, n8); + energy += c7 / distance_n8; + force += n8 * c7 / (distance_n8 * distance); } - if (const double c9 = _coefficients[8]; c9 != 0.0) + if (const double c9 = _coefficients.at(8); c9 != 0.0) { - const double cexp10 = _coefficients[9]; - const double rExp11 = _coefficients[10]; + const double cexp10 = _coefficients.at(9); + const double rExp11 = _coefficients.at(10); const double helper = ::exp(cexp10 * (distance - rExp11)); @@ -113,10 +116,10 @@ std::pair GuffPair::calculate(const double distance) const force += c9 * cexp10 * helper / ((1 + helper) * (1 + helper)); } - if (const double c12 = _coefficients[11]; c12 != 0.0) + if (const double c12 = _coefficients.at(11); c12 != 0.0) { - const double cexp13 = _coefficients[12]; - const double rExp14 = _coefficients[13]; + const double cexp13 = _coefficients.at(12); + const double rExp14 = _coefficients.at(13); const double helper = ::exp(cexp13 * (distance - rExp14)); @@ -124,11 +127,11 @@ std::pair GuffPair::calculate(const double distance) const force += c12 * cexp13 * helper / ((1 + helper) * (1 + helper)); } - if (const double c15 = _coefficients[14]; c15 != 0.0) + if (const double c15 = _coefficients.at(14); c15 != 0.0) { - const double cexp16 = _coefficients[15]; - const double rExp17 = _coefficients[16]; - const double n18 = _coefficients[17]; + const double cexp16 = _coefficients.at(15); + const double rExp17 = _coefficients.at(16); + const double n18 = _coefficients.at(17); const double distance_n18 = ::pow(distance - rExp17, n18); const double helper = c15 * ::exp(cexp16 * distance_n18); @@ -137,11 +140,11 @@ std::pair GuffPair::calculate(const double distance) const force += -cexp16 * n18 * distance_n18 / (distance - rExp17) * helper; } - if (const double c19 = _coefficients[18]; c19 != 0.0) + if (const double c19 = _coefficients.at(18); c19 != 0.0) { - const double cexp20 = _coefficients[19]; - const double rExp21 = _coefficients[20]; - const double n22 = _coefficients[21]; + const double cexp20 = _coefficients.at(19); + const double rExp21 = _coefficients.at(20); + const double n22 = _coefficients.at(21); const double distance_n22 = ::pow(distance - rExp21, n22); const double helper = c19 * ::exp(cexp20 * distance_n22); @@ -155,10 +158,3 @@ std::pair GuffPair::calculate(const double distance) const return {energy, force}; } - -/** - * @brief get the coefficients of the GuffPair - * - * @return std::vector - */ -std::vector GuffPair::getCoefficients() const { return _coefficients; } diff --git a/src/potential/nonCoulomb/lennardJonesPair.cpp b/src/potential/nonCoulomb/lennardJonesPair.cpp index 8ed9dc8cf..e4e4454ef 100644 --- a/src/potential/nonCoulomb/lennardJonesPair.cpp +++ b/src/potential/nonCoulomb/lennardJonesPair.cpp @@ -22,10 +22,9 @@ #include "lennardJonesPair.hpp" -#include "mathUtilities.hpp" // for compare +#include "strongTypes.hpp" using namespace potential; -using namespace utilities; /** * @brief Construct a new Lennard Jones Pair:: Lennard Jones Pair object @@ -33,19 +32,16 @@ using namespace utilities; * @param vanDerWaalsType1 * @param vanDerWaalsType2 * @param cutOff - * @param c6 - * @param c12 + * @param params */ LennardJonesPair::LennardJonesPair( - const size_t vanDerWaalsType1, - const size_t vanDerWaalsType2, - const double cutOff, - const double c6, - const double c12 + const ExtVdwType vanDerWaalsType1, + const ExtVdwType vanDerWaalsType2, + const double cutOff, + const LJParams ¶ms ) : NonCoulombPair(vanDerWaalsType1, vanDerWaalsType2, cutOff), - _c6(c6), - _c12(c12) + _params(params) { } @@ -56,12 +52,8 @@ LennardJonesPair::LennardJonesPair( * @param c6 * @param c12 */ -LennardJonesPair::LennardJonesPair( - const double cutOff, - const double c6, - const double c12 -) - : NonCoulombPair(cutOff), _c6(c6), _c12(c12) +LennardJonesPair::LennardJonesPair(const double cutOff, const LJParams ¶ms) + : NonCoulombPair(cutOff), _params(params) { } @@ -75,13 +67,12 @@ LennardJonesPair::LennardJonesPair( * @param c12 */ LennardJonesPair::LennardJonesPair( - const double cutOff, - const double energyCutoff, - const double forceCutoff, - const double c6, - const double c12 + const double cutOff, + const double energyCutoff, + const double forceCutoff, + const LJParams ¶ms ) - : NonCoulombPair(cutOff, energyCutoff, forceCutoff), _c6(c6), _c12(c12) + : NonCoulombPair(cutOff, energyCutoff, forceCutoff), _params(params) { } @@ -96,8 +87,7 @@ bool LennardJonesPair::operator==(const LennardJonesPair &other) const { auto isEqual = true; isEqual = isEqual && NonCoulombPair::operator==(other); - isEqual = isEqual && compare(_c6, other._c6); - isEqual = isEqual && compare(_c12, other._c12); + isEqual = isEqual && _params == other._params; return isEqual; } @@ -116,30 +106,16 @@ std::pair LennardJonesPair::calculate( const auto distanceSixth = distanceThird * distanceThird; const auto distanceTwelfth = distanceSixth * distanceSixth; - auto energy = _c12 / distanceTwelfth; - energy += _c6 / distanceSixth; + auto energy = _params.c12 / distanceTwelfth; + energy += _params.c6 / distanceSixth; energy -= _energyCutOff; energy -= _forceCutOff * (_radialCutOff - distance); // NOLINTBEGIN(cppcoreguidelines-avoid-magic-numbers) - auto force = 12.0 * _c12 / (distanceTwelfth * distance); - force += 6.0 * _c6 / (distanceSixth * distance); + auto force = 12.0 * _params.c12 / (distanceTwelfth * distance); + force += 6.0 * _params.c6 / (distanceSixth * distance); force -= _forceCutOff; // NOLINTEND(cppcoreguidelines-avoid-magic-numbers) return {energy, force}; } - -/** - * @brief get the c6 and c12 coefficients - * - * @return double - */ -double LennardJonesPair::getC6() const { return _c6; } - -/** - * @brief get the c6 and c12 coefficients - * - * @return double - */ -double LennardJonesPair::getC12() const { return _c12; } diff --git a/src/potential/nonCoulomb/morsePair.cpp b/src/potential/nonCoulomb/morsePair.cpp index da8b82b43..bcaa280ba 100644 --- a/src/potential/nonCoulomb/morsePair.cpp +++ b/src/potential/nonCoulomb/morsePair.cpp @@ -24,10 +24,7 @@ #include // for exp -#include "mathUtilities.hpp" // for compare - using namespace potential; -using namespace utilities; /** * @brief Construct a new Morse Pair:: Morse Pair object @@ -35,22 +32,18 @@ using namespace utilities; * @param vanDerWaalsType1 * @param vanDerWaalsType2 * @param cutOff - * @param dissociationEnergy - * @param wellWidth - * @param equilibriumDistance + * @param params */ MorsePair::MorsePair( - const size_t vanDerWaalsType1, - const size_t vanDerWaalsType2, - const double cutOff, - const double dissociationEnergy, - const double wellWidth, - const double equilibriumDistance + const ExtVdwType vanDerWaalsType1, + const ExtVdwType vanDerWaalsType2, + const double cutOff, + const MorseParams ¶ms ) : NonCoulombPair(vanDerWaalsType1, vanDerWaalsType2, cutOff), - _dissociationEnergy(dissociationEnergy), - _wellWidth(wellWidth), - _equilibriumDistance(equilibriumDistance){}; + _params(params) +{ +} /** * @brief Construct a new Morse Pair:: Morse Pair object @@ -60,16 +53,10 @@ MorsePair::MorsePair( * @param wellWidth * @param equilibriumDistance */ -MorsePair::MorsePair( - const double cutOff, - const double dissociationEnergy, - const double wellWidth, - const double equilibriumDistance -) - : NonCoulombPair(cutOff), - _dissociationEnergy(dissociationEnergy), - _wellWidth(wellWidth), - _equilibriumDistance(equilibriumDistance){}; +MorsePair::MorsePair(const double cutOff, const MorseParams ¶ms) + : NonCoulombPair(cutOff), _params(params) +{ +} /** * @brief Construct a new Morse Pair:: Morse Pair object @@ -77,22 +64,17 @@ MorsePair::MorsePair( * @param cutOff * @param energyCutoff * @param forceCutoff - * @param dissociationEnergy - * @param wellWidth - * @param equilibriumDistance + * @param params */ MorsePair::MorsePair( - const double cutOff, - const double energyCutoff, - const double forceCutoff, - const double dissociationEnergy, - const double wellWidth, - const double equilibriumDistance + const double cutOff, + const double energyCutoff, + const double forceCutoff, + const MorseParams ¶ms ) - : NonCoulombPair(cutOff, energyCutoff, forceCutoff), - _dissociationEnergy(dissociationEnergy), - _wellWidth(wellWidth), - _equilibriumDistance(equilibriumDistance){}; + : NonCoulombPair(cutOff, energyCutoff, forceCutoff), _params(params) +{ +} /** * @brief operator overload for the comparison of two MorsePair objects @@ -106,9 +88,7 @@ bool MorsePair::operator==(const MorsePair &other) const auto isEq = true; isEq = isEq && NonCoulombPair::operator==(other); - isEq = isEq && compare(_dissociationEnergy, other._dissociationEnergy); - isEq = isEq && compare(_wellWidth, other._wellWidth); - isEq = isEq && compare(_equilibriumDistance, other._equilibriumDistance); + isEq = isEq && _params == other._params; return isEq; } @@ -121,41 +101,18 @@ bool MorsePair::operator==(const MorsePair &other) const */ std::pair MorsePair::calculate(const double distance) const { - const auto deltaEquilibrium = distance - _equilibriumDistance; - const auto expTerm = std::exp(-_wellWidth * deltaEquilibrium); - const auto oneMinusExpTerm = 1.0 - expTerm; + const auto deltaEquilibrium = distance - _params.equilibriumDistance; + const auto expTerm = std::exp(-_params.wellWidth * deltaEquilibrium); + const auto oneMinusExpTerm = 1.0 - expTerm; - auto energy = _dissociationEnergy * oneMinusExpTerm * oneMinusExpTerm; - energy -= _energyCutOff; - energy -= _forceCutOff * (_radialCutOff - distance); + auto energy = + _params.dissociationEnergy * oneMinusExpTerm * oneMinusExpTerm; + energy -= _energyCutOff; + energy -= _forceCutOff * (_radialCutOff - distance); - auto force = -2.0 * _dissociationEnergy * _wellWidth; + auto force = -2.0 * _params.dissociationEnergy * _params.wellWidth; force *= expTerm * oneMinusExpTerm; force -= _forceCutOff; return {energy, force}; } - -/** - * @brief get the dissociation energy - * - * @return double - */ -double MorsePair::getDissociationEnergy() const { return _dissociationEnergy; } - -/** - * @brief get the well width - * - * @return double - */ -double MorsePair::getWellWidth() const { return _wellWidth; } - -/** - * @brief get the equilibrium distance - * - * @return double - */ -double MorsePair::getEquilibriumDistance() const -{ - return _equilibriumDistance; -} \ No newline at end of file diff --git a/src/potential/nonCoulomb/nonCoulombPair.cpp b/src/potential/nonCoulomb/nonCoulombPair.cpp index f9195884a..291c901a5 100644 --- a/src/potential/nonCoulomb/nonCoulombPair.cpp +++ b/src/potential/nonCoulomb/nonCoulombPair.cpp @@ -35,20 +35,22 @@ using namespace utilities; * @param cutOff */ NonCoulombPair::NonCoulombPair( - const size_t vanDerWaalsType1, - const size_t vanDerWaalsType2, - const double cutOff + const ExtVdwType vanDerWaalsType1, + const ExtVdwType vanDerWaalsType2, + const double cutOff ) : _vanDerWaalsType1(vanDerWaalsType1), _vanDerWaalsType2(vanDerWaalsType2), - _radialCutOff(cutOff){}; + _radialCutOff(cutOff) +{ +} /** * @brief Construct a new Non Coulomb Pair:: Non Coulomb Pair object * * @param cutOff */ -NonCoulombPair::NonCoulombPair(const double cutOff) : _radialCutOff(cutOff){}; +NonCoulombPair::NonCoulombPair(const double cutOff) : _radialCutOff(cutOff) {} /** * @brief Construct a new Non Coulomb Pair:: Non Coulomb Pair object @@ -64,7 +66,9 @@ NonCoulombPair::NonCoulombPair( ) : _radialCutOff(cutoff), _energyCutOff(energyCutoff), - _forceCutOff(forceCutoff){}; + _forceCutOff(forceCutoff) +{ +} /** * @brief operator overload for the comparison of two NonCoulombPair objects @@ -103,7 +107,7 @@ bool NonCoulombPair::operator==(const NonCoulombPair &other) const * * @param internalType1 */ -void NonCoulombPair::setInternalType1(const size_t internalType1) +void NonCoulombPair::setInternalType1(const VdwType internalType1) { _internalType1 = internalType1; } @@ -113,7 +117,7 @@ void NonCoulombPair::setInternalType1(const size_t internalType1) * * @param internalType2 */ -void NonCoulombPair::setInternalType2(const size_t internalType2) +void NonCoulombPair::setInternalType2(const VdwType internalType2) { _internalType2 = internalType2; } @@ -157,30 +161,36 @@ void NonCoulombPair::setForceCutOff(const double forceCutoff) /** * @brief get van der Waals type 1 * - * @return size_t + * @return ExtVdwType */ -size_t NonCoulombPair::getVanDerWaalsType1() const { return _vanDerWaalsType1; } +ExtVdwType NonCoulombPair::getVanDerWaalsType1() const +{ + return _vanDerWaalsType1; +} /** * @brief get van der Waals type 2 * - * @return size_t + * @return ExtVdwType */ -size_t NonCoulombPair::getVanDerWaalsType2() const { return _vanDerWaalsType2; } +ExtVdwType NonCoulombPair::getVanDerWaalsType2() const +{ + return _vanDerWaalsType2; +} /** * @brief get internal type 1 * - * @return size_t + * @return VdwType */ -size_t NonCoulombPair::getInternalType1() const { return _internalType1; } +VdwType NonCoulombPair::getInternalType1() const { return _internalType1; } /** * @brief get internal type 2 * - * @return size_t + * @return VdwType */ -size_t NonCoulombPair::getInternalType2() const { return _internalType2; } +VdwType NonCoulombPair::getInternalType2() const { return _internalType2; } /** * @brief get energy cut off diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt index e214d195c..c8381964a 100644 --- a/src/utilities/CMakeLists.txt +++ b/src/utilities/CMakeLists.txt @@ -2,6 +2,7 @@ add_library(utilities executablePath.cpp stringUtilities.cpp mathUtilities.cpp + strongTypes.cpp ) target_include_directories(utilities diff --git a/src/utilities/strongTypes.cpp b/src/utilities/strongTypes.cpp new file mode 100644 index 000000000..703db3768 --- /dev/null +++ b/src/utilities/strongTypes.cpp @@ -0,0 +1,66 @@ +/***************************************************************************** + + + PQ + Copyright (C) 2023-now Jakob Gamper + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + + +******************************************************************************/ + +#include "strongTypes.hpp" + +#include "mathUtilities.hpp" + +/** + * @brief operator overload for the comparison of two LJParams objects + * + * @param other + * @return true + * @return false + */ +bool LJParams::operator==(const LJParams &other) const +{ + return utilities::compare(c6, other.c6) && + utilities::compare(c12, other.c12); +} + +/** + * @brief compare two MorseParams objects for equality + * + * @param other + * @return true + * @return false + */ +bool MorseParams::operator==(const MorseParams &other) const +{ + return utilities::compare(dissociationEnergy, other.dissociationEnergy) && + utilities::compare(wellWidth, other.wellWidth) && + utilities::compare(equilibriumDistance, other.equilibriumDistance); +} + +/** + * @brief compare two BuckinghamParams objects for equality + * + * @param other + * @return true + * @return false + */ +bool BuckinghamParams::operator==(const BuckinghamParams &other) const +{ + return utilities::compare(scaling, other.scaling) && + utilities::compare(dRho, other.dRho) && + utilities::compare(c6, other.c6); +} diff --git a/tests/include/testUtils/testNonCoulombPairUtils.hpp b/tests/include/testUtils/testNonCoulombPairUtils.hpp new file mode 100644 index 000000000..ef520c5bc --- /dev/null +++ b/tests/include/testUtils/testNonCoulombPairUtils.hpp @@ -0,0 +1,71 @@ +/***************************************************************************** + + + PQ + Copyright (C) 2023-now Jakob Gamper + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + + +******************************************************************************/ + +#ifndef _TEST_MORSE_PAIR_UTILS_ +#define _TEST_MORSE_PAIR_UTILS_ + +#include "buckinghamPair.hpp" +#include "guffPair.hpp" +#include "lennardJonesPair.hpp" +#include "morsePair.hpp" + +/** + * @brief struct TestMorsePairUtils + * + */ +struct TestMorsePairUtils +{ + static const MorseParams& params(const potential::MorsePair* morsePair); +}; + +/** + * @brief struct TestLJPairUtils + * + */ +struct TestLJPairUtils +{ + static const LJParams& params(const potential::LennardJonesPair* ljPair); +}; + +/** + * @brief struct TestBuckinghamPairUtils + * + */ +struct TestBuckinghamPairUtils +{ + static const BuckinghamParams& params( + const potential::BuckinghamPair* buckPair + ); +}; + +/** + * @brief struct TestGuffPairUtils + * + */ +struct TestGuffPairUtils +{ + static const std::array& coeffs( + const potential::GuffPair* guffPair + ); +}; + +#endif diff --git a/tests/src/forceField/testAngleForceField.cpp b/tests/src/forceField/testAngleForceField.cpp index 75754ca51..93092d7a9 100644 --- a/tests/src/forceField/testAngleForceField.cpp +++ b/tests/src/forceField/testAngleForceField.cpp @@ -58,11 +58,10 @@ TEST_F(TestAngleForceField, calculateEnergyAndForces) auto coulombPotential = potential::CoulombShiftedPotential(10.0); auto nonCoulombPair = potential::LennardJonesPair( - static_cast(1), - static_cast(1), + ExtVdwType(1), + ExtVdwType(1), 5.0, - 2.0, - 4.0 + LJParams{.c6 = 2.0, .c12 = 4.0} ); setNonCoulombPairsMatrix( linearAlgebra::Matrix>(2, 2) @@ -86,9 +85,9 @@ TEST_F(TestAngleForceField, calculateEnergyAndForces) atom2->setForce({0.0, 0.0, 0.0}); atom3->setForce({0.0, 0.0, 0.0}); - atom1->setInternalGlobalVDWType(0); - atom2->setInternalGlobalVDWType(1); - atom3->setInternalGlobalVDWType(1); + atom1->setInternalGlobalVDWType(VdwType{0}); + atom2->setInternalGlobalVDWType(VdwType{1}); + atom3->setInternalGlobalVDWType(VdwType{1}); atom1->setAtomType(0); atom2->setAtomType(1); diff --git a/tests/src/forceField/testBondForceField.cpp b/tests/src/forceField/testBondForceField.cpp index 51dcf5fe8..c61b9bf51 100644 --- a/tests/src/forceField/testBondForceField.cpp +++ b/tests/src/forceField/testBondForceField.cpp @@ -22,9 +22,8 @@ #include // for EXPECT_NEAR, Test, InitGoogleTest, RUN_ALL_TESTS -#include // for sqrt -#include // for size_t -#include // for shared_ptr, allocator +#include // for sqrt +#include // for shared_ptr, allocator #include "../potential/nonCoulomb/testForceFieldNonCoulomb.hpp" #include "atom.hpp" // for Atom @@ -38,6 +37,7 @@ #include "molecule.hpp" // for Molecule #include "physicalData.hpp" // for PhysicalData #include "simulationBox.hpp" // for SimulationBox +#include "strongTypes.hpp" namespace potential { @@ -57,11 +57,10 @@ TEST_F(TestBondForceField, calculateEnergyAndForces) auto coulombPotential = potential::CoulombShiftedPotential(10.0); auto nonCoulombPair = potential::LennardJonesPair( - static_cast(0), - static_cast(1), + ExtVdwType(0), + ExtVdwType(1), 5.0, - 2.0, - 4.0 + LJParams{.c6 = 2.0, .c12 = 4.0} ); setNonCoulombPairsMatrix( linearAlgebra::Matrix>(2, 2) @@ -80,8 +79,8 @@ TEST_F(TestBondForceField, calculateEnergyAndForces) atom2->setPosition({1.0, 2.0, 3.0}); atom1->setForce({0.0, 0.0, 0.0}); atom2->setForce({0.0, 0.0, 0.0}); - atom1->setInternalGlobalVDWType(0); - atom2->setInternalGlobalVDWType(1); + atom1->setInternalGlobalVDWType(VdwType{0}); + atom2->setInternalGlobalVDWType(VdwType{1}); atom1->setAtomType(0); atom2->setAtomType(1); atom1->setPartialCharge(1.0); diff --git a/tests/src/forceField/testDihedralForceField.cpp b/tests/src/forceField/testDihedralForceField.cpp index dcf9e2c54..2226c90a2 100644 --- a/tests/src/forceField/testDihedralForceField.cpp +++ b/tests/src/forceField/testDihedralForceField.cpp @@ -22,9 +22,8 @@ #include // for EXPECT_NEAR, Test, InitGoogleTest, RUN_ALL_TESTS -#include // for sqrt -#include // for size_t -#include // for shared_ptr, allocator +#include // for sqrt +#include // for shared_ptr, allocator #include "../potential/nonCoulomb/testForceFieldNonCoulomb.hpp" #include "atom.hpp" // for Atom @@ -59,11 +58,10 @@ TEST_F(TestDihedralForceField, calculateEnergyAndForces) auto coulombPotential = potential::CoulombShiftedPotential(20.0); auto nonCoulombPair = potential::LennardJonesPair( - static_cast(0), - static_cast(1), + ExtVdwType(0), + ExtVdwType(1), 15.0, - 2.0, - 4.0 + LJParams{.c6 = 2.0, .c12 = 4.0} ); setNonCoulombPairsMatrix( linearAlgebra::Matrix>(2, 2) @@ -90,10 +88,10 @@ TEST_F(TestDihedralForceField, calculateEnergyAndForces) atom3->setForce({0.0, 0.0, 0.0}); atom4->setForce({0.0, 0.0, 0.0}); - atom1->setInternalGlobalVDWType(0); - atom2->setInternalGlobalVDWType(1); - atom3->setInternalGlobalVDWType(0); - atom4->setInternalGlobalVDWType(1); + atom1->setInternalGlobalVDWType(VdwType{0}); + atom2->setInternalGlobalVDWType(VdwType{1}); + atom3->setInternalGlobalVDWType(VdwType{0}); + atom4->setInternalGlobalVDWType(VdwType{1}); atom1->setAtomType(0); atom2->setAtomType(1); diff --git a/tests/src/forceField/testForceField.cpp b/tests/src/forceField/testForceField.cpp index ee5d5a244..cc5651c47 100644 --- a/tests/src/forceField/testForceField.cpp +++ b/tests/src/forceField/testForceField.cpp @@ -22,8 +22,7 @@ #include // for Test, CmpHelperNE, TestInfo -#include // for size_t -#include // for shared_ptr, allocator +#include // for shared_ptr, allocator #include "../potential/nonCoulomb/testForceFieldNonCoulomb.hpp" #include "angleForceField.hpp" // for AngleForceField @@ -194,11 +193,10 @@ TEST_F(TestForceField, calculateBondedInteractions) auto coulombPotential = potential::CoulombShiftedPotential(20.0); auto nonCoulombPair = potential::LennardJonesPair( - static_cast(0), - static_cast(1), + ExtVdwType(0), + ExtVdwType(1), 15.0, - 2.0, - 4.0 + LJParams{.c6 = 2.0, .c12 = 4.0} ); setNonCoulombPairsMatrix( linearAlgebra::Matrix>(2, 2) @@ -225,10 +223,10 @@ TEST_F(TestForceField, calculateBondedInteractions) atom3->setForce({0.0, 0.0, 0.0}); atom4->setForce({0.0, 0.0, 0.0}); - atom1->setInternalGlobalVDWType(0); - atom2->setInternalGlobalVDWType(1); - atom3->setInternalGlobalVDWType(0); - atom4->setInternalGlobalVDWType(1); + atom1->setInternalGlobalVDWType(VdwType{0}); + atom2->setInternalGlobalVDWType(VdwType{1}); + atom3->setInternalGlobalVDWType(VdwType{0}); + atom4->setInternalGlobalVDWType(VdwType{1}); atom1->setAtomType(0); atom2->setAtomType(1); @@ -315,11 +313,10 @@ TEST_F(TestForceField, correctLinker) auto coulombPotential = potential::CoulombShiftedPotential(10.0); auto nonCoulombPair = potential::LennardJonesPair( - static_cast(0), - static_cast(1), + ExtVdwType(0), + ExtVdwType(1), 5.0, - 2.0, - 4.0 + LJParams{.c6 = 2.0, .c12 = 4.0} ); setNonCoulombPairsMatrix( linearAlgebra::Matrix>(2, 2) @@ -333,8 +330,8 @@ TEST_F(TestForceField, correctLinker) atom1->setForce({0.0, 0.0, 0.0}); atom2->setForce({0.0, 0.0, 0.0}); - atom1->setInternalGlobalVDWType(0); - atom2->setInternalGlobalVDWType(1); + atom1->setInternalGlobalVDWType(VdwType{0}); + atom2->setInternalGlobalVDWType(VdwType{1}); atom1->setAtomType(0); atom2->setAtomType(1); atom1->setPartialCharge(1.0); diff --git a/tests/src/input/parameterFileReader/CMakeLists.txt b/tests/src/input/parameterFileReader/CMakeLists.txt index a56b2a9c6..d9a3aefb1 100644 --- a/tests/src/input/parameterFileReader/CMakeLists.txt +++ b/tests/src/input/parameterFileReader/CMakeLists.txt @@ -23,6 +23,7 @@ foreach(source_file ${source_files}) parameterFileReader gtest gmock + testUtils pq_test_main ) add_test( diff --git a/tests/src/input/parameterFileReader/testNonCoulombicTypesSection.cpp b/tests/src/input/parameterFileReader/testNonCoulombicTypesSection.cpp index a57874c40..80e1a3359 100644 --- a/tests/src/input/parameterFileReader/testNonCoulombicTypesSection.cpp +++ b/tests/src/input/parameterFileReader/testNonCoulombicTypesSection.cpp @@ -22,15 +22,17 @@ #include // for EXPECT_EQ, TestInfo (ptr only) -#include "buckinghamPair.hpp" // for BuckinghamPair -#include "engine.hpp" // for Engine -#include "exceptions.hpp" // for ParameterFileException -#include "forceFieldNonCoulomb.hpp" // for ForceFieldNonCoulomb -#include "gtest/gtest.h" // for Message, TestPartResult, tes... -#include "lennardJonesPair.hpp" // for LennardJonesPair -#include "morsePair.hpp" // for MorsePair -#include "nonCoulombicsSection.hpp" // for NonCoulombicsSection -#include "potentialSettings.hpp" // for PotentialSettings +#include "buckinghamPair.hpp" // for BuckinghamPair +#include "engine.hpp" // for Engine +#include "exceptions.hpp" // for ParameterFileException +#include "forceFieldNonCoulomb.hpp" // for ForceFieldNonCoulomb +#include "gtest/gtest.h" // for Message, TestPartResult, tes... +#include "lennardJonesPair.hpp" // for LennardJonesPair +#include "morsePair.hpp" // for MorsePair +#include "nonCoulombicsSection.hpp" // for NonCoulombicsSection +#include "potentialSettings.hpp" // for PotentialSettings +#include "strongTypes.hpp" +#include "testNonCoulombPairUtils.hpp" #include "testParameterFileSection.hpp" // for TestParameterFileSection #include "throwWithMessage.hpp" // for ASSERT_THROW_MSG @@ -53,10 +55,10 @@ TEST_F(TestParameterFileSection, processSectionLennardJones) const auto *pairVector = potential.getNonCoulombPairsVector()[0].get(); const auto *pair = dynamic_cast(pairVector); - EXPECT_EQ(pair->getVanDerWaalsType1(), 0); - EXPECT_EQ(pair->getVanDerWaalsType2(), 1); - EXPECT_EQ(pair->getC6(), 1.22); - EXPECT_EQ(pair->getC12(), 234.3); + EXPECT_EQ(pair->getVanDerWaalsType1(), ExtVdwType{0}); + EXPECT_EQ(pair->getVanDerWaalsType2(), ExtVdwType{1}); + EXPECT_EQ(TestLJPairUtils::params(pair).c6, 1.22); + EXPECT_EQ(TestLJPairUtils::params(pair).c12, 234.3); EXPECT_EQ(pair->getRadialCutOff(), 324.3); lineElements = {"0", "1", "1.22", "234.3"}; @@ -65,10 +67,10 @@ TEST_F(TestParameterFileSection, processSectionLennardJones) const auto *pairVector2 = potential.getNonCoulombPairsVector()[1].get(); auto *pair2 = dynamic_cast(pairVector2); - EXPECT_EQ(pair2->getVanDerWaalsType1(), 0); - EXPECT_EQ(pair2->getVanDerWaalsType2(), 1); - EXPECT_EQ(pair2->getC6(), 1.22); - EXPECT_EQ(pair2->getC12(), 234.3); + EXPECT_EQ(pair2->getVanDerWaalsType1(), ExtVdwType{0}); + EXPECT_EQ(pair2->getVanDerWaalsType2(), ExtVdwType{1}); + EXPECT_EQ(TestLJPairUtils::params(pair2).c6, 1.22); + EXPECT_EQ(TestLJPairUtils::params(pair2).c12, 234.3); EXPECT_EQ(pair2->getRadialCutOff(), 12.5); lineElements = {"1", "2", "1.0", "0", "2", "3.3"}; @@ -92,11 +94,11 @@ TEST_F(TestParameterFileSection, processSectionBuckingham) const auto *pairVector = potential.getNonCoulombPairsVector()[0].get(); const auto *pair = dynamic_cast(pairVector); - EXPECT_EQ(pair->getVanDerWaalsType1(), 0); - EXPECT_EQ(pair->getVanDerWaalsType2(), 1); - EXPECT_EQ(pair->getA(), 1.22); - EXPECT_EQ(pair->getDRho(), 234.3); - EXPECT_EQ(pair->getC6(), 324.3); + EXPECT_EQ(pair->getVanDerWaalsType1(), ExtVdwType{0}); + EXPECT_EQ(pair->getVanDerWaalsType2(), ExtVdwType{1}); + EXPECT_EQ(TestBuckinghamPairUtils::params(pair).scaling, 1.22); + EXPECT_EQ(TestBuckinghamPairUtils::params(pair).dRho, 234.3); + EXPECT_EQ(TestBuckinghamPairUtils::params(pair).c6, 324.3); EXPECT_EQ(pair->getRadialCutOff(), 435.0); lineElements = {"0", "1", "1.22", "234.3", "324.3"}; @@ -105,11 +107,11 @@ TEST_F(TestParameterFileSection, processSectionBuckingham) const auto *pairVector2 = potential.getNonCoulombPairsVector()[1].get(); const auto *pair2 = dynamic_cast(pairVector2); - EXPECT_EQ(pair2->getVanDerWaalsType1(), 0); - EXPECT_EQ(pair2->getVanDerWaalsType2(), 1); - EXPECT_EQ(pair2->getA(), 1.22); - EXPECT_EQ(pair2->getDRho(), 234.3); - EXPECT_EQ(pair2->getC6(), 324.3); + EXPECT_EQ(pair2->getVanDerWaalsType1(), ExtVdwType{0}); + EXPECT_EQ(pair2->getVanDerWaalsType2(), ExtVdwType{1}); + EXPECT_EQ(TestBuckinghamPairUtils::params(pair2).scaling, 1.22); + EXPECT_EQ(TestBuckinghamPairUtils::params(pair2).dRho, 234.3); + EXPECT_EQ(TestBuckinghamPairUtils::params(pair2).c6, 324.3); EXPECT_EQ(pair2->getRadialCutOff(), 12.5); lineElements = {"1", "2", "1.0", "0", "2", "3.3", "345"}; @@ -133,11 +135,13 @@ TEST_F(TestParameterFileSection, processSectionMorse) auto *pair = dynamic_cast( potential.getNonCoulombPairsVector()[0].get() ); - EXPECT_EQ(pair->getVanDerWaalsType1(), 0); - EXPECT_EQ(pair->getVanDerWaalsType2(), 1); - EXPECT_EQ(pair->getDissociationEnergy(), 1.22); - EXPECT_EQ(pair->getWellWidth(), 234.3); - EXPECT_EQ(pair->getEquilibriumDistance(), 324.3); + EXPECT_EQ(pair->getVanDerWaalsType1(), ExtVdwType{0}); + EXPECT_EQ(pair->getVanDerWaalsType2(), ExtVdwType{1}); + + const auto &morseParams = TestMorsePairUtils::params(pair); + EXPECT_EQ(morseParams.dissociationEnergy, 1.22); + EXPECT_EQ(morseParams.wellWidth, 234.3); + EXPECT_EQ(morseParams.equilibriumDistance, 324.3); EXPECT_EQ(pair->getRadialCutOff(), 435.0); lineElements = {"0", "1", "1.22", "234.3", "324.3"}; @@ -146,11 +150,13 @@ TEST_F(TestParameterFileSection, processSectionMorse) auto *pair2 = dynamic_cast( potential.getNonCoulombPairsVector()[1].get() ); - EXPECT_EQ(pair2->getVanDerWaalsType1(), 0); - EXPECT_EQ(pair2->getVanDerWaalsType2(), 1); - EXPECT_EQ(pair2->getDissociationEnergy(), 1.22); - EXPECT_EQ(pair2->getWellWidth(), 234.3); - EXPECT_EQ(pair2->getEquilibriumDistance(), 324.3); + EXPECT_EQ(pair2->getVanDerWaalsType1(), ExtVdwType{0}); + EXPECT_EQ(pair2->getVanDerWaalsType2(), ExtVdwType{1}); + + const auto &morseParams2 = TestMorsePairUtils::params(pair2); + EXPECT_EQ(morseParams2.dissociationEnergy, 1.22); + EXPECT_EQ(morseParams2.wellWidth, 234.3); + EXPECT_EQ(morseParams2.equilibriumDistance, 324.3); EXPECT_EQ(pair2->getRadialCutOff(), 12.5); lineElements = {"1", "2", "1.0", "0", "2", "3.3", "345"}; diff --git a/tests/src/input/testGuffDatReader.cpp b/tests/src/input/testGuffDatReader.cpp index d25aeb68c..894bbddc6 100644 --- a/tests/src/input/testGuffDatReader.cpp +++ b/tests/src/input/testGuffDatReader.cpp @@ -30,8 +30,8 @@ #include // for string, basic_string, char_traits #include // for vector -#include "buckinghamPair.hpp" // for BuckinghamPair -#include "constants/internalConversionFactors.hpp" // for _COULOMB_PREFACTOR_ +#include "buckinghamPair.hpp" // for BuckinghamPair +#include "constants.hpp" // for _COULOMB_PREFACTOR_ #include "defaults.hpp" // for _NUMBER_OF_GUFF_ENTRIES_ #include "engine.hpp" // for Engine #include "exceptions.hpp" // for GuffDatException, UserInputException @@ -42,7 +42,9 @@ #include "morsePair.hpp" // for MorsePair #include "potentialSettings.hpp" // for PotentialSettings, string #include "settings.hpp" // for Settings -#include "throwWithMessage.hpp" // for EXPECT_THROW_MSG +#include "strongTypes.hpp" +#include "testNonCoulombPairUtils.hpp" +#include "throwWithMessage.hpp" // for EXPECT_THROW_MSG using namespace input::guffdat; using namespace potential; @@ -182,12 +184,15 @@ TEST_F(TestGuffDatReader, parseLine) _engine->getPotential()->getNonCoulombPotential() ); const auto &pair = dynamic_cast( - *potential.getNonCoulPair({1, 2, 1, 0}).get() + *potential.getNonCoulPair({1, 2, 1, 0}, {VdwType{1}, VdwType{0}}).get() ); EXPECT_EQ( pair, - LennardJonesPair(PotentialSettings::getCoulombRadiusCutOff(), 2.0, 3.0) + LennardJonesPair( + PotentialSettings::getCoulombRadiusCutOff(), + LJParams{2.0, 3.0} + ) ); } @@ -197,14 +202,15 @@ TEST_F(TestGuffDatReader, addLennardJonesPair) _guffDatReader->addLennardJonesPair(1, 2, 0, 0, {1.0, 2.0, 3.0}, 10.0); - const auto &pair = - dynamic_cast(*(_engine->getPotential() - ->getNonCoulombPotential() - .getNonCoulPair({1, 2, 0, 0}) - .get())); + const auto &pair = dynamic_cast( + *(_engine->getPotential() + ->getNonCoulombPotential() + .getNonCoulPair({1, 2, 0, 0}, {VdwType{0}, VdwType{0}}) + .get()) + ); - EXPECT_EQ(pair.getC6(), 1.0); - EXPECT_EQ(pair.getC12(), 3.0); + EXPECT_EQ(TestLJPairUtils::params(&pair).c6, 1.0); + EXPECT_EQ(TestLJPairUtils::params(&pair).c12, 3.0); EXPECT_EQ(pair.getRadialCutOff(), 10.0); // FIXME: Does not work on macos using EXPECT_EQ @@ -218,11 +224,12 @@ TEST_F(TestGuffDatReader, addLennardJonesPair) 6.0 / ::pow(10.0, 7) + 12.0 * 3.0 / ::pow(10.0, 13) ); - const auto &pair2 = - dynamic_cast(*(_engine->getPotential() - ->getNonCoulombPotential() - .getNonCoulPair({2, 1, 0, 0}) - .get())); + const auto &pair2 = dynamic_cast( + *(_engine->getPotential() + ->getNonCoulombPotential() + .getNonCoulPair({2, 1, 0, 0}, {VdwType{0}, VdwType{1}}) + .get()) + ); EXPECT_EQ(pair, pair2); } @@ -233,15 +240,16 @@ TEST_F(TestGuffDatReader, addBuckinghamPair) _guffDatReader->addBuckinghamPair(1, 2, 0, 0, {1.0, 2.0, 3.0}, 10.0); - const auto &pair = - dynamic_cast(*(_engine->getPotential() - ->getNonCoulombPotential() - .getNonCoulPair({1, 2, 0, 0}) - .get())); + const auto &pair = dynamic_cast( + *(_engine->getPotential() + ->getNonCoulombPotential() + .getNonCoulPair({1, 2, 0, 0}, {VdwType{0}, VdwType{0}}) + .get()) + ); - EXPECT_EQ(pair.getA(), 1.0); - EXPECT_EQ(pair.getDRho(), 2.0); - EXPECT_EQ(pair.getC6(), 3.0); + EXPECT_EQ(TestBuckinghamPairUtils::params(&pair).scaling, 1.0); + EXPECT_EQ(TestBuckinghamPairUtils::params(&pair).dRho, 2.0); + EXPECT_EQ(TestBuckinghamPairUtils::params(&pair).c6, 3.0); EXPECT_EQ(pair.getRadialCutOff(), 10.0); EXPECT_EQ( pair.getEnergyCutOff(), @@ -252,11 +260,12 @@ TEST_F(TestGuffDatReader, addBuckinghamPair) -2.0 * ::exp(10.0 * 2.0) + 6.0 * 3.0 / ::pow(10.0, 7) ); - const auto &pair2 = - dynamic_cast(*(_engine->getPotential() - ->getNonCoulombPotential() - .getNonCoulPair({2, 1, 0, 0}) - .get())); + const auto &pair2 = dynamic_cast( + *(_engine->getPotential() + ->getNonCoulombPotential() + .getNonCoulPair({2, 1, 0, 0}, {VdwType{0}, VdwType{1}}) + .get()) + ); EXPECT_EQ(pair, pair2); } @@ -267,15 +276,18 @@ TEST_F(TestGuffDatReader, addMorsePair) _guffDatReader->addMorsePair(1, 2, 0, 0, {1.0, 2.0, 3.0}, 10.0); - const auto &pair = - dynamic_cast(*(_engine->getPotential() - ->getNonCoulombPotential() - .getNonCoulPair({1, 2, 0, 0}) - .get())); + const auto &pair = dynamic_cast( + *(_engine->getPotential() + ->getNonCoulombPotential() + .getNonCoulPair({1, 2, 0, 0}, {VdwType{0}, VdwType{0}}) + .get()) + ); + + const auto ¶ms = TestMorsePairUtils::params(&pair); - EXPECT_EQ(pair.getDissociationEnergy(), 1.0); - EXPECT_EQ(pair.getWellWidth(), 2.0); - EXPECT_EQ(pair.getEquilibriumDistance(), 3.0); + EXPECT_EQ(params.dissociationEnergy, 1.0); + EXPECT_EQ(params.wellWidth, 2.0); + EXPECT_EQ(params.equilibriumDistance, 3.0); EXPECT_EQ(pair.getRadialCutOff(), 10.0); EXPECT_EQ( pair.getEnergyCutOff(), @@ -287,11 +299,12 @@ TEST_F(TestGuffDatReader, addMorsePair) 2.0 ); - const auto &pair2 = - dynamic_cast(*(_engine->getPotential() - ->getNonCoulombPotential() - .getNonCoulPair({2, 1, 0, 0}) - .get())); + const auto &pair2 = dynamic_cast( + *(_engine->getPotential() + ->getNonCoulombPotential() + .getNonCoulPair({2, 1, 0, 0}, {VdwType{0}, VdwType{1}}) + .get()) + ); EXPECT_EQ(pair, pair2); } @@ -305,18 +318,19 @@ TEST_F(TestGuffDatReader, addGuffPair) 2, 0, 0, - {1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1}, + {1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2}, 10.0 ); - const auto &pair = - dynamic_cast(*(_engine->getPotential() - ->getNonCoulombPotential() - .getNonCoulPair({1, 2, 0, 0}) - .get())); + const auto &pair = dynamic_cast( + *(_engine->getPotential() + ->getNonCoulombPotential() + .getNonCoulPair({1, 2, 0, 0}, {VdwType{0}, VdwType{0}}) + .get()) + ); EXPECT_THAT( - pair.getCoefficients(), + TestGuffPairUtils::coeffs(&pair), testing::ElementsAre( 1, 2, @@ -339,22 +353,22 @@ TEST_F(TestGuffDatReader, addGuffPair) 1, 1, 1, - 2, - 1 + 2 ) ); EXPECT_EQ(pair.getRadialCutOff(), 10.0); EXPECT_EQ(pair.getEnergyCutOff(), 3.0121946291700612e+35); EXPECT_EQ(pair.getForceCutOff(), -5.4219503325061099e+36); - const auto &pair2 = - dynamic_cast(*(_engine->getPotential() - ->getNonCoulombPotential() - .getNonCoulPair({2, 1, 0, 0}) - .get())); + const auto &pair2 = dynamic_cast( + *(_engine->getPotential() + ->getNonCoulombPotential() + .getNonCoulPair({2, 1, 0, 0}, {VdwType{0}, VdwType{1}}) + .get()) + ); EXPECT_THAT( - pair2.getCoefficients(), + TestGuffPairUtils::coeffs(&pair2), testing::ElementsAre( 1, 2, @@ -377,70 +391,109 @@ TEST_F(TestGuffDatReader, addGuffPair) 1, 1, 1, - 2, - 1 + 2 ) ); EXPECT_EQ(pair2.getRadialCutOff(), 10.0); - EXPECT_EQ(pair.getEnergyCutOff(), 3.0121946291700612e+35); - EXPECT_EQ(pair.getForceCutOff(), -5.4219503325061099e+36); + EXPECT_EQ(pair2.getEnergyCutOff(), 3.0121946291700612e+35); + EXPECT_EQ(pair2.getForceCutOff(), -5.4219503325061099e+36); } TEST_F(TestGuffDatReader, addNonCoulombPair) { const auto &guffCoefficients = - std::vector({1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, - 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1}); + std::array{ + 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2 + }; _guffDatReader->setupGuffMaps(); PotentialSettings::setNonCoulombType("lj"); - _guffDatReader->addNonCoulombPair(1, 2, 0, 0, guffCoefficients, 10.0); + _guffDatReader->addNonCoulombPair( + 1, + 2, + 0, + 0, + {guffCoefficients.begin(), guffCoefficients.end()}, + 10.0 + ); EXPECT_NO_THROW( - [[maybe_unused]] const auto &dummy = - dynamic_cast(*(_engine->getPotential() - ->getNonCoulombPotential() - .getNonCoulPair({2, 1, 0, 0}) - .get())) + [[maybe_unused]] const auto &dummy = dynamic_cast( + *(_engine->getPotential() + ->getNonCoulombPotential() + .getNonCoulPair({2, 1, 0, 0}, {VdwType{0}, VdwType{1}}) + .get()) + ) ); PotentialSettings::setNonCoulombType("buck"); - _guffDatReader->addNonCoulombPair(1, 2, 0, 0, guffCoefficients, 10.0); + _guffDatReader->addNonCoulombPair( + 1, + 2, + 0, + 0, + {guffCoefficients.begin(), guffCoefficients.end()}, + 10.0 + ); EXPECT_NO_THROW( - [[maybe_unused]] const auto &dummy = - dynamic_cast(*(_engine->getPotential() - ->getNonCoulombPotential() - .getNonCoulPair({2, 1, 0, 0}) - .get())) + [[maybe_unused]] const auto &dummy = dynamic_cast( + *(_engine->getPotential() + ->getNonCoulombPotential() + .getNonCoulPair({2, 1, 0, 0}, {VdwType{0}, VdwType{1}}) + .get()) + ) ); PotentialSettings::setNonCoulombType("morse"); - _guffDatReader->addNonCoulombPair(1, 2, 0, 0, guffCoefficients, 10.0); + _guffDatReader->addNonCoulombPair( + 1, + 2, + 0, + 0, + {guffCoefficients.begin(), guffCoefficients.end()}, + 10.0 + ); EXPECT_NO_THROW( - [[maybe_unused]] const auto &dummy = - dynamic_cast(*(_engine->getPotential() - ->getNonCoulombPotential() - .getNonCoulPair({2, 1, 0, 0}) - .get())) + [[maybe_unused]] const auto &dummy = dynamic_cast( + *(_engine->getPotential() + ->getNonCoulombPotential() + .getNonCoulPair({2, 1, 0, 0}, {VdwType{0}, VdwType{1}}) + .get()) + ) ); PotentialSettings::setNonCoulombType("guff"); - _guffDatReader->addNonCoulombPair(1, 2, 0, 0, guffCoefficients, 10.0); + _guffDatReader->addNonCoulombPair( + 1, + 2, + 0, + 0, + {guffCoefficients.begin(), guffCoefficients.end()}, + 10.0 + ); EXPECT_NO_THROW( - [[maybe_unused]] const auto &dummy = - dynamic_cast(*(_engine->getPotential() - ->getNonCoulombPotential() - .getNonCoulPair({2, 1, 0, 0}) - .get())) + [[maybe_unused]] const auto &dummy = dynamic_cast( + *(_engine->getPotential() + ->getNonCoulombPotential() + .getNonCoulPair({2, 1, 0, 0}, {VdwType{1}, VdwType{2}}) + .get()) + ) ); PotentialSettings::setNonCoulombType("lj_9_12"); EXPECT_THROW_MSG( - _guffDatReader->addNonCoulombPair(1, 2, 0, 0, guffCoefficients, 10.0), + _guffDatReader->addNonCoulombPair( + 1, + 2, + 0, + 0, + {guffCoefficients.begin(), guffCoefficients.end()}, + 10.0 + ), UserInputException, std::format( "Invalid nonCoulombic type {} given", diff --git a/tests/src/intraNonBonded/testIntraNonBonded.cpp b/tests/src/intraNonBonded/testIntraNonBonded.cpp index 1d49d3852..32ae6fb36 100644 --- a/tests/src/intraNonBonded/testIntraNonBonded.cpp +++ b/tests/src/intraNonBonded/testIntraNonBonded.cpp @@ -22,10 +22,9 @@ #include // for Test, EXPECT_EQ, TestInfo -#include // for size_t -#include // for format -#include // for shared_ptr, allocator, make_shared -#include // for vector +#include // for format +#include // for shared_ptr, allocator, make_shared +#include // for vector #include "../potential/nonCoulomb/testForceFieldNonCoulomb.hpp" #include "atom.hpp" // for Atom @@ -42,7 +41,8 @@ #include "physicalData.hpp" // for PhysicalData #include "potentialSettings.hpp" // for PotentialSettings #include "simulationBox.hpp" // for SimulationBox -#include "throwWithMessage.hpp" // for EXPECT_THROW_MSG +#include "strongTypes.hpp" +#include "throwWithMessage.hpp" // for EXPECT_THROW_MSG namespace potential { @@ -182,8 +182,8 @@ TEST_F(TestIntraNonBonded, calculate) atom2->setPosition({0.0, 0.0, 11.0}); atom1->setForce({0.0, 0.0, 0.0}); atom2->setForce({0.0, 0.0, 0.0}); - atom1->setInternalGlobalVDWType(0); - atom2->setInternalGlobalVDWType(1); + atom1->setInternalGlobalVDWType(VdwType{0}); + atom2->setInternalGlobalVDWType(VdwType{1}); atom1->setAtomType(0); atom2->setAtomType(1); atom1->setPartialCharge(0.5); @@ -206,11 +206,10 @@ TEST_F(TestIntraNonBonded, calculate) ); auto nonCoulombPair = potential::LennardJonesPair( - static_cast(0), - static_cast(1), + ExtVdwType(0), + ExtVdwType(1), 10.0, - 2.0, - 3.0 + LJParams{.c6 = 2.0, .c12 = 3.0} ); setNonCoulombPairsMatrix(0, 1, nonCoulombPair); setNonCoulombPairsMatrix(1, 0, nonCoulombPair); diff --git a/tests/src/intraNonBonded/testIntraNonBondedMap.cpp b/tests/src/intraNonBonded/testIntraNonBondedMap.cpp index 9dd1e3dc5..8f779a074 100644 --- a/tests/src/intraNonBonded/testIntraNonBondedMap.cpp +++ b/tests/src/intraNonBonded/testIntraNonBondedMap.cpp @@ -22,9 +22,8 @@ #include // for Test, EXPECT_NEAR, InitGoogleTest, RUN_ALL. -#include // for size_t -#include // for shared_ptr, allocator -#include // for vector +#include // for shared_ptr, allocator +#include // for vector #include "../potential/nonCoulomb/testForceFieldNonCoulomb.hpp" #include "atom.hpp" // for Atom @@ -39,6 +38,7 @@ #include "physicalData.hpp" // for PhysicalData #include "potentialSettings.hpp" // for PotentialSettings #include "simulationBox.hpp" // for SimulationBox +#include "strongTypes.hpp" namespace potential { @@ -64,8 +64,8 @@ TEST_F(TestIntraNonBondedMap, calculateSingleInteractionAndCalculate) atom2->setPosition({0.0, 0.0, 11.0}); atom1->setForce({0.0, 0.0, 0.0}); atom2->setForce({0.0, 0.0, 0.0}); - atom1->setInternalGlobalVDWType(0); - atom2->setInternalGlobalVDWType(1); + atom1->setInternalGlobalVDWType(VdwType{0}); + atom2->setInternalGlobalVDWType(VdwType{1}); atom1->setAtomType(0); atom2->setAtomType(1); atom1->setPartialCharge(0.5); @@ -88,11 +88,10 @@ TEST_F(TestIntraNonBondedMap, calculateSingleInteractionAndCalculate) ); auto nonCoulombPair = potential::LennardJonesPair( - static_cast(0), - static_cast(1), + ExtVdwType(0), + ExtVdwType(1), 10.0, - 2.0, - 3.0 + LJParams{.c6 = 2.0, .c12 = 3.0} ); setNonCoulombPairsMatrix(0, 1, nonCoulombPair); setNonCoulombPairsMatrix(1, 0, nonCoulombPair); diff --git a/tests/src/molsys/testSimulationBox.cpp b/tests/src/molsys/testSimulationBox.cpp index a78fcdea1..93bbe3aa3 100644 --- a/tests/src/molsys/testSimulationBox.cpp +++ b/tests/src/molsys/testSimulationBox.cpp @@ -31,8 +31,9 @@ #include "exceptions.hpp" // for ManostatException, RstFileException #include "gtest/gtest.h" // for Message, TestPartResult, AssertionRe... #include "potentialSettings.hpp" // for PotentialSettings -#include "throwWithMessage.hpp" // for throwWithMessage -#include "vectorNear.hpp" // for EXPECT_VECTOR_NEAR +#include "strongTypes.hpp" +#include "throwWithMessage.hpp" // for throwWithMessage +#include "vectorNear.hpp" // for EXPECT_VECTOR_NEAR /** * @brief tests numberOfAtoms function @@ -236,12 +237,12 @@ TEST_F(TestSimulationBox, setupExternalToInternalGlobalVdwTypesMap) molsys::MoleculeType molecule1(1); molsys::MoleculeType molecule2(2); - molecule1.addExternalGlobalVDWType(1); - molecule1.addExternalGlobalVDWType(3); - molecule1.addExternalGlobalVDWType(5); + molecule1.addExternalGlobalVDWType(ExtVdwType{1}); + molecule1.addExternalGlobalVDWType(ExtVdwType{3}); + molecule1.addExternalGlobalVDWType(ExtVdwType{5}); - molecule2.addExternalGlobalVDWType(3); - molecule2.addExternalGlobalVDWType(5); + molecule2.addExternalGlobalVDWType(ExtVdwType{3}); + molecule2.addExternalGlobalVDWType(ExtVdwType{5}); simulationBox.addMoleculeType(molecule1); simulationBox.addMoleculeType(molecule2); @@ -251,13 +252,16 @@ TEST_F(TestSimulationBox, setupExternalToInternalGlobalVdwTypesMap) EXPECT_EQ(simulationBox.getExternalGlobalVdwTypes().size(), 3); EXPECT_EQ( simulationBox.getExternalGlobalVdwTypes(), - std::vector({1, 3, 5}) + std::vector({ExtVdwType{1}, ExtVdwType{3}, ExtVdwType{5}}) ); - EXPECT_EQ(simulationBox.getExternalToInternalGlobalVDWTypes().size(), 3); - EXPECT_EQ(simulationBox.getExternalToInternalGlobalVDWTypes().at(1), 0); - EXPECT_EQ(simulationBox.getExternalToInternalGlobalVDWTypes().at(3), 1); - EXPECT_EQ(simulationBox.getExternalToInternalGlobalVDWTypes().at(5), 2); + const auto &externalToInternalMap = + simulationBox.getExternalToInternalGlobalVDWTypes(); + + EXPECT_EQ(externalToInternalMap.size(), 3); + EXPECT_EQ(externalToInternalMap.at(ExtVdwType{1}), VdwType{0}); + EXPECT_EQ(externalToInternalMap.at(ExtVdwType{3}), VdwType{1}); + EXPECT_EQ(externalToInternalMap.at(ExtVdwType{5}), VdwType{2}); } /** @@ -564,14 +568,14 @@ TEST_F(TestSimulationBox, assignsInternalVdwTypesToAtoms) { molsys::SimulationBox simBox; molsys::MoleculeType type(1); - type.addExternalGlobalVDWType(4); - type.addExternalGlobalVDWType(9); + type.addExternalGlobalVDWType(ExtVdwType{4}); + type.addExternalGlobalVDWType(ExtVdwType{9}); simBox.addMoleculeType(type); auto atom1 = std::make_shared(); auto atom2 = std::make_shared(); - atom1->setExternalGlobalVDWType(4); - atom2->setExternalGlobalVDWType(9); + atom1->setExternalGlobalVDWType(ExtVdwType{4}); + atom2->setExternalGlobalVDWType(ExtVdwType{9}); molsys::Molecule molecule(1); molecule.setNumberOfAtoms(2); @@ -581,8 +585,9 @@ TEST_F(TestSimulationBox, assignsInternalVdwTypesToAtoms) simBox.setupExternalToInternalGlobalVdwTypesMap(); - EXPECT_EQ(simBox.getMolecule(0).getAtom(0).getInternalGlobalVDWType(), 0); - EXPECT_EQ(simBox.getMolecule(0).getAtom(1).getInternalGlobalVDWType(), 1); + auto &moleculeResult = simBox.getMolecule(0); + EXPECT_EQ(moleculeResult.getAtom(0).getInternalGlobalVDWType(), VdwType{0}); + EXPECT_EQ(moleculeResult.getAtom(1).getInternalGlobalVDWType(), VdwType{1}); } TEST_F(TestSimulationBox, forceMetricsAndAtomStateUpdates) diff --git a/tests/src/potential/nonCoulomb/CMakeLists.txt b/tests/src/potential/nonCoulomb/CMakeLists.txt index 0bc9793e4..c77c4c75b 100644 --- a/tests/src/potential/nonCoulomb/CMakeLists.txt +++ b/tests/src/potential/nonCoulomb/CMakeLists.txt @@ -18,6 +18,7 @@ foreach(source_file ${source_files}) PRIVATE PQ_input potential + nonCoulombPotential gtest gmock pq_test_main diff --git a/tests/src/potential/nonCoulomb/testBuckinghamPair.cpp b/tests/src/potential/nonCoulomb/testBuckinghamPair.cpp index f5b225a06..b97eb8d83 100644 --- a/tests/src/potential/nonCoulomb/testBuckinghamPair.cpp +++ b/tests/src/potential/nonCoulomb/testBuckinghamPair.cpp @@ -22,9 +22,8 @@ #include // for Test, CmpHelperFloatingPointEQ, EXPECT_EQ -#include // for pow, exp -#include // for size_t -#include // for vector +#include // for pow, exp +#include // for vector #include "buckinghamPair.hpp" // for BuckinghamPair #include "gtest/gtest.h" // for AssertionResult, Message, TestPartResult @@ -37,38 +36,71 @@ using namespace potential; */ TEST(TestBuckinghamPair, equalsOperator) { - const size_t vdwType1 = 0; - const size_t vdwType2 = 1; - const size_t vdwType3 = 2; - const auto nonCoulombPair1 = - BuckinghamPair(vdwType1, vdwType2, 1.0, 2.0, 3.0, 4.0); - - const auto nonCoulombPair2 = - BuckinghamPair(vdwType1, vdwType2, 1.0, 2.0, 3.0, 4.0); + const ExtVdwType vdwType1{0}; + const ExtVdwType vdwType2{1}; + const ExtVdwType vdwType3{2}; + + const auto nonCoulombPair1 = BuckinghamPair( + vdwType1, + vdwType2, + 1.0, + BuckinghamParams{.scaling = 2.0, .dRho = 3.0, .c6 = 4.0} + ); + + const auto nonCoulombPair2 = BuckinghamPair( + vdwType1, + vdwType2, + 1.0, + BuckinghamParams{.scaling = 2.0, .dRho = 3.0, .c6 = 4.0} + ); EXPECT_TRUE(nonCoulombPair1 == nonCoulombPair2); - const auto nonCoulombPair3 = - BuckinghamPair(vdwType2, vdwType1, 1.0, 2.0, 3.0, 4.0); + const auto nonCoulombPair3 = BuckinghamPair( + vdwType2, + vdwType1, + 1.0, + BuckinghamParams{.scaling = 2.0, .dRho = 3.0, .c6 = 4.0} + ); EXPECT_TRUE(nonCoulombPair1 == nonCoulombPair3); - const auto nonCoulombPair4 = - BuckinghamPair(vdwType1, vdwType3, 1.0, 2.0, 3.0, 4.0); + const auto nonCoulombPair4 = BuckinghamPair( + vdwType1, + vdwType3, + 1.0, + BuckinghamParams{.scaling = 2.0, .dRho = 3.0, .c6 = 4.0} + ); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair4); - const auto nonCoulombPair5 = - BuckinghamPair(vdwType1, vdwType2, 2.0, 2.0, 3.0, 4.0); + const auto nonCoulombPair5 = BuckinghamPair( + vdwType1, + vdwType2, + 2.0, + BuckinghamParams{.scaling = 2.0, .dRho = 3.0, .c6 = 4.0} + ); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair5); - const auto nonCoulombPair6 = - BuckinghamPair(vdwType1, vdwType2, 1.0, 3.0, 3.0, 4.0); + const auto nonCoulombPair6 = BuckinghamPair( + vdwType1, + vdwType2, + 1.0, + BuckinghamParams{.scaling = 3.0, .dRho = 3.0, .c6 = 4.0} + ); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair6); - const auto nonCoulombPair7 = - BuckinghamPair(vdwType1, vdwType2, 1.0, 2.0, 4.0, 4.0); + const auto nonCoulombPair7 = BuckinghamPair( + vdwType1, + vdwType2, + 1.0, + BuckinghamParams{.scaling = 2.0, .dRho = 4.0, .c6 = 4.0} + ); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair7); - const auto nonCoulombPair8 = - BuckinghamPair(vdwType1, vdwType2, 1.0, 2.0, 3.0, 5.0); + const auto nonCoulombPair8 = BuckinghamPair( + vdwType1, + vdwType2, + 1.0, + BuckinghamParams{.scaling = 2.0, .dRho = 3.0, .c6 = 5.0} + ); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair8); } @@ -87,9 +119,11 @@ TEST(TestBuckinghamPair, calculateEnergyAndForces) rncCutoff, energyCutoff, forceCutoff, - coefficients[0], - coefficients[1], - coefficients[2] + BuckinghamParams{ + .scaling = coefficients[0], + .dRho = coefficients[1], + .c6 = coefficients[2] + } ); auto distance = 2.0; diff --git a/tests/src/potential/nonCoulomb/testForceFieldNonCoulomb.cpp b/tests/src/potential/nonCoulomb/testForceFieldNonCoulomb.cpp index 6627871a9..8c413bc21 100644 --- a/tests/src/potential/nonCoulomb/testForceFieldNonCoulomb.cpp +++ b/tests/src/potential/nonCoulomb/testForceFieldNonCoulomb.cpp @@ -24,8 +24,6 @@ #include // for Test, EXPECT_EQ, TestInfo -#include // for size_t -#include // for map #include // for make_shared, shared_ptr #include // for optional, nullopt #include // for move @@ -37,7 +35,8 @@ #include "lennardJonesPair.hpp" // for LennardJonesPair #include "matrix.hpp" // for Matrix #include "nonCoulombPair.hpp" // for NonCoulombPair -#include "throwWithMessage.hpp" // for EXPECT_THROW_MSG +#include "strongTypes.hpp" +#include "throwWithMessage.hpp" // for EXPECT_THROW_MSG TEST_F(TestNonCoulombPotentialFF, copyConstructorCopiesOwnedMatrix) { @@ -45,11 +44,10 @@ TEST_F(TestNonCoulombPotentialFF, copyConstructorCopiesOwnedMatrix) linearAlgebra::Matrix>(1) ); const auto pair = potential::LennardJonesPair( - static_cast(1), - static_cast(1), + ExtVdwType(1), + ExtVdwType(1), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ); setNonCoulombPairsMatrix(0, 0, pair); _nonCoulombPotential->setNonCoulombPairsVector( @@ -65,11 +63,10 @@ TEST_F(TestNonCoulombPotentialFF, copyConstructorCopiesOwnedMatrix) ); const auto replacement = potential::LennardJonesPair( - static_cast(1), - static_cast(1), + ExtVdwType(1), + ExtVdwType(1), 3.0, - 2.0, - 1.0 + LJParams{.c6 = 2.0, .c12 = 1.0} ); setNonCoulombPairsMatrix(*_nonCoulombPotential, 0, 0, replacement); EXPECT_NE( @@ -84,11 +81,10 @@ TEST_F(TestNonCoulombPotentialFF, copyAssignmentCopiesOwnedMatrix) linearAlgebra::Matrix>(1) ); const auto pair = potential::LennardJonesPair( - static_cast(1), - static_cast(1), + ExtVdwType(1), + ExtVdwType(1), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ); setNonCoulombPairsMatrix(0, 0, pair); @@ -96,7 +92,7 @@ TEST_F(TestNonCoulombPotentialFF, copyAssignmentCopiesOwnedMatrix) copy = *_nonCoulombPotential; const auto matrixElement = getNonCoulombPairsMatrix(copy)(0, 0); - const auto *self = © + const auto* self = © copy = *self; EXPECT_EQ(getNonCoulombPairsMatrix(copy)(0, 0), matrixElement); } @@ -107,11 +103,10 @@ TEST_F(TestNonCoulombPotentialFF, moveOperationsTransferOwnedMatrix) linearAlgebra::Matrix>(1) ); const auto pair = potential::LennardJonesPair( - static_cast(1), - static_cast(1), + ExtVdwType(1), + ExtVdwType(1), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ); setNonCoulombPairsMatrix(0, 0, pair); @@ -132,45 +127,38 @@ TEST_F(TestNonCoulombPotentialFF, determineInternalGlobalVdwTypes) { _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(5), + ExtVdwType(1), + ExtVdwType(5), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(2), + ExtVdwType(1), + ExtVdwType(2), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); - std::map externalToInternalTypes({{1, 0}, {2, 1}, {5, 2}}); + std::unordered_map externalToInternalTypes( + {{ExtVdwType{1}, VdwType{0}}, + {ExtVdwType{2}, VdwType{1}}, + {ExtVdwType{5}, VdwType{2}}} + ); _nonCoulombPotential->determineInternalGlobalVdwTypes( externalToInternalTypes ); - EXPECT_EQ( - _nonCoulombPotential->getNonCoulombPairsVector()[0]->getInternalType1(), - 0 - ); - EXPECT_EQ( - _nonCoulombPotential->getNonCoulombPairsVector()[0]->getInternalType2(), - 2 - ); - EXPECT_EQ( - _nonCoulombPotential->getNonCoulombPairsVector()[1]->getInternalType1(), - 0 - ); - EXPECT_EQ( - _nonCoulombPotential->getNonCoulombPairsVector()[1]->getInternalType2(), - 1 - ); + const auto& nonCoulPairsVec = + _nonCoulombPotential->getNonCoulombPairsVector(); + + EXPECT_EQ(nonCoulPairsVec[0]->getInternalType1(), VdwType{0}); + EXPECT_EQ(nonCoulPairsVec[0]->getInternalType2(), VdwType{2}); + EXPECT_EQ(nonCoulPairsVec[1]->getInternalType1(), VdwType{0}); + EXPECT_EQ(nonCoulPairsVec[1]->getInternalType2(), VdwType{1}); } /** @@ -180,23 +168,21 @@ TEST_F(TestNonCoulombPotentialFF, determineInternalGlobalVdwTypes) TEST_F(TestNonCoulombPotentialFF, fillDiagOfNonCoulPairsMatrix) { auto nonCoulombicPair1 = potential::LennardJonesPair( - static_cast(1), - static_cast(1), + ExtVdwType(1), + ExtVdwType(1), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ); - nonCoulombicPair1.setInternalType1(0); - nonCoulombicPair1.setInternalType2(0); + nonCoulombicPair1.setInternalType1(VdwType{0}); + nonCoulombicPair1.setInternalType2(VdwType{0}); auto nonCoulombicPair2 = potential::LennardJonesPair( - static_cast(9), - static_cast(9), + ExtVdwType(9), + ExtVdwType(9), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ); - nonCoulombicPair2.setInternalType1(9); - nonCoulombicPair2.setInternalType2(9); + nonCoulombicPair2.setInternalType1(VdwType{9}); + nonCoulombicPair2.setInternalType2(VdwType{9}); std::vector> diagonalElements = { std::make_shared(nonCoulombicPair1), @@ -207,10 +193,10 @@ TEST_F(TestNonCoulombPotentialFF, fillDiagOfNonCoulPairsMatrix) EXPECT_EQ(getNonCoulombPairsMatrix().rows(), 2); EXPECT_EQ(getNonCoulombPairsMatrix().cols(), 2); - EXPECT_EQ(getNonCoulombPairsMatrix()(0, 0)->getInternalType1(), 0); - EXPECT_EQ(getNonCoulombPairsMatrix()(0, 0)->getInternalType2(), 0); - EXPECT_EQ(getNonCoulombPairsMatrix()(1, 1)->getInternalType1(), 9); - EXPECT_EQ(getNonCoulombPairsMatrix()(1, 1)->getInternalType2(), 9); + EXPECT_EQ(getNonCoulombPairsMatrix()(0, 0)->getInternalType1(), VdwType{0}); + EXPECT_EQ(getNonCoulombPairsMatrix()(0, 0)->getInternalType2(), VdwType{0}); + EXPECT_EQ(getNonCoulombPairsMatrix()(1, 1)->getInternalType1(), VdwType{9}); + EXPECT_EQ(getNonCoulombPairsMatrix()(1, 1)->getInternalType2(), VdwType{9}); } /** @@ -225,34 +211,39 @@ TEST_F( { _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(5), + ExtVdwType(1), + ExtVdwType(5), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(2), + ExtVdwType(1), + ExtVdwType(2), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); // these two lines were already tested in // TestPotential_determineInternalGlobalVdwTypes - std::map externalToInternalTypes({{1, 0}, {2, 1}, {5, 2}}); + std::unordered_map externalToInternalTypes( + {{ExtVdwType{1}, VdwType{0}}, + {ExtVdwType{2}, VdwType{1}}, + {ExtVdwType{5}, VdwType{2}}} + ); _nonCoulombPotential->determineInternalGlobalVdwTypes( externalToInternalTypes ); auto nonCoulombicPair = - _nonCoulombPotential->findNonCoulPairByInternalTypes(0, 2); - EXPECT_EQ((*nonCoulombicPair)->getInternalType1(), 0); - EXPECT_EQ((*nonCoulombicPair)->getInternalType2(), 2); + _nonCoulombPotential->findNonCoulPairByInternalTypes( + VdwType{0}, + VdwType{2} + ); + EXPECT_EQ((*nonCoulombicPair)->getInternalType1(), VdwType{0}); + EXPECT_EQ((*nonCoulombicPair)->getInternalType2(), VdwType{2}); } /** @@ -267,32 +258,37 @@ TEST_F( { _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(5), + ExtVdwType(1), + ExtVdwType(5), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(2), + ExtVdwType(1), + ExtVdwType(2), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); // these two lines were already tested in // TestPotential_determineInternalGlobalVdwTypes - std::map externalToInternalTypes({{1, 0}, {2, 1}, {5, 2}}); + std::unordered_map externalToInternalTypes( + {{ExtVdwType{1}, VdwType{0}}, + {ExtVdwType{2}, VdwType{1}}, + {ExtVdwType{5}, VdwType{2}}} + ); _nonCoulombPotential->determineInternalGlobalVdwTypes( externalToInternalTypes ); auto nonCoulombicPair = - _nonCoulombPotential->findNonCoulPairByInternalTypes(0, 3); + _nonCoulombPotential->findNonCoulPairByInternalTypes( + VdwType{0}, + VdwType{3} + ); EXPECT_EQ(nonCoulombicPair, std::nullopt); } @@ -308,45 +304,51 @@ TEST_F( { _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(5), + ExtVdwType(1), + ExtVdwType(5), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(5), + ExtVdwType(1), + ExtVdwType(5), 2.0, - 5.0, - 1.0 + LJParams{.c6 = 5.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(2), + ExtVdwType(1), + ExtVdwType(2), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); // these two lines were already tested in // TestPotential_determineInternalGlobalVdwTypes - std::map externalToInternalTypes({{1, 0}, {2, 1}, {5, 2}}); + std::unordered_map externalToInternalTypes( + {{ExtVdwType{1}, VdwType{0}}, + {ExtVdwType{2}, VdwType{1}}, + {ExtVdwType{5}, VdwType{2}}} + ); _nonCoulombPotential->determineInternalGlobalVdwTypes( externalToInternalTypes ); EXPECT_THROW_MSG( [[maybe_unused]] const auto dummy = - _nonCoulombPotential->findNonCoulPairByInternalTypes(0, 2), + _nonCoulombPotential + ->findNonCoulPairByInternalTypes(VdwType{0}, VdwType{2}), exc::ParameterFileException, - "Non coulombic pair with global van der waals types 1 and 5 is defined " - "twice in the parameter file." + std::format( + "Non coulombic pair with global van der waals types {} and {} is " + "defined twice in the parameter file.", + ExtVdwType{1}.toString(), + ExtVdwType{5}.toString() + ) ); } @@ -362,17 +364,20 @@ TEST_F( { _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(5), + ExtVdwType(1), + ExtVdwType(5), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); // these two lines were already tested in // TestPotential_determineInternalGlobalVdwTypes - std::map externalToInternalTypes({{1, 0}, {2, 1}, {5, 2}}); + std::unordered_map externalToInternalTypes( + {{ExtVdwType{1}, VdwType{0}}, + {ExtVdwType{2}, VdwType{1}}, + {ExtVdwType{5}, VdwType{2}}} + ); _nonCoulombPotential->determineInternalGlobalVdwTypes( externalToInternalTypes ); @@ -400,35 +405,36 @@ TEST_F( { _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(2), + ExtVdwType(1), + ExtVdwType(2), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(5), + ExtVdwType(1), + ExtVdwType(5), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(2), - static_cast(5), + ExtVdwType(2), + ExtVdwType(5), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); // these two lines were already tested in // TestPotential_determineInternalGlobalVdwTypes - std::map externalToInternalTypes({{1, 0}, {2, 1}, {5, 2}}); + std::unordered_map externalToInternalTypes( + {{ExtVdwType{1}, VdwType{0}}, + {ExtVdwType{2}, VdwType{1}}, + {ExtVdwType{5}, VdwType{2}}} + ); _nonCoulombPotential->determineInternalGlobalVdwTypes( externalToInternalTypes ); @@ -437,10 +443,10 @@ TEST_F( ); _nonCoulombPotential->fillOffDiagOfNonCoulPairsMatrix(); - EXPECT_EQ(getNonCoulombPairsMatrix()(0, 1)->getInternalType1(), 0); - EXPECT_EQ(getNonCoulombPairsMatrix()(0, 1)->getInternalType2(), 1); - EXPECT_EQ(getNonCoulombPairsMatrix()(1, 0)->getInternalType1(), 0); - EXPECT_EQ(getNonCoulombPairsMatrix()(1, 0)->getInternalType2(), 1); + EXPECT_EQ(getNonCoulombPairsMatrix()(0, 1)->getInternalType1(), VdwType{0}); + EXPECT_EQ(getNonCoulombPairsMatrix()(0, 1)->getInternalType2(), VdwType{1}); + EXPECT_EQ(getNonCoulombPairsMatrix()(1, 0)->getInternalType1(), VdwType{0}); + EXPECT_EQ(getNonCoulombPairsMatrix()(1, 0)->getInternalType2(), VdwType{1}); } /** @@ -455,35 +461,36 @@ TEST_F( { _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(2), - static_cast(1), + ExtVdwType(2), + ExtVdwType(1), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(5), - static_cast(1), + ExtVdwType(5), + ExtVdwType(1), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(5), - static_cast(2), + ExtVdwType(5), + ExtVdwType(2), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); // these two lines were already tested in // TestPotential_determineInternalGlobalVdwTypes - std::map externalToInternalTypes({{1, 0}, {2, 1}, {5, 2}}); + std::unordered_map externalToInternalTypes( + {{ExtVdwType{1}, VdwType{0}}, + {ExtVdwType{2}, VdwType{1}}, + {ExtVdwType{5}, VdwType{2}}} + ); _nonCoulombPotential->determineInternalGlobalVdwTypes( externalToInternalTypes ); @@ -492,10 +499,10 @@ TEST_F( ); _nonCoulombPotential->fillOffDiagOfNonCoulPairsMatrix(); - EXPECT_EQ(getNonCoulombPairsMatrix()(0, 1)->getInternalType1(), 1); - EXPECT_EQ(getNonCoulombPairsMatrix()(0, 1)->getInternalType2(), 0); - EXPECT_EQ(getNonCoulombPairsMatrix()(1, 0)->getInternalType1(), 1); - EXPECT_EQ(getNonCoulombPairsMatrix()(1, 0)->getInternalType2(), 0); + EXPECT_EQ(getNonCoulombPairsMatrix()(0, 1)->getInternalType1(), VdwType{1}); + EXPECT_EQ(getNonCoulombPairsMatrix()(0, 1)->getInternalType2(), VdwType{0}); + EXPECT_EQ(getNonCoulombPairsMatrix()(1, 0)->getInternalType1(), VdwType{1}); + EXPECT_EQ(getNonCoulombPairsMatrix()(1, 0)->getInternalType2(), VdwType{0}); } /** @@ -510,44 +517,44 @@ TEST_F( { _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(2), - static_cast(1), + ExtVdwType(2), + ExtVdwType(1), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(2), + ExtVdwType(1), + ExtVdwType(2), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(5), + ExtVdwType(1), + ExtVdwType(5), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(2), - static_cast(5), + ExtVdwType(2), + ExtVdwType(5), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); // these two lines were already tested in // TestPotential_determineInternalGlobalVdwTypes - std::map externalToInternalTypes({{1, 0}, {2, 1}, {5, 2}}); + std::unordered_map externalToInternalTypes( + {{ExtVdwType{1}, VdwType{0}}, + {ExtVdwType{2}, VdwType{1}}, + {ExtVdwType{5}, VdwType{2}}} + ); _nonCoulombPotential->determineInternalGlobalVdwTypes( externalToInternalTypes ); @@ -556,10 +563,10 @@ TEST_F( ); _nonCoulombPotential->fillOffDiagOfNonCoulPairsMatrix(); - EXPECT_EQ(getNonCoulombPairsMatrix()(0, 1)->getInternalType1(), 0); - EXPECT_EQ(getNonCoulombPairsMatrix()(0, 1)->getInternalType2(), 1); - EXPECT_EQ(getNonCoulombPairsMatrix()(1, 0)->getInternalType1(), 0); - EXPECT_EQ(getNonCoulombPairsMatrix()(1, 0)->getInternalType2(), 1); + EXPECT_EQ(getNonCoulombPairsMatrix()(0, 1)->getInternalType1(), VdwType{0}); + EXPECT_EQ(getNonCoulombPairsMatrix()(0, 1)->getInternalType2(), VdwType{1}); + EXPECT_EQ(getNonCoulombPairsMatrix()(1, 0)->getInternalType1(), VdwType{0}); + EXPECT_EQ(getNonCoulombPairsMatrix()(1, 0)->getInternalType2(), VdwType{1}); } /** @@ -574,26 +581,28 @@ TEST_F( { _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(2), + ExtVdwType(1), + ExtVdwType(2), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(2), - static_cast(1), + ExtVdwType(2), + ExtVdwType(1), 5.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); // these two lines were already tested in // TestPotential_determineInternalGlobalVdwTypes - std::map externalToInternalTypes({{1, 0}, {2, 1}, {5, 2}}); + std::unordered_map externalToInternalTypes( + {{ExtVdwType{1}, VdwType{0}}, + {ExtVdwType{2}, VdwType{1}}, + {ExtVdwType{5}, VdwType{2}}} + ); _nonCoulombPotential->determineInternalGlobalVdwTypes( externalToInternalTypes ); @@ -604,8 +613,15 @@ TEST_F( EXPECT_THROW_MSG( _nonCoulombPotential->fillOffDiagOfNonCoulPairsMatrix(), exc::ParameterFileException, - "Non-coulombic pairs with global van der Waals types 1, 2 and 2, 1 in " - "the parameter file have different parameters" + std::format( + "Non-coulombic pairs with global van der Waals types {}, {} and " + "{}, {} in " + "the parameter file have different parameters", + ExtVdwType{1}.toString(), + ExtVdwType{2}.toString(), + ExtVdwType{2}.toString(), + ExtVdwType{1}.toString() + ) ); } @@ -617,44 +633,44 @@ TEST_F(TestNonCoulombPotentialFF, getSelfInteractionNonCoulPairs) { _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(5), + ExtVdwType(1), + ExtVdwType(5), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(1), - static_cast(2), + ExtVdwType(1), + ExtVdwType(2), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(2), - static_cast(2), + ExtVdwType(2), + ExtVdwType(2), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); _nonCoulombPotential->addNonCoulombicPair( std::make_shared( - static_cast(5), - static_cast(5), + ExtVdwType(5), + ExtVdwType(5), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ) ); // these two lines were already tested in // TestPotential_determineInternalGlobalVdwTypes - std::map externalToInternalTypes({{1, 0}, {2, 1}, {5, 2}}); + std::unordered_map externalToInternalTypes( + {{ExtVdwType{1}, VdwType{0}}, + {ExtVdwType{2}, VdwType{1}}, + {ExtVdwType{5}, VdwType{2}}} + ); _nonCoulombPotential->determineInternalGlobalVdwTypes( externalToInternalTypes ); @@ -674,72 +690,72 @@ TEST_F(TestNonCoulombPotentialFF, sortNonCoulombicsPairs) auto vector = std::vector>(); auto pair1 = std::make_shared( - static_cast(1), - static_cast(1), + ExtVdwType(1), + ExtVdwType(1), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ); - pair1->setInternalType1(1); - pair1->setInternalType2(5); + pair1->setInternalType1(VdwType{1}); + pair1->setInternalType2(VdwType{5}); vector.push_back(pair1); auto pair2 = std::make_shared( - static_cast(2), - static_cast(2), + ExtVdwType(2), + ExtVdwType(2), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ); - pair2->setInternalType1(2); - pair2->setInternalType2(2); + pair2->setInternalType1(VdwType{2}); + pair2->setInternalType2(VdwType{2}); vector.push_back(pair2); auto pair3 = std::make_shared( - static_cast(2), - static_cast(3), + ExtVdwType(2), + ExtVdwType(3), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ); - pair3->setInternalType1(2); - pair3->setInternalType2(3); + pair3->setInternalType1(VdwType{2}); + pair3->setInternalType2(VdwType{3}); vector.push_back(pair3); auto pair4 = std::make_shared( - static_cast(1), - static_cast(4), + ExtVdwType(1), + ExtVdwType(4), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ); - pair4->setInternalType1(1); - pair4->setInternalType2(4); + pair4->setInternalType1(VdwType{1}); + pair4->setInternalType2(VdwType{4}); vector.push_back(pair4); _nonCoulombPotential->sortNonCoulombicsPairs(vector); - EXPECT_EQ(vector[0]->getInternalType1(), 1); - EXPECT_EQ(vector[0]->getInternalType2(), 4); - EXPECT_EQ(vector[1]->getInternalType1(), 1); - EXPECT_EQ(vector[1]->getInternalType2(), 5); - EXPECT_EQ(vector[2]->getInternalType1(), 2); - EXPECT_EQ(vector[2]->getInternalType2(), 2); - EXPECT_EQ(vector[3]->getInternalType1(), 2); - EXPECT_EQ(vector[3]->getInternalType2(), 3); + EXPECT_EQ(vector[0]->getInternalType1(), VdwType{1}); + EXPECT_EQ(vector[0]->getInternalType2(), VdwType{4}); + EXPECT_EQ(vector[1]->getInternalType1(), VdwType{1}); + EXPECT_EQ(vector[1]->getInternalType2(), VdwType{5}); + EXPECT_EQ(vector[2]->getInternalType1(), VdwType{2}); + EXPECT_EQ(vector[2]->getInternalType2(), VdwType{2}); + EXPECT_EQ(vector[3]->getInternalType1(), VdwType{2}); + EXPECT_EQ(vector[3]->getInternalType2(), VdwType{3}); auto pair5 = std::make_shared( - static_cast(1), - static_cast(1), + ExtVdwType(1), + ExtVdwType(1), 2.0, - 1.0, - 1.0 + LJParams{.c6 = 1.0, .c12 = 1.0} ); - pair5->setInternalType1(1); - pair5->setInternalType2(5); + pair5->setInternalType1(VdwType{1}); + pair5->setInternalType2(VdwType{5}); vector.push_back(pair5); EXPECT_THROW_MSG( _nonCoulombPotential->sortNonCoulombicsPairs(vector), exc::ParameterFileException, - "Non-coulombic pairs with global van der Waals types 1 and 1 in the " - "parameter file are defined twice" + std::format( + "Non-coulombic pairs with global van der Waals types {} and {} in " + "the " + "parameter file are defined twice", + ExtVdwType{1}.toString(), + ExtVdwType{1}.toString() + ) ); } diff --git a/tests/src/potential/nonCoulomb/testGuffPair.cpp b/tests/src/potential/nonCoulomb/testGuffPair.cpp index 706307c1b..a1fad447f 100644 --- a/tests/src/potential/nonCoulomb/testGuffPair.cpp +++ b/tests/src/potential/nonCoulomb/testGuffPair.cpp @@ -22,11 +22,10 @@ #include // for Test, CmpHelperFloatingPointEQ, InitGo... -#include // for ::pow, ::exp -#include // for vector, allocator +#include // for ::pow, ::exp -#include "gtest/gtest.h" // for Message, TestPartResult -#include "guffPair.hpp" // for GuffPair +#include "gtest/gtest.h" // for Message, TestPartResult +#include "guffPair.hpp" // for GuffPair using namespace potential; @@ -37,9 +36,10 @@ using namespace potential; TEST(TestGuffPair, calculateEnergyAndForces) { const auto coefficients = - std::vector{2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, - 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, - 18.0, 19.0, 20.0, 21.0, 22.0, 23.0}; + std::array{ + 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, + 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0 + }; const auto rncCutoff = 3.0; const double energyCutoff = 1.0; const double forceCutoff = 2.0; @@ -121,14 +121,16 @@ TEST(TestGuffPair, calculateWithSparseCoefficients) // and c19 blocks would otherwise compute pow(distance - rExp17/21, n) // with rExp17 == rExp21 == distance, which is 0^0 = NaN-prone; the gate // keeps the result finite. - auto coefficients = std::vector(22, 0.0); - coefficients[0] = 4.0; // c1 - coefficients[1] = 6.0; // n2 - coefficients[8] = 1.5; // c9 - coefficients[9] = 2.0; // cexp10 - coefficients[10] = 1.0; // rExp11 - coefficients[16] = distance; // rExp17 - would make distance - rExp17 == 0 - coefficients[17] = 3.0; // n18 + auto coefficients = std::array{}; + coefficients.fill(0.0); + coefficients[0] = 4.0; // c1 + coefficients[1] = 6.0; // n2 + coefficients[8] = 1.5; // c9 + coefficients[9] = 2.0; // cexp10 + coefficients[10] = 1.0; // rExp11 + coefficients[16] = + distance; // rExp17 - would make distance - rExp17 == 0 + coefficients[17] = 3.0; // n18 coefficients[20] = distance; // rExp21 coefficients[21] = 3.0; // n22 @@ -158,4 +160,4 @@ TEST(TestGuffPair, calculateWithSparseCoefficients) EXPECT_DOUBLE_EQ(force, expectedForce); EXPECT_FALSE(std::isnan(energy)); EXPECT_FALSE(std::isnan(force)); -} \ No newline at end of file +} diff --git a/tests/src/potential/nonCoulomb/testLennardJonesPair.cpp b/tests/src/potential/nonCoulomb/testLennardJonesPair.cpp index f733cce35..b091c2218 100644 --- a/tests/src/potential/nonCoulomb/testLennardJonesPair.cpp +++ b/tests/src/potential/nonCoulomb/testLennardJonesPair.cpp @@ -22,9 +22,8 @@ #include // for Test, CmpHelperFloatingPointEQ, EXPECT_EQ -#include // for pow -#include // for size_t -#include // for vector +#include // for pow +#include // for vector #include "gtest/gtest.h" // for AssertionResult, Message, TestPartResult #include "lennardJonesPair.hpp" // for LennardJonesPair @@ -32,39 +31,42 @@ using namespace potential; /** - * @brief tests the equals operator of BuckinghamPair + * @brief tests the equals operator of LennardJonesPair * */ TEST(TestLennardJonesPair, equalsOperator) { - const size_t vdwType1 = 0; - const size_t vdwType2 = 1; - const size_t vdwType3 = 2; - const auto nonCoulombPair1 = - LennardJonesPair(vdwType1, vdwType2, 1.0, 2.0, 3.0); + const ExtVdwType vdwType1{0}; + const ExtVdwType vdwType2{1}; + const ExtVdwType vdwType3{2}; + const auto ljParams1 = LJParams{.c6 = 2.0, .c12 = 3.0}; + const auto ljParams2 = LJParams{.c6 = 3.0, .c12 = 3.0}; + const auto ljParams3 = LJParams{.c6 = 2.0, .c12 = 4.0}; + const auto nonCoulombPair1 = + LennardJonesPair(vdwType1, vdwType2, 1.0, ljParams1); const auto nonCoulombPair2 = - LennardJonesPair(vdwType1, vdwType2, 1.0, 2.0, 3.0); + LennardJonesPair(vdwType1, vdwType2, 1.0, ljParams1); EXPECT_TRUE(nonCoulombPair1 == nonCoulombPair2); const auto nonCoulombPair3 = - LennardJonesPair(vdwType2, vdwType1, 1.0, 2.0, 3.0); + LennardJonesPair(vdwType2, vdwType1, 1.0, ljParams1); EXPECT_TRUE(nonCoulombPair1 == nonCoulombPair3); const auto nonCoulombPair4 = - LennardJonesPair(vdwType1, vdwType3, 1.0, 2.0, 3.0); + LennardJonesPair(vdwType1, vdwType3, 1.0, ljParams1); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair4); const auto nonCoulombPair5 = - LennardJonesPair(vdwType1, vdwType2, 2.0, 2.0, 3.0); + LennardJonesPair(vdwType1, vdwType2, 2.0, ljParams1); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair5); const auto nonCoulombPair6 = - LennardJonesPair(vdwType1, vdwType2, 1.0, 3.0, 3.0); + LennardJonesPair(vdwType1, vdwType2, 1.0, ljParams2); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair6); const auto nonCoulombPair7 = - LennardJonesPair(vdwType1, vdwType2, 1.0, 2.0, 4.0); + LennardJonesPair(vdwType1, vdwType2, 1.0, ljParams3); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair7); } @@ -83,8 +85,7 @@ TEST(TestLennardJonesPair, calculateEnergyAndForces) rncCutoff, energyCutoff, forceCutoff, - coefficients[0], - coefficients[1] + LJParams{.c6 = coefficients[0], .c12 = coefficients[1]} ); auto distance = 2.0; diff --git a/tests/src/potential/nonCoulomb/testMorsePair.cpp b/tests/src/potential/nonCoulomb/testMorsePair.cpp index fa317471e..e490eb06b 100644 --- a/tests/src/potential/nonCoulomb/testMorsePair.cpp +++ b/tests/src/potential/nonCoulomb/testMorsePair.cpp @@ -22,12 +22,12 @@ #include // for Test, CmpHelperFloatingPointEQ, EXPECT_EQ -#include // for exp, pow -#include // for size_t -#include // for vector +#include // for exp, pow +#include // for vector #include "gtest/gtest.h" // for AssertionResult, Message, TestPartResult #include "morsePair.hpp" // for MorsePair +#include "strongTypes.hpp" using namespace potential; @@ -37,38 +37,103 @@ using namespace potential; */ TEST(TestMorsePair, equalsOperator) { - const size_t vdwType1 = 0; - const size_t vdwType2 = 1; - const size_t vdwType3 = 2; - const auto nonCoulombPair1 = - MorsePair(vdwType1, vdwType2, 1.0, 2.0, 3.0, 4.0); - - const auto nonCoulombPair2 = - MorsePair(vdwType1, vdwType2, 1.0, 2.0, 3.0, 4.0); + const ExtVdwType vdwType1{0}; + const ExtVdwType vdwType2{1}; + const ExtVdwType vdwType3{2}; + + const auto nonCoulombPair1 = MorsePair( + vdwType1, + vdwType2, + 1.0, + MorseParams{ + .dissociationEnergy = 2.0, + .wellWidth = 3.0, + .equilibriumDistance = 4.0 + } + ); + + const auto nonCoulombPair2 = MorsePair( + vdwType1, + vdwType2, + 1.0, + MorseParams{ + .dissociationEnergy = 2.0, + .wellWidth = 3.0, + .equilibriumDistance = 4.0 + } + ); EXPECT_TRUE(nonCoulombPair1 == nonCoulombPair2); - const auto nonCoulombPair3 = - MorsePair(vdwType2, vdwType1, 1.0, 2.0, 3.0, 4.0); + const auto nonCoulombPair3 = MorsePair( + vdwType2, + vdwType1, + 1.0, + MorseParams{ + .dissociationEnergy = 2.0, + .wellWidth = 3.0, + .equilibriumDistance = 4.0 + } + ); EXPECT_TRUE(nonCoulombPair1 == nonCoulombPair3); - const auto nonCoulombPair4 = - MorsePair(vdwType1, vdwType3, 1.0, 2.0, 3.0, 4.0); + const auto nonCoulombPair4 = MorsePair( + vdwType1, + vdwType3, + 1.0, + MorseParams{ + .dissociationEnergy = 2.0, + .wellWidth = 3.0, + .equilibriumDistance = 4.0 + } + ); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair4); - const auto nonCoulombPair5 = - MorsePair(vdwType1, vdwType2, 2.0, 2.0, 3.0, 4.0); + const auto nonCoulombPair5 = MorsePair( + vdwType1, + vdwType2, + 2.0, + MorseParams{ + .dissociationEnergy = 2.0, + .wellWidth = 3.0, + .equilibriumDistance = 4.0 + } + ); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair5); - const auto nonCoulombPair6 = - MorsePair(vdwType1, vdwType2, 1.0, 3.0, 3.0, 4.0); + const auto nonCoulombPair6 = MorsePair( + vdwType1, + vdwType2, + 1.0, + MorseParams{ + .dissociationEnergy = 3.0, + .wellWidth = 3.0, + .equilibriumDistance = 4.0 + } + ); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair6); - const auto nonCoulombPair7 = - MorsePair(vdwType1, vdwType2, 1.0, 2.0, 4.0, 4.0); + const auto nonCoulombPair7 = MorsePair( + vdwType1, + vdwType2, + 1.0, + MorseParams{ + .dissociationEnergy = 2.0, + .wellWidth = 4.0, + .equilibriumDistance = 4.0 + } + ); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair7); - const auto nonCoulombPair8 = - MorsePair(vdwType1, vdwType2, 1.0, 2.0, 3.0, 5.0); + const auto nonCoulombPair8 = MorsePair( + vdwType1, + vdwType2, + 1.0, + MorseParams{ + .dissociationEnergy = 2.0, + .wellWidth = 3.0, + .equilibriumDistance = 5.0 + } + ); EXPECT_FALSE(nonCoulombPair1 == nonCoulombPair8); } @@ -87,9 +152,11 @@ TEST(TestMorsePair, calculateEnergyAndForces) rncCutoff, energyCutoff, forceCutoff, - coefficients[0], - coefficients[1], - coefficients[2] + MorseParams{ + .dissociationEnergy = coefficients[0], + .wellWidth = coefficients[1], + .equilibriumDistance = coefficients[2] + } ); const auto distance = 2.0; diff --git a/tests/src/potential/testBruteForceCellListEquivalence.cpp b/tests/src/potential/testBruteForceCellListEquivalence.cpp index 3edb2b3ca..75dcf6bcf 100644 --- a/tests/src/potential/testBruteForceCellListEquivalence.cpp +++ b/tests/src/potential/testBruteForceCellListEquivalence.cpp @@ -39,6 +39,7 @@ #include "potentialCellList.hpp" #include "potentialSettings.hpp" #include "simulationBox.hpp" +#include "strongTypes.hpp" using linearAlgebra::Vec3D; using molsys::Atom; @@ -111,7 +112,7 @@ namespace atom->setAtomType(0); atom->setExternalAtomType(p.molType); atom->setPartialCharge(p.molType == 1 ? 0.5 : -0.3); - atom->setInternalGlobalVDWType(0); + atom->setInternalGlobalVDWType(VdwType{0}); atom->setForceToZero(); Molecule molecule; @@ -141,8 +142,7 @@ namespace const auto pair = std::make_shared( kCoulombCutOff, - /*c6=*/-1.0, - /*c12=*/1.0 + LJParams{.c6 = -1.0, .c12 = 1.0} ); for (size_t m1 = 1; m1 <= 2; ++m1) diff --git a/tests/src/potential/testPairPotentialDerivatives.cpp b/tests/src/potential/testPairPotentialDerivatives.cpp index b3a602035..927d813f4 100644 --- a/tests/src/potential/testPairPotentialDerivatives.cpp +++ b/tests/src/potential/testPairPotentialDerivatives.cpp @@ -22,14 +22,13 @@ #include -#include - #include "buckinghamPair.hpp" #include "coulombShiftedPotential.hpp" #include "coulombWolf.hpp" #include "guffPair.hpp" #include "lennardJonesPair.hpp" #include "morsePair.hpp" +#include "strongTypes.hpp" namespace { @@ -52,7 +51,7 @@ namespace ) { const auto [energy, force] = calculate(distance); - (void)energy; + (void) energy; const auto energyDerivative = centralDifference( [&calculate](const double r) { return calculate(r).first; }, @@ -62,20 +61,22 @@ namespace EXPECT_NEAR(force, -energyDerivative, tolerance); } - std::vector buildGuffCoefficients() + std::array buildGuffCoefficients() { - return { - 0.5, 2.0, -0.2, 4.0, 0.1, 6.0, -0.05, 8.0, - 0.4, 1.2, 1.3, -0.3, -1.1, 2.5, 0.02, -0.5, - 0.3, 2.0, -0.01, -0.25, -0.4, 2.0 - }; + return {0.5, 2.0, -0.2, 4.0, 0.1, 6.0, -0.05, 8.0, + 0.4, 1.2, 1.3, -0.3, -1.1, 2.5, 0.02, -0.5, + 0.3, 2.0, -0.01, -0.25, -0.4, 2.0}; } } // namespace TEST(TestPairPotentialDerivatives, LennardJonesForceIsNegativeEnergyDerivative) { - const auto potential = - potential::LennardJonesPair(4.0, 0.15, -0.2, -1.0, 1.5); + const auto potential = potential::LennardJonesPair( + 4.0, + 0.15, + -0.2, + LJParams{.c6 = -1.0, .c12 = 1.5} + ); expectForceIsNegativeEnergyDerivative( [&potential](const double r) { return potential.calculate(r); }, @@ -86,8 +87,12 @@ TEST(TestPairPotentialDerivatives, LennardJonesForceIsNegativeEnergyDerivative) TEST(TestPairPotentialDerivatives, BuckinghamForceIsNegativeEnergyDerivative) { - const auto potential = - potential::BuckinghamPair(4.0, 0.25, -0.1, 2.0, -1.1, -0.4); + const auto potential = potential::BuckinghamPair( + 4.0, + 0.25, + -0.1, + BuckinghamParams{.scaling = 2.0, .dRho = -1.1, .c6 = -0.4} + ); expectForceIsNegativeEnergyDerivative( [&potential](const double r) { return potential.calculate(r); }, @@ -98,8 +103,16 @@ TEST(TestPairPotentialDerivatives, BuckinghamForceIsNegativeEnergyDerivative) TEST(TestPairPotentialDerivatives, MorseForceIsNegativeEnergyDerivative) { - const auto potential = - potential::MorsePair(4.0, 0.3, -0.2, 2.5, 1.4, 1.1); + const auto potential = potential::MorsePair( + 4.0, + 0.3, + -0.2, + MorseParams{ + .dissociationEnergy = 2.5, + .wellWidth = 1.4, + .equilibriumDistance = 1.1 + } + ); expectForceIsNegativeEnergyDerivative( [&potential](const double r) { return potential.calculate(r); }, @@ -120,7 +133,10 @@ TEST(TestPairPotentialDerivatives, GuffForceIsNegativeEnergyDerivative) ); } -TEST(TestPairPotentialDerivatives, ShiftedCoulombForceIsNegativeEnergyDerivative) +TEST( + TestPairPotentialDerivatives, + ShiftedCoulombForceIsNegativeEnergyDerivative +) { const auto potential = potential::CoulombShiftedPotential(4.0); const auto chargeProduct = 0.75; @@ -170,10 +186,19 @@ TEST(TestPairPotentialDerivatives, NonCoulombShiftedPairsAreZeroAtCutoff) constexpr auto cutoff = 4.0; const auto lennardJonesUnshifted = - potential::LennardJonesPair(cutoff, -1.0, 1.5); - const auto buckinghamUnshifted = - potential::BuckinghamPair(cutoff, 2.0, -1.1, -0.4); - const auto morseUnshifted = potential::MorsePair(cutoff, 2.5, 1.4, 1.1); + potential::LennardJonesPair(cutoff, LJParams{.c6 = -1.0, .c12 = 1.5}); + const auto buckinghamUnshifted = potential::BuckinghamPair( + cutoff, + BuckinghamParams{.scaling = 2.0, .dRho = -1.1, .c6 = -0.4} + ); + const auto morseUnshifted = potential::MorsePair( + cutoff, + MorseParams{ + .dissociationEnergy = 2.5, + .wellWidth = 1.4, + .equilibriumDistance = 1.1 + } + ); const auto guffUnshifted = potential::GuffPair(cutoff, buildGuffCoefficients()); @@ -190,24 +215,23 @@ TEST(TestPairPotentialDerivatives, NonCoulombShiftedPairsAreZeroAtCutoff) cutoff, ljEnergyCutoff, ljForceCutoff, - -1.0, - 1.5 + LJParams{.c6 = -1.0, .c12 = 1.5} ); const auto buckingham = potential::BuckinghamPair( cutoff, buckEnergyCutoff, buckForceCutoff, - 2.0, - -1.1, - -0.4 + BuckinghamParams{.scaling = 2.0, .dRho = -1.1, .c6 = -0.4} ); const auto morse = potential::MorsePair( cutoff, morseEnergyCutoff, morseForceCutoff, - 2.5, - 1.4, - 1.1 + MorseParams{ + .dissociationEnergy = 2.5, + .wellWidth = 1.4, + .equilibriumDistance = 1.1 + } ); const auto guff = potential::GuffPair( cutoff, diff --git a/tests/src/setup/testPotentialSetup.cpp b/tests/src/setup/testPotentialSetup.cpp index f614c66ef..d6069baac 100644 --- a/tests/src/setup/testPotentialSetup.cpp +++ b/tests/src/setup/testPotentialSetup.cpp @@ -22,8 +22,7 @@ #include // for TestInfo (ptr only), EXPECT_EQ -#include // for size_t -#include // for make_shared +#include // for make_shared #include "coulombReactionField.hpp" // for CoulombReactionField #include "coulombShiftedPotential.hpp" // for CoulombShiftedPotential @@ -37,7 +36,8 @@ #include "moleculeType.hpp" // for MoleculeType #include "potentialSettings.hpp" // for PotentialSettings #include "potentialSetup.hpp" // for PotentialSetup, setupPotential -#include "testSetup.hpp" // for TestSetup +#include "strongTypes.hpp" +#include "testSetup.hpp" // for TestSetup #include "testUtils.hpp" #include "throwWithMessage.hpp" // for EXPECT_THROW_MSG @@ -126,8 +126,8 @@ TEST_F(TestSetup, setupNonCoulombicPairs) PotentialSetup potentialSetup(*_engine); auto molecule = molsys::MoleculeType(1); - molecule.addExternalGlobalVDWType(0); - molecule.addExternalGlobalVDWType(1); + molecule.addExternalGlobalVDWType(ExtVdwType{0}); + molecule.addExternalGlobalVDWType(ExtVdwType{1}); _engine->getSimulationBox().addMoleculeType(molecule); @@ -142,13 +142,17 @@ TEST_F(TestSetup, setupNonCoulombicPairs) _engine->getPotential()->getNonCoulombPotential() ); - const auto zero = static_cast(0); - const auto one = static_cast(1); - - auto nonCoulombPair1 = LennardJonesPair(zero, zero, 10.0, 2.0, 3.0); - auto nonCoulombPair2 = LennardJonesPair(one, zero, 10.0, 2.0, 3.0); - auto nonCoulombPair3 = LennardJonesPair(zero, one, 10.0, 2.0, 3.0); - auto nonCoulombPair4 = LennardJonesPair(one, one, 10.0, 2.0, 3.0); + const auto zero = ExtVdwType(0); + const auto one = ExtVdwType(1); + + auto nonCoulombPair1 = + LennardJonesPair(zero, zero, 10.0, LJParams{.c6 = 2.0, .c12 = 3.0}); + auto nonCoulombPair2 = + LennardJonesPair(one, zero, 10.0, LJParams{.c6 = 2.0, .c12 = 3.0}); + auto nonCoulombPair3 = + LennardJonesPair(zero, one, 10.0, LJParams{.c6 = 2.0, .c12 = 3.0}); + auto nonCoulombPair4 = + LennardJonesPair(one, one, 10.0, LJParams{.c6 = 2.0, .c12 = 3.0}); nonCoulombPotential.addNonCoulombicPair( std::make_shared(nonCoulombPair1) diff --git a/tests/src/setup/testSimulationBoxSetup.cpp b/tests/src/setup/testSimulationBoxSetup.cpp index 808885beb..cd9715189 100644 --- a/tests/src/setup/testSimulationBoxSetup.cpp +++ b/tests/src/setup/testSimulationBoxSetup.cpp @@ -40,7 +40,8 @@ #include "simulationBox.hpp" // for SimulationBox #include "simulationBoxSettings.hpp" // for SimulationBoxSettings #include "simulationBoxSetup.hpp" // for SimulationBoxSetup, setupSimulationBox -#include "testSetup.hpp" // for TestSetup +#include "strongTypes.hpp" +#include "testSetup.hpp" // for TestSetup using setup::molsys::SimulationBoxSetup; @@ -162,9 +163,9 @@ TEST_F(TestSetup, setExternalVDWTypes) ::molsys::MoleculeType moleculeType(1); moleculeType.setNumberOfAtoms(3); - moleculeType.addExternalGlobalVDWType(0); - moleculeType.addExternalGlobalVDWType(1); - moleculeType.addExternalGlobalVDWType(2); + moleculeType.addExternalGlobalVDWType(ExtVdwType{0}); + moleculeType.addExternalGlobalVDWType(ExtVdwType{1}); + moleculeType.addExternalGlobalVDWType(ExtVdwType{2}); _engine->getSimulationBox().addMolecule(molecule); _engine->getSimulationBox().addMolecule(qmMolecule); @@ -177,21 +178,11 @@ TEST_F(TestSetup, setExternalVDWTypes) SimulationBoxSetup simulationBoxSetup(*_engine); simulationBoxSetup.setExternalVDWTypes(); - EXPECT_EQ( - _engine->getSimulationBox().getMolecules()[0].getExternalGlobalVDWTypes( - )[0], - 0 - ); - EXPECT_EQ( - _engine->getSimulationBox().getMolecules()[0].getExternalGlobalVDWTypes( - )[1], - 1 - ); - EXPECT_EQ( - _engine->getSimulationBox().getMolecules()[0].getExternalGlobalVDWTypes( - )[2], - 2 - ); + const auto &moleculeResult = _engine->getSimulationBox().getMolecules()[0]; + + EXPECT_EQ(moleculeResult.getExternalGlobalVDWTypes()[0], ExtVdwType{0}); + EXPECT_EQ(moleculeResult.getExternalGlobalVDWTypes()[1], ExtVdwType{1}); + EXPECT_EQ(moleculeResult.getExternalGlobalVDWTypes()[2], ExtVdwType{2}); } TEST_F(TestSetup, setPartialCharges) @@ -224,18 +215,11 @@ TEST_F(TestSetup, setPartialCharges) SimulationBoxSetup simulationBoxSetup(*_engine); simulationBoxSetup.setPartialCharges(); - EXPECT_EQ( - _engine->getSimulationBox().getMolecules()[0].getPartialCharges()[0], - 0.0 - ); - EXPECT_EQ( - _engine->getSimulationBox().getMolecules()[0].getPartialCharges()[1], - 1.0 - ); - EXPECT_EQ( - _engine->getSimulationBox().getMolecules()[0].getPartialCharges()[2], - 2.0 - ); + const auto &moleculeResult = _engine->getSimulationBox().getMolecules()[0]; + + EXPECT_EQ(moleculeResult.getPartialCharges()[0], 0.0); + EXPECT_EQ(moleculeResult.getPartialCharges()[1], 1.0); + EXPECT_EQ(moleculeResult.getPartialCharges()[2], 2.0); } TEST_F(TestSetup, testSetAtomMasses) @@ -256,18 +240,11 @@ TEST_F(TestSetup, testSetAtomMasses) SimulationBoxSetup simulationBoxSetup(*_engine); simulationBoxSetup.setAtomMasses(); - EXPECT_DOUBLE_EQ( - _engine->getSimulationBox().getMolecules()[0].getAtomMass(0), - 12.0107 - ); - EXPECT_DOUBLE_EQ( - _engine->getSimulationBox().getMolecules()[0].getAtomMass(1), - 1.00794 - ); - EXPECT_DOUBLE_EQ( - _engine->getSimulationBox().getMolecules()[0].getAtomMass(2), - 15.9994 - ); + const auto &moleculeResult = _engine->getSimulationBox().getMolecules()[0]; + + EXPECT_DOUBLE_EQ(moleculeResult.getAtomMass(0), 12.0107); + EXPECT_DOUBLE_EQ(moleculeResult.getAtomMass(1), 1.00794); + EXPECT_DOUBLE_EQ(moleculeResult.getAtomMass(2), 15.9994); } TEST_F(TestSetup, testSetAtomMassesThrowsError) @@ -310,18 +287,11 @@ TEST_F(TestSetup, testSetAtomicNumbers) SimulationBoxSetup simulationBoxSetup(*_engine); simulationBoxSetup.setAtomicNumbers(); - EXPECT_EQ( - _engine->getSimulationBox().getMolecules()[0].getAtomicNumber(0), - AtomNumber{6} - ); - EXPECT_EQ( - _engine->getSimulationBox().getMolecules()[0].getAtomicNumber(1), - AtomNumber{1} - ); - EXPECT_EQ( - _engine->getSimulationBox().getMolecules()[0].getAtomicNumber(2), - AtomNumber{8} - ); + const auto &moleculeResult = _engine->getSimulationBox().getMolecules()[0]; + + EXPECT_EQ(moleculeResult.getAtomicNumber(0), AtomNumber{6}); + EXPECT_EQ(moleculeResult.getAtomicNumber(1), AtomNumber{1}); + EXPECT_EQ(moleculeResult.getAtomicNumber(2), AtomNumber{8}); } TEST_F(TestSetup, testSetAtomicNumbersThrowsError) @@ -561,9 +531,9 @@ TEST_F(TestSetup, testFullSetup) moleculeType1.addExternalAtomType(0); moleculeType1.addExternalAtomType(0); moleculeType1.addExternalAtomType(0); - moleculeType1.addExternalGlobalVDWType(0); - moleculeType1.addExternalGlobalVDWType(1); - moleculeType1.addExternalGlobalVDWType(2); + moleculeType1.addExternalGlobalVDWType(ExtVdwType{0}); + moleculeType1.addExternalGlobalVDWType(ExtVdwType{1}); + moleculeType1.addExternalGlobalVDWType(ExtVdwType{2}); moleculeType2.setNumberOfAtoms(2); moleculeType2.addAtomName("H"); @@ -574,8 +544,8 @@ TEST_F(TestSetup, testFullSetup) moleculeType2.addAtomType(0); moleculeType2.addExternalAtomType(0); moleculeType2.addExternalAtomType(0); - moleculeType2.addExternalGlobalVDWType(0); - moleculeType2.addExternalGlobalVDWType(1); + moleculeType2.addExternalGlobalVDWType(ExtVdwType{0}); + moleculeType2.addExternalGlobalVDWType(ExtVdwType{1}); _engine->getSimulationBox().addMoleculeType(moleculeType1); _engine->getSimulationBox().addMoleculeType(moleculeType2); diff --git a/tests/src/testUtils/CMakeLists.txt b/tests/src/testUtils/CMakeLists.txt index 4b8f1f630..17d85ada1 100644 --- a/tests/src/testUtils/CMakeLists.txt +++ b/tests/src/testUtils/CMakeLists.txt @@ -1,5 +1,6 @@ add_library(testUtils testUtils.cpp + testNonCoulombPairUtils.cpp ) target_include_directories(testUtils diff --git a/tests/src/testUtils/testNonCoulombPairUtils.cpp b/tests/src/testUtils/testNonCoulombPairUtils.cpp new file mode 100644 index 000000000..0d1071667 --- /dev/null +++ b/tests/src/testUtils/testNonCoulombPairUtils.cpp @@ -0,0 +1,78 @@ +/***************************************************************************** + + + PQ + Copyright (C) 2023-now Jakob Gamper + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + + +******************************************************************************/ + +#include "testNonCoulombPairUtils.hpp" + +#include "lennardJonesPair.hpp" +#include "morsePair.hpp" + +/** + * @brief Get the MorseParams from a MorsePair. + * + * @param morsePair pointer to the MorsePair object + * @return const potential::MorseParams& reference to the MorseParams + */ +const MorseParams& TestMorsePairUtils::params( + const potential::MorsePair* morsePair +) +{ + return morsePair->_params; +} + +/** + * @brief Get the LJParams from a LennardJonesPair. + * + * @param ljPair pointer to the LennardJonesPair object + * @return const LJParams& reference to the LJParams + */ +const LJParams& TestLJPairUtils::params( + const potential::LennardJonesPair* ljPair +) +{ + return ljPair->_params; +} + +/** + * @brief Get the BuckinghamParams from a BuckinghamPair. + * + * @param buckPair pointer to the BuckinghamPair object + * @return const BuckinghamParams& reference to the BuckinghamParams + */ +const BuckinghamParams& TestBuckinghamPairUtils::params( + const potential::BuckinghamPair* buckPair +) +{ + return buckPair->_params; +} + +/** + * @brief Get the coefficients from a GuffPair. + * + * @param guffPair pointer to the GuffPair object + * @return const std::array& reference + * to the coefficients + */ +const std::array& TestGuffPairUtils:: + coeffs(const potential::GuffPair* guffPair) +{ + return guffPair->_coefficients; +} diff --git a/tests/src/waterModels/testWaterModels.cpp b/tests/src/waterModels/testWaterModels.cpp index b354ce8d8..162efeff0 100644 --- a/tests/src/waterModels/testWaterModels.cpp +++ b/tests/src/waterModels/testWaterModels.cpp @@ -104,7 +104,7 @@ namespace oxygen->setQMCharge(-0.9); oxygen->setPosition(origin); oxygen->setAtomType(0); - oxygen->setInternalGlobalVDWType(0); + oxygen->setInternalGlobalVDWType(VdwType{0}); oxygen->setForceToZero(); h1->setAtomicNumber(AtomNumber{1}); @@ -112,7 +112,7 @@ namespace h1->setQMCharge(0.45); h1->setPosition(origin + Vec3D{geometry.oh1, 0.0, 0.0}); h1->setAtomType(1); - h1->setInternalGlobalVDWType(0); + h1->setInternalGlobalVDWType(VdwType{0}); h1->setForceToZero(); h2->setAtomicNumber(AtomNumber{1}); @@ -127,7 +127,7 @@ namespace } ); h2->setAtomType(1); - h2->setInternalGlobalVDWType(0); + h2->setInternalGlobalVDWType(VdwType{0}); h2->setForceToZero(); Molecule water; @@ -219,8 +219,7 @@ namespace const auto pair = std::make_shared( kCutOff, - /*c6=*/-1.0, - /*c12=*/1.0 + LJParams{.c6 = -1.0, .c12 = 1.0} ); for (size_t mol1 = 1; mol1 <= 2; ++mol1) @@ -443,7 +442,10 @@ TEST(InterWater, PairEvaluatorsApplySymmetricAndOneWayForces) atom2.setForceToZero(); const auto coulomb = std::make_shared(kCutOff); - const LennardJonesPair nonCoulomb(kCutOff, -1.0, 1.0); + const LennardJonesPair nonCoulomb( + kCutOff, + LJParams{.c6 = -1.0, .c12 = 1.0} + ); ExposedInterWaterStrategy strategy; EXPECT_DOUBLE_EQ(nonCoulomb.getRadialCutOff(), kCutOff); @@ -525,11 +527,18 @@ TEST(InterWater, NonOxygenOnlyStateInitializesEveryPair) PotentialSettings::setCoulombRadiusCutOff(kCutOff); PotentialSettings::setNonCoulombRadiusCutOff(kCutOff); - auto oxygenOxygen = std::make_unique(kCutOff, -1.0, 1.0); - auto oxygenHydrogen = - std::make_unique(kCutOff, -1.0, 1.0); - auto hydrogenHydrogen = - std::make_unique(kCutOff, -1.0, 1.0); + auto oxygenOxygen = std::make_unique( + kCutOff, + LJParams{.c6 = -1.0, .c12 = 1.0} + ); + auto oxygenHydrogen = std::make_unique( + kCutOff, + LJParams{.c6 = -1.0, .c12 = 1.0} + ); + auto hydrogenHydrogen = std::make_unique( + kCutOff, + LJParams{.c6 = -1.0, .c12 = 1.0} + ); const auto *oxygenOxygenView = oxygenOxygen.get(); const auto *oxygenHydrogenView = oxygenHydrogen.get(); const auto *hydrogenHydrogenView = hydrogenHydrogen.get(); @@ -652,14 +661,14 @@ TEST(PotentialTemplates, QmChargesAndOneWayInteractions) atom1.setPartialCharge(-0.8); atom1.setQMCharge(-0.9); atom1.setAtomType(0); - atom1.setInternalGlobalVDWType(0); + atom1.setInternalGlobalVDWType(VdwType{0}); atom1.setForceToZero(); Atom atom2; atom2.setPosition({1.2, 0.1, 0.0}); atom2.setPartialCharge(0.4); atom2.setAtomType(0); - atom2.setInternalGlobalVDWType(0); + atom2.setInternalGlobalVDWType(VdwType{0}); atom2.setForceToZero(); HybridSettings::setUseQMCharges(true);