diff --git a/source/api_c/include/c_api.h b/source/api_c/include/c_api.h index a97ab6b79a..509895b8ce 100644 --- a/source/api_c/include/c_api.h +++ b/source/api_c/include/c_api.h @@ -34,6 +34,17 @@ typedef struct DP_DeepPot DP_DeepPot; **/ extern DP_DeepPot* DP_NewDeepPot(const char* c_model); +/** + * @brief DP constructor with initialization. + * + * @param c_model The name of the frozen model file. + * @param gpu_rank The rank of the GPU. + * @param c_file_content The content of the model file. + * @return DP_DeepPot* A pointer to the deep potential. + */ +extern DP_DeepPot* DP_NewDeepPotWithParam( + const char* c_model, const int gpu_rank, const char* c_file_content); + /** * @brief Evaluate the energy, force and virial by using a DP. (double version) * @param[in] dp The DP to use. @@ -279,6 +290,17 @@ typedef struct DP_DeepTensor DP_DeepTensor; **/ extern DP_DeepTensor* DP_NewDeepTensor(const char* c_model); +/** + * @brief Deep Tensor constructor with initialization. + * + * @param c_model The name of the frozen model file. + * @param gpu_rank The rank of the GPU. + * @param c_name_scope The name scope. + * @return DP_DeepTensor* A pointer to the deep tensor. + */ +extern DP_DeepTensor* DP_NewDeepTensorWithParam( + const char* c_model, const int gpu_rank, const char* c_name_scope); + /** * @brief Evaluate the tensor by using a DP. (double version) * @param[in] dt The Deep Tensor to use. @@ -525,6 +547,120 @@ int* DP_DeepTensorGetSelTypes(DP_DeepTensor* dt); */ int DP_DeepTensorGetNumbSelTypes(DP_DeepTensor* dt); +/** +* @brief The dipole charge modifier. +**/ +typedef struct DP_DipoleChargeModifier DP_DipoleChargeModifier; + +/** +* @brief Dipole charge modifier constructor with initialization. +* @param[in] c_model The name of the frozen model file. +* @returns A pointer to the dipole charge modifier. +**/ +extern DP_DipoleChargeModifier* DP_NewDipoleChargeModifier(const char* c_model); + +/** + * @brief Dipole charge modifier constructor with initialization. + * + * @param c_model The name of the frozen model file. + * @param gpu_rank The rank of the GPU. + * @param c_name_scope The name scope. + * @return DP_DipoleChargeModifier* A pointer to the dipole charge modifier. + */ +extern DP_DipoleChargeModifier* DP_NewDipoleChargeModifierWithParam( + const char* c_model, const int gpu_rank, const char* c_name_scope); + +/** +* @brief Evaluate the force and virial correction by using a dipole charge modifier with the neighbor list. (double version) +* @param[in] dcm The dipole charge modifier to use. +* @param[in] natoms The number of atoms. +* @param[in] coord The coordinates of atoms. The array should be of size natoms x 3. +* @param[in] atype The atom types. The array should contain natoms ints. +* @param[in] cell The cell of the region. The array should be of size 9. Pass NULL if pbc is not used. +* @param[in] pairs The pairs of atoms. The list should contain npairs pairs of ints. +* @param[in] npairs The number of pairs. +* @param[in] delef_ The electric field on each atom. The array should be of size nframes x natoms x 3. +* @param[in] nghost The number of ghost atoms. +* @param[in] nlist The neighbor list. +* @param[out] dfcorr_ Output force correction. The array should be of size natoms x 3. +* @param[out] dvcorr_ Output virial correction. The array should be of size 9. +* @warning The output arrays should be allocated before calling this function. Pass NULL if not required. + **/ +extern void DP_DipoleChargeModifierComputeNList ( + DP_DipoleChargeModifier* dcm, + const int natom, + const double* coord, + const int* atype, + const double* cell, + const int* pairs, + const int npairs, + const double* delef_, + const int nghost, + const DP_Nlist* nlist, + double* dfcorr_, + double* dvcorr_ + ); + + +/** +* @brief Evaluate the force and virial correction by using a dipole charge modifier with the neighbor list. (float version) +* @param[in] dcm The dipole charge modifier to use. +* @param[in] natoms The number of atoms. +* @param[in] coord The coordinates of atoms. The array should be of size natoms x 3. +* @param[in] atype The atom types. The array should contain natoms ints. +* @param[in] cell The cell of the region. The array should be of size 9. Pass NULL if pbc is not used. +* @param[in] pairs The pairs of atoms. The list should contain npairs pairs of ints. +* @param[in] npairs The number of pairs. +* @param[in] delef_ The electric field on each atom. The array should be of size nframes x natoms x 3. +* @param[in] nghost The number of ghost atoms. +* @param[in] nlist The neighbor list. +* @param[out] dfcorr_ Output force correction. The array should be of size natoms x 3. +* @param[out] dvcorr_ Output virial correction. The array should be of size 9. +* @warning The output arrays should be allocated before calling this function. Pass NULL if not required. + **/ +extern void DP_DipoleChargeModifierComputeNListf ( + DP_DipoleChargeModifier* dcm, + const int natom, + const float* coord, + const int* atype, + const float* cell, + const int* pairs, + const int npairs, + const float* delef_, + const int nghost, + const DP_Nlist* nlist, + float* dfcorr_, + float* dvcorr_ + ); + +/** + * @brief Get the type map of a DipoleChargeModifier. + * @param[in] dcm The DipoleChargeModifier to use. + * @return The cutoff radius. +*/ +double DP_DipoleChargeModifierGetCutoff(DP_DipoleChargeModifier* dt); + +/** + * @brief Get the type map of a DipoleChargeModifier. + * @param[in] dcm The DipoleChargeModifier to use. + * @return The number of types of the DipoleChargeModifier. +*/ +int DP_DipoleChargeModifierGetNumbTypes(DP_DipoleChargeModifier* dt); + +/** + * @brief Get sel types of a DipoleChargeModifier. + * @param[in] dcm The DipoleChargeModifier to use. + * @return The sel types +*/ +int* DP_DipoleChargeModifierGetSelTypes(DP_DipoleChargeModifier* dt); + +/** + * @brief Get the number of sel types of a DipoleChargeModifier. + * @param[in] dcm The DipoleChargeModifier to use. + * @return The number of sel types +*/ +int DP_DipoleChargeModifierGetNumbSelTypes(DP_DipoleChargeModifier* dt); + /** * @brief Convert PBtxt to PB. * @param[in] c_pbtxt The name of the PBtxt file. diff --git a/source/api_c/include/c_api_internal.h b/source/api_c/include/c_api_internal.h index ec4eba1b17..43a4081f38 100644 --- a/source/api_c/include/c_api_internal.h +++ b/source/api_c/include/c_api_internal.h @@ -1,6 +1,7 @@ #include "neighbor_list.h" #include "DeepPot.h" #include "DeepTensor.h" +#include "DataModifier.h" struct DP_Nlist { DP_Nlist(deepmd::InputNlist& nl); @@ -25,3 +26,9 @@ struct DP_DeepTensor { deepmd::DeepTensor dt; }; + +struct DP_DipoleChargeModifier { + DP_DipoleChargeModifier(deepmd::DipoleChargeModifier& dcm); + + deepmd::DipoleChargeModifier dcm; +}; diff --git a/source/api_c/include/deepmd.hpp b/source/api_c/include/deepmd.hpp index 1007f92ef2..02ef27d331 100644 --- a/source/api_c/include/deepmd.hpp +++ b/source/api_c/include/deepmd.hpp @@ -347,6 +347,60 @@ inline void _DP_DeepTensorComputeNList( DP_DeepTensorComputeNListf(dt, natom, coord, atype, cell, nghost, nlist, global_tensor, force, virial, atomic_tensor, atomic_virial, size_at); } +template +inline void _DP_DipoleChargeModifierComputeNList( + DP_DipoleChargeModifier* dcm, + const int natom, + const FPTYPE* coord, + const int* atype, + const FPTYPE* cell, + const int* pairs, + const int npairs, + const FPTYPE* delef_, + const int nghost, + const DP_Nlist* nlist, + FPTYPE* dfcorr_, + FPTYPE* dvcorr_ +); + +template <> +inline void _DP_DipoleChargeModifierComputeNList( + DP_DipoleChargeModifier* dcm, + const int natom, + const double* coord, + const int* atype, + const double* cell, + const int* pairs, + const int npairs, + const double* delef_, + const int nghost, + const DP_Nlist* nlist, + double* dfcorr_, + double* dvcorr_ +) +{ + DP_DipoleChargeModifierComputeNList(dcm, natom, coord, atype, cell, pairs, npairs, delef_, nghost, nlist, dfcorr_, dvcorr_); +} + +template <> +inline void _DP_DipoleChargeModifierComputeNList( + DP_DipoleChargeModifier* dcm, + const int natom, + const float* coord, + const int* atype, + const float* cell, + const int* pairs, + const int npairs, + const float* delef_, + const int nghost, + const DP_Nlist* nlist, + float* dfcorr_, + float* dvcorr_ +) +{ + DP_DipoleChargeModifierComputeNListf(dcm, natom, coord, atype, cell, pairs, npairs, delef_, nghost, nlist, dfcorr_, dvcorr_); +} + namespace deepmd { namespace hpp @@ -431,23 +485,27 @@ namespace deepmd /** * @brief DP constructor with initialization. * @param[in] model The name of the frozen model file. + * @param[in] gpu_rank The GPU rank. + * @param[in] file_content The content of the frozen model file. **/ - DeepPot(const std::string &model) : dp(nullptr) + DeepPot(const std::string &model, const int &gpu_rank = 0, const std::string &file_content = "") : dp(nullptr) { - init(model); + init(model, gpu_rank, file_content); }; /** * @brief Initialize the DP. * @param[in] model The name of the frozen model file. + * @param[in] gpu_rank The GPU rank. + * @param[in] file_content The content of the frozen model file. **/ - void init(const std::string &model) + void init(const std::string &model, const int &gpu_rank = 0, const std::string &file_content = "") { if (dp) { std::cerr << "WARNING: deepmd-kit should not be initialized twice, do nothing at the second call of initializer" << std::endl; return; } - dp = DP_NewDeepPot(model.c_str()); + dp = DP_NewDeepPotWithParam(model.c_str(), gpu_rank, file_content.c_str()); }; /** @@ -849,22 +907,22 @@ namespace deepmd * @brief DeepTensor constructor with initialization. * @param[in] model The name of the frozen model file. **/ - DeepTensor(const std::string &model) : dt(nullptr) + DeepTensor(const std::string &model, const int &gpu_rank = 0, const std::string &name_scope = "") : dt(nullptr) { - init(model); + init(model, gpu_rank, name_scope); }; /** * @brief Initialize the DeepTensor. * @param[in] model The name of the frozen model file. **/ - void init(const std::string &model) + void init(const std::string &model, const int &gpu_rank = 0, const std::string &name_scope = "") { if (dt) { std::cerr << "WARNING: deepmd-kit should not be initialized twice, do nothing at the second call of initializer" << std::endl; return; } - dt = DP_NewDeepTensor(model.c_str()); + dt = DP_NewDeepTensorWithParam(model.c_str(), gpu_rank, name_scope.c_str()); odim = output_dim(); nsel_types = DP_DeepTensorGetNumbSelTypes(dt); }; @@ -1161,5 +1219,112 @@ namespace deepmd int odim; int nsel_types; }; + + class DipoleChargeModifier + { + public: + /** + * @brief DipoleChargeModifier constructor without initialization. + **/ + DipoleChargeModifier() : dcm(nullptr) {}; + ~DipoleChargeModifier(){}; + /** + * @brief DipoleChargeModifier constructor with initialization. + * @param[in] model The name of the frozen model file. + * @param[in] gpu_rank The rank of the GPU to be used. + * @param[in] name_scope The name scope of the model. + **/ + DipoleChargeModifier(const std::string &model, const int &gpu_rank = 0, const std::string &name_scope = "") : dcm(nullptr) + { + init(model, gpu_rank, name_scope); + }; + /** + * @brief Initialize the DipoleChargeModifier. + * @param[in] model The name of the frozen model file. + * @param[in] gpu_rank The rank of the GPU to be used. + * @param[in] name_scope The name scope of the model. + **/ + void init(const std::string &model, const int &gpu_rank = 0, const std::string &name_scope = "") + { + if (dcm) + { + std::cerr << "WARNING: deepmd-kit should not be initialized twice, do nothing at the second call of initializer" << std::endl; + return; + } + dcm = DP_NewDipoleChargeModifierWithParam(model.c_str(), gpu_rank, name_scope.c_str()); + nsel_types = DP_DipoleChargeModifierGetNumbSelTypes(dcm); + }; + /** + * @brief Evaluate the force and virial correction by using this dipole charge modifier. + * @param[out] dfcorr_ The force correction on each atom. + * @param[out] dvcorr_ The virial correction. + * @param[in] dcoord_ The coordinates of atoms. The array should be of size natoms x 3. + * @param[in] datype_ The atom types. The list should contain natoms ints. + * @param[in] dbox The cell of the region. The array should be of size 9. + * @param[in] pairs The pairs of atoms. The list should contain npairs pairs of ints. + * @param[in] delef_ The electric field on each atom. The array should be of size natoms x 3. + * @param[in] nghost The number of ghost atoms. + * @param[in] lmp_list The neighbor list. + **/ + template + void compute (std::vector & dfcorr_, + std::vector & dvcorr_, + const std::vector & dcoord_, + const std::vector & datype_, + const std::vector & dbox, + const std::vector> & pairs, + const std::vector & delef_, + const int nghost, + const InputNlist & lmp_list) + { + unsigned int natoms = datype_.size(); + assert(natoms * 3 == dcoord_.size()); + if (!dbox.empty()) + { + assert(dbox.size() == 9); + } + const VALUETYPE *dcoord = &dcoord_[0]; + const VALUETYPE *dbox_ = !dbox.empty() ? &dbox[0] : nullptr; + const int *datype = &datype_[0]; + const int npairs = pairs.size(); + const int *dpairs = reinterpret_cast(&pairs[0]); + const VALUETYPE *delef = &delef_[0]; + + dfcorr_.resize(natoms * 3); + dvcorr_.resize(9); + VALUETYPE *dfcorr = &dfcorr_[0]; + VALUETYPE *dvcorr = &dvcorr_[0]; + + _DP_DipoleChargeModifierComputeNList(dcm, natoms, dcoord, datype, dbox_, dpairs, npairs, delef, nghost, lmp_list.nl, dfcorr, dvcorr); + }; + /** + * @brief Get the cutoff radius. + * @return The cutoff radius. + **/ + double cutoff() const + { + assert(dcm); + return DP_DipoleChargeModifierGetCutoff(dcm); + }; + /** + * @brief Get the number of types. + * @return The number of types. + **/ + int numb_types() const + { + assert(dcm); + return DP_DipoleChargeModifierGetNumbTypes(dcm); + }; + + std::vector sel_types() const + { + int* sel_types_arr = DP_DipoleChargeModifierGetSelTypes(dcm); + std::vector sel_types_vec = std::vector(sel_types_arr, sel_types_arr + nsel_types); + return sel_types_vec; + } + private: + DP_DipoleChargeModifier *dcm; + int nsel_types; + }; } } diff --git a/source/api_c/src/c_api.cc b/source/api_c/src/c_api.cc index 5fbf40bcd2..396b9629a4 100644 --- a/source/api_c/src/c_api.cc +++ b/source/api_c/src/c_api.cc @@ -7,6 +7,7 @@ #include "common.h" #include "DeepPot.h" #include "DeepTensor.h" +#include "DataModifier.h" extern "C" { @@ -32,6 +33,15 @@ DP_DeepPot* DP_NewDeepPot(const char* c_model) { return new_dp; } +DP_DeepPot* DP_NewDeepPotWithParam( + const char* c_model, const int gpu_rank, const char* c_file_content) { + std::string model(c_model); + std::string file_content(c_file_content); + deepmd::DeepPot dp(model, gpu_rank, file_content); + DP_DeepPot* new_dp = new DP_DeepPot(dp); + return new_dp; +} + DP_DeepPotModelDevi::DP_DeepPotModelDevi(deepmd::DeepPotModelDevi& dp) : dp(dp) {} @@ -52,6 +62,34 @@ DP_DeepTensor* DP_NewDeepTensor(const char* c_model) { return new_dt; } +DP_DeepTensor* DP_NewDeepTensorWithParam( + const char* c_model, const int gpu_rank, const char* c_name_scope) { + std::string model(c_model); + std::string name_scope(c_name_scope); + deepmd::DeepTensor dt(model, gpu_rank, name_scope); + DP_DeepTensor* new_dt = new DP_DeepTensor(dt); + return new_dt; +} + +DP_DipoleChargeModifier::DP_DipoleChargeModifier(deepmd::DipoleChargeModifier& dcm) + : dcm(dcm) {} + +DP_DipoleChargeModifier* DP_NewDipoleChargeModifier(const char* c_model) { + std::string model(c_model); + deepmd::DipoleChargeModifier dcm(model); + DP_DipoleChargeModifier* new_dcm = new DP_DipoleChargeModifier(dcm); + return new_dcm; +} + +DP_DipoleChargeModifier* DP_NewDipoleChargeModifierWithParam( + const char* c_model, const int gpu_rank, const char* c_name_scope) { + std::string model(c_model); + std::string name_scope(c_name_scope); + deepmd::DipoleChargeModifier dcm(model, gpu_rank, name_scope); + DP_DipoleChargeModifier* new_dcm = new DP_DipoleChargeModifier(dcm); + return new_dcm; +} + } // extern "C" template @@ -536,6 +574,77 @@ void DP_DeepTensorComputeNList_variant ( int* size_at ); +template +inline +void DP_DipoleChargeModifierComputeNList_variant ( + DP_DipoleChargeModifier* dcm, + const int natoms, + const VALUETYPE* coord, + const int* atype, + const VALUETYPE* cell, + const int* pairs, + const int npairs, + const VALUETYPE* delef, + const int nghost, + const DP_Nlist* nlist, + VALUETYPE* dfcorr_, + VALUETYPE* dvcorr_ + ){ + // init C++ vectors from C arrays + std::vector coord_(coord, coord+natoms*3); + std::vector atype_(atype, atype+natoms); + std::vector cell_; + if (cell) { + // pbc + cell_.assign(cell, cell+9); + } + // pairs + std::vector > pairs_; + for (int i = 0; i < npairs; i++) { + pairs_.push_back(std::make_pair(pairs[i*2], pairs[i*2+1])); + } + std::vector delef_(delef, delef+natoms*3); + std::vector df, dv; + + dcm->dcm.compute(df, dv, coord_, atype_, cell_, pairs_, delef_, nghost, nlist->nl); + // copy from C++ vectors to C arrays, if not NULL pointer + if(dfcorr_) std::copy(df.begin(), df.end(), dfcorr_); + if(dvcorr_) std::copy(dv.begin(), dv.end(), dvcorr_); +} + +template +void DP_DipoleChargeModifierComputeNList_variant ( + DP_DipoleChargeModifier* dcm, + const int natoms, + const double* coord, + const int* atype, + const double* cell, + const int* pairs, + const int npairs, + const double* delef, + const int nghost, + const DP_Nlist* nlist, + double* dfcorr_, + double* dvcorr_ + ); + +template +void DP_DipoleChargeModifierComputeNList_variant ( + DP_DipoleChargeModifier* dcm, + const int natoms, + const float* coord, + const int* atype, + const float* cell, + const int* pairs, + const int npairs, + const float* delef, + const int nghost, + const DP_Nlist* nlist, + float* dfcorr_, + float* dvcorr_ + ); + + extern "C" { void DP_DeepPotCompute ( @@ -828,6 +937,64 @@ int DP_DeepTensorGetNumbSelTypes( return dt->dt.sel_types().size(); } +void DP_DipoleChargeModifierComputeNList ( + DP_DipoleChargeModifier* dcm, + const int natom, + const double* coord, + const int* atype, + const double* cell, + const int* pairs, + const int npairs, + const double* delef_, + const int nghost, + const DP_Nlist* nlist, + double* dfcorr_, + double* dvcorr_ + ){ + DP_DipoleChargeModifierComputeNList_variant(dcm, natom, coord, atype, cell, pairs, npairs, delef_, nghost, nlist, dfcorr_, dvcorr_); +} + +void DP_DipoleChargeModifierComputeNListf ( + DP_DipoleChargeModifier* dcm, + const int natom, + const float* coord, + const int* atype, + const float* cell, + const int* pairs, + const int npairs, + const float* delef_, + const int nghost, + const DP_Nlist* nlist, + float* dfcorr_, + float* dvcorr_ + ){ + DP_DipoleChargeModifierComputeNList_variant(dcm, natom, coord, atype, cell, pairs, npairs, delef_, nghost, nlist, dfcorr_, dvcorr_); +} + +double DP_DipoleChargeModifierGetCutoff( + DP_DipoleChargeModifier* dcm + ) { + return dcm->dcm.cutoff(); +} + +int DP_DipoleChargeModifierGetNumbTypes( + DP_DipoleChargeModifier* dcm + ) { + return dcm->dcm.numb_types(); +} + +int* DP_DipoleChargeModifierGetSelTypes( + DP_DipoleChargeModifier* dcm + ) { + return (int*) &(dcm->dcm.sel_types())[0]; +} + +int DP_DipoleChargeModifierGetNumbSelTypes( + DP_DipoleChargeModifier* dcm + ) { + return dcm->dcm.sel_types().size(); +} + void DP_ConvertPbtxtToPb( const char* c_pbtxt, const char* c_pb diff --git a/source/api_c/tests/CMakeLists.txt b/source/api_c/tests/CMakeLists.txt index 7d864ee295..5abbe8b6e6 100644 --- a/source/api_c/tests/CMakeLists.txt +++ b/source/api_c/tests/CMakeLists.txt @@ -11,7 +11,7 @@ set_target_properties( add_executable( runUnitTests_c ${TEST_SRC} ) target_link_libraries(runUnitTests_c PRIVATE GTest::gtest_main ${LIB_DEEPMD_C} rt coverage_config) -target_link_libraries(runUnitTests_c PRIVATE ${LIB_DEEPMD}) +target_link_libraries(runUnitTests_c PRIVATE ${LIB_DEEPMD} ${LIB_DEEPMD_CC}) target_precompile_headers(runUnitTests_c PRIVATE test_utils.h [["deepmd.hpp"]]) add_test( runUnitTests_c runUnitTests_c ) set_target_properties( diff --git a/source/api_c/tests/test_dipolecharge.cc b/source/api_c/tests/test_dipolecharge.cc new file mode 100644 index 0000000000..db15d62334 --- /dev/null +++ b/source/api_c/tests/test_dipolecharge.cc @@ -0,0 +1,242 @@ +#include +#include +#include +#include +#include +#include "deepmd.hpp" +#include "test_utils.h" +#include "common.h" +#include "ewald.h" +#include "region.h" + +template +class TestDipoleCharge : public ::testing::Test +{ +protected: + std::vector coord = { + 4.6067455554, 8.8719311819, 6.3886531197, + 4.0044515745, 4.2449530507, 7.7902855220, + 2.6453069446, 0.8772647726, 1.2804446790, + 1.1445332290, 0.0067366438, 1.8606485070, + 7.1002867706, 5.0325506787, 3.1805888348, + 4.5352891138, 7.7389683929, 9.4260970128, + 2.1833238914, 9.0916071034, 7.2299906064, + 4.1040157820, 1.0496745045, 5.4748315591, + }; + std::vector atype = { + 0,3,2,1,3,4,1,4 + }; + std::vector box = { + 10., 0., 0., 0., 10., 0., 0., 0., 10. + }; + std::vector expected_e = { + 3.671081837126222158e+00 + }; + std::vector expected_f = { + 8.786854427753210128e-01,-1.590752486903602159e-01,-2.709225006303785932e-01,-4.449513960033193438e-01,-1.564291540964127813e-01,2.139031741772115178e-02,1.219699614140521193e+00,-5.580358618499958734e-02,-3.878662478349682585e-01,-1.286685244990778854e+00,1.886475802950296488e-01,3.904450515493615437e-01,1.605017382138404849e-02,2.138016869742287995e-01,-2.617514921203008965e-02,2.877081057057793712e-01,-3.846449683844421763e-01,3.048855616906603894e-02,-9.075632811311897807e-01,-6.509653472431625731e-03,2.302010972126376787e-01,2.370565856822822726e-01,3.600133435593881881e-01,1.243887532859055609e-02 + }; + std::vector expected_v = { + 3.714071471995848417e-01,6.957130186032146613e-01,-1.158289779017217302e+00,6.957130186032139951e-01,-1.400130091653774933e+01,-3.631620234653316626e-01,-1.158289779017217302e+00,-3.631620234653316626e-01,3.805077486043773050e+00 + }; + std::vector charge_map = { + 1., 1., 1., 1., 1., -1., -3. + }; + int natoms; + int ntypes; + std::vector type_asso; + double expected_tot_e; + std::vectorexpected_tot_v; + + deepmd::hpp::DeepTensor dp; + deepmd::hpp::DipoleChargeModifier dm; + + void SetUp() override { + std::string file_name = "../../tests/infer/dipolecharge_e.pbtxt"; + std::string model = "dipolecharge_e.pb"; + deepmd::hpp::convert_pbtxt_to_pb(file_name, model); + dp.init(model, 0, "dipole_charge"); + dm.init(model, 0, "dipole_charge"); + + natoms = atype.size(); + ntypes = 5; + type_asso.resize(ntypes, -1); + type_asso[1] = 5; + type_asso[3] = 6; + + EXPECT_EQ(natoms * 3, expected_f.size()); + EXPECT_EQ(9, expected_v.size()); + }; + + void TearDown() override { + remove( "dipolecharge_e.pb" ) ; + }; +}; + +static bool +_in_vec(const int & value, + const std::vector & vec) +{ + // naive impl. + for(int ii = 0; ii < vec.size(); ++ii){ + if(value == vec[ii]) return true; + } + return false; +} + +TYPED_TEST_SUITE(TestDipoleCharge, ValueTypes); + +TYPED_TEST(TestDipoleCharge, cpu_lmp_nlist) +{ + using VALUETYPE = TypeParam; + std::vector& coord = this->coord; + std::vector& atype = this->atype; + std::vector& box = this->box; + std::vector& expected_e = this->expected_e; + std::vector& expected_f = this->expected_f; + std::vector& expected_v = this->expected_v; + std::vector& charge_map = this->charge_map; + int& natoms = this->natoms; + int& ntypes = this->ntypes; + std::vector& type_asso = this->type_asso; + double& expected_tot_e = this->expected_tot_e; + std::vector&expected_tot_v = this->expected_tot_v; + deepmd::hpp::DeepTensor& dp = this->dp; + deepmd::hpp::DipoleChargeModifier& dm = this->dm; + // build nlist + // float rc = dp.cutoff(); + float rc = 4.0; + int nloc = coord.size() / 3; + std::vector coord_cpy; + std::vector atype_cpy, mapping; + std::vector > nlist_data; + _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping, + coord, atype, box, rc); + int nall = coord_cpy.size() / 3; + int nghost = nall - nloc; + std::vector ilist(nloc), numneigh(nloc); + std::vector firstneigh(nloc); + deepmd::hpp::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]); + convert_nlist(inlist, nlist_data); + + // evaluate dipole + std::vector dipole, dipole_recd(nloc*3, 0.0); + dp.compute(dipole, coord_cpy, atype_cpy, box, nall-nloc, inlist); + + // add virtual atoms to the system + // // a lot of mappings + std::vector sel_types = dp.sel_types(); + std::vector sel_fwd, sel_bwd; + int sel_nghost; + deepmd::select_by_type(sel_fwd, sel_bwd, sel_nghost, coord_cpy, atype_cpy, nghost, sel_types); + int sel_nall = sel_bwd.size(); + int sel_nloc = sel_nall - sel_nghost; + std::vector sel_atype(sel_bwd.size()); + deepmd::select_map(sel_atype, atype, sel_fwd, 1); + // Yixiao: because the deeptensor already return the correct order, the following map is no longer needed + // deepmd::AtomMap nnp_map(sel_atype.begin(), sel_atype.begin() + sel_nloc); + // const std::vector & sort_fwd_map(nnp_map.get_fwd_map()); + + // // add coords + std::vector add_coord; + std::vector add_atype; + std::vector> pairs; + for(int ii = 0; ii < nloc; ++ii){ + if(_in_vec(atype[ii], sel_types)){ + // Yixiao: the sort map is no longer needed + // int res_idx = sort_fwd_map[sel_fwd[ii]]; + int res_idx = sel_fwd[ii]; + std::vector tmp_coord(3); + for(int dd = 0; dd < 3; ++dd){ + tmp_coord[dd] = coord[ii*3+dd] + dipole[res_idx*3+dd]; + dipole_recd[ii*3+dd] = dipole[res_idx*3+dd]; + } + pairs.push_back(std::pair(ii, add_atype.size()+atype.size())); + // std::cout << ii << " " + // << atype[ii] << " " + // << res_idx << " " + // << type_asso[atype[ii]] << " " + // << " pair " + // << pairs.back().first << " " << pairs.back().second << " " + // << std::endl; + add_coord.insert(add_coord.end(), tmp_coord.begin(), tmp_coord.end()); + add_atype.push_back(type_asso[atype[ii]]); + } + } + coord.insert(coord.end(), add_coord.begin(), add_coord.end()); + atype.insert(atype.end(), add_atype.begin(), add_atype.end()); + nloc = atype.size(); + EXPECT_EQ(atype.size()*3, coord.size()); + + // get charge value + std::vector charge(nloc); + for(int ii = 0; ii < nloc; ++ii){ + charge[ii] = charge_map[atype[ii]]; + } + + // compute the recp part of the ele interaction + VALUETYPE eener; + std::vector eforce, evirial; + deepmd::Region region; + init_region_cpu(region, &box[0]); + deepmd::EwaldParameters eparam; + eparam.beta = 0.2; + eparam.spacing = 4; + ewald_recp(eener, eforce, evirial, coord, charge, region, eparam); + + EXPECT_LT(fabs(eener - expected_e[0]), 1e-6); + EXPECT_EQ(eforce.size(), coord.size()); + EXPECT_EQ(evirial.size(), 9); + + // extend the system with virtual atoms, and build nlist + _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping, + coord, atype, box, rc); + nall = coord_cpy.size() / 3; + nghost = nall - nloc; + ilist.resize(nloc); + numneigh.resize(nloc); + firstneigh.resize(nloc); + inlist.inum = nloc; + inlist.ilist = &ilist[0]; + inlist.numneigh = &numneigh[0]; + inlist.firstneigh = &firstneigh[0]; + convert_nlist(inlist, nlist_data); + + // compute force and virial + std::vector force_, force, virial; + dm.compute(force_, virial, coord_cpy, atype_cpy, box, pairs, eforce, nghost, inlist); + // for(int ii = 0; ii < force_.size(); ++ii){ + // std::cout << force_[ii] << " " ; + // } + // std::cout << std::endl; + _fold_back(force, force_, mapping, nloc, nall, 3); + + // compare force + EXPECT_EQ(force.size(), nloc*3); + // note nloc > expected_f.size(), because nloc contains virtual atoms. + for(int ii = 0; ii < expected_f.size(); ++ii){ + EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-6); + } + + // add recp virial and viral corr to virial + // virial = virial_recp + virial_dipolecharge + virial_corr + for (int dd0 = 0; dd0 < 3; ++dd0){ + for (int dd1 = 0; dd1 < 3; ++dd1){ + virial[dd0*3+dd1] += evirial[dd0*3+dd1]; + } + } + for(int ii = 0; ii < pairs.size(); ++ii){ + int idx0 = pairs[ii].first; + int idx1 = pairs[ii].second; + for (int dd0 = 0; dd0 < 3; ++dd0){ + for (int dd1 = 0; dd1 < 3; ++dd1){ + virial[dd0*3+dd1] -= eforce[idx1*3+dd0] * dipole_recd[idx0*3+dd1]; + } + } + } + // compare virial + EXPECT_EQ(virial.size(), 3*3); + for(int ii = 0; ii < expected_v.size(); ++ii){ + EXPECT_LT(fabs(virial[ii] - expected_v[ii]), 1e-5); + } +} + diff --git a/source/api_cc/include/DataModifier.h b/source/api_cc/include/DataModifier.h index 9c936939a1..50cba7fb5d 100644 --- a/source/api_cc/include/DataModifier.h +++ b/source/api_cc/include/DataModifier.h @@ -3,19 +3,53 @@ #include "DeepPot.h" namespace deepmd{ +/** +* @brief Dipole charge modifier. +**/ class DipoleChargeModifier { public: + /** + * @brief Dipole charge modifier without initialization. + **/ DipoleChargeModifier(); + /** + * @brief Dipole charge modifier without initialization. + * @param[in] model The name of the frozen model file. + * @param[in] gpu_rank The GPU rank. Default is 0. + * @param[in] name_scope The name scope. + **/ DipoleChargeModifier(const std::string & model, const int & gpu_rank = 0, const std::string & name_scope = ""); ~DipoleChargeModifier (); + /** + * @brief Initialize the dipole charge modifier. + * @param[in] model The name of the frozen model file. + * @param[in] gpu_rank The GPU rank. Default is 0. + * @param[in] name_scope The name scope. + **/ void init (const std::string & model, const int & gpu_rank = 0, const std::string & name_scope = ""); + /** + * @brief Print the DP summary to the screen. + * @param[in] pre The prefix to each line. + **/ void print_summary(const std::string &pre) const; public: + /** + * @brief Evaluate the force and virial correction by using this dipole charge modifier. + * @param[out] dfcorr_ The force correction on each atom. + * @param[out] dvcorr_ The virial correction. + * @param[in] dcoord_ The coordinates of atoms. The array should be of size natoms x 3. + * @param[in] datype_ The atom types. The list should contain natoms ints. + * @param[in] dbox The cell of the region. The array should be of size 9. + * @param[in] pairs The pairs of atoms. The list should contain npairs pairs of ints. + * @param[in] delef_ The electric field on each atom. The array should be of size natoms x 3. + * @param[in] nghost The number of ghost atoms. + * @param[in] lmp_list The neighbor list. + **/ template void compute (std::vector & dfcorr_, std::vector & dvcorr_, @@ -26,8 +60,20 @@ class DipoleChargeModifier const std::vector & delef_, const int nghost, const InputNlist & lmp_list); + /** + * @brief Get cutoff radius. + * @return double cutoff radius. + */ double cutoff () const {assert(inited); return rcut;}; + /** + * @brief Get the number of atom types. + * @return int number of atom types. + */ int numb_types () const {assert(inited); return ntypes;}; + /** + * @brief Get the list of sel types. + * @return The list of sel types. + */ std::vector sel_types () const {assert(inited); return sel_type;}; private: tensorflow::Session* session; diff --git a/source/api_cc/include/DeepTensor.h b/source/api_cc/include/DeepTensor.h index e21f3d25fb..dc02f41bef 100644 --- a/source/api_cc/include/DeepTensor.h +++ b/source/api_cc/include/DeepTensor.h @@ -162,6 +162,10 @@ class DeepTensor * @return The output dimension. **/ int output_dim () const {assert(inited); return odim;}; + /** + * @brief Get the list of sel types. + * @return The list of sel types. + */ const std::vector & sel_types () const {assert(inited); return sel_type;}; private: tensorflow::Session* session; diff --git a/source/api_cc/src/DataModifier.cc b/source/api_cc/src/DataModifier.cc index 06c04e3a79..5dd5f0e1fd 100644 --- a/source/api_cc/src/DataModifier.cc +++ b/source/api_cc/src/DataModifier.cc @@ -17,7 +17,7 @@ DipoleChargeModifier(const std::string & model, : inited (false), name_scope(name_scope_), graph_def(new GraphDef()) { - init(model, gpu_rank); + init(model, gpu_rank, name_scope_); } DipoleChargeModifier:: diff --git a/source/api_cc/src/DeepTensor.cc b/source/api_cc/src/DeepTensor.cc index 92efda3a67..427e77d067 100644 --- a/source/api_cc/src/DeepTensor.cc +++ b/source/api_cc/src/DeepTensor.cc @@ -17,7 +17,7 @@ DeepTensor(const std::string & model, : inited (false), name_scope(name_scope_), graph_def(new GraphDef()) { - init(model, gpu_rank); + init(model, gpu_rank, name_scope_); } DeepTensor::~DeepTensor() {