diff --git a/source/md/CMakeLists.txt b/source/md/CMakeLists.txt deleted file mode 100644 index dfd1c547be..0000000000 --- a/source/md/CMakeLists.txt +++ /dev/null @@ -1,82 +0,0 @@ -# md -set(MAKE_FF_AD FALSE) - -list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake/) -find_package(xdrfile REQUIRED) - -list(APPEND MD_INCLUDE_PATH "include") -list(APPEND MD_INCLUDE_PATH ${XDRFILE_INCLUDE_DIRS}) - -file(GLOB MD_SRC src/*.cc src/*.cpp) - -set(MDNN_SOURCE_FILES mdnn.cc) -if(MAKE_FF_AD) - set(MDAD_SOURCE_FILES mdad.cc) - set(MDFF_SOURCE_FILES mdff.cc) -endif() - -function(_add_md_variant variant_name prec_def) - set(libname "${LIB_DEEPMD_NATIVE}${variant_name}") - set(dp_mdnn_name "dp_mdnn${variant_name}") - set(dp_mdff_name "dp_mdff${variant_name}") - set(dp_mdad_name "dp_mdad${variant_name}") - - add_library(${libname} SHARED ${MD_SRC}) - target_link_libraries(${libname} PRIVATE ${LIB_DEEPMD}) - target_include_directories(${libname} PUBLIC ${MD_INCLUDE_PATH}) - set_target_properties(${libname} PROPERTIES COMPILE_DEFINITIONS ${prec_def} - INSTALL_RPATH "$ORIGIN") - - add_executable(${dp_mdnn_name} ${MDNN_SOURCE_FILES}) - if(MAKE_FF_AD) - add_executable(${dp_mdff_name} ${MDFF_SOURCE_FILES}) - add_executable(${dp_mdad_name} ${MDAD_SOURCE_FILES}) - endif() - - # link: libdeepmd_native libdeepmd_cc libxdr - target_link_libraries( - ${dp_mdnn_name} PRIVATE ${libname} ${LIB_DEEPMD_CC}${variant_name} - ${XDRFILE_LIBRARIES}) - target_include_directories(${dp_mdnn_name} - PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../3rdparty/) - if(MAKE_FF_AD) - target_link_libraries( - ${dp_mdad_name} PRIVATE ${libname} ${LIB_DEEPMD_CC}${variant_name} - ${XDRFILE_LIBRARIES}) - target_include_directories(${dp_mdad_name} - PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../3rdparty/) - target_link_libraries( - ${dp_mdff_name} PRIVATE ${libname} ${LIB_DEEPMD_CC}${variant_name} - ${XDRFILE_LIBRARIES}) - target_include_directories(${dp_mdff_name} - PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../3rdparty/) - endif() - - set_target_properties( - ${dp_mdnn_name} - PROPERTIES COMPILE_DEFINITIONS ${prec_def} - LINK_FLAGS "-Wl,-rpath,'$ORIGIN'/../lib -Wl,-z,defs" - INSTALL_RPATH "$ORIGIN/../lib:${TensorFlow_LIBRARY_PATH}") - if(MAKE_FF_AD) - set_target_properties( - ${dp_mdad_name} - PROPERTIES COMPILE_DEFINITIONS ${prec_def} - LINK_FLAGS "-Wl,-rpath,'$ORIGIN'/../lib -Wl,-z,defs" - INSTALL_RPATH "$ORIGIN/../lib:${TensorFlow_LIBRARY_PATH}") - set_target_properties( - ${dp_mdff_name} - PROPERTIES COMPILE_DEFINITIONS ${prec_def} - LINK_FLAGS "-Wl,-rpath,'$ORIGIN'/../lib -Wl,-z,defs" - INSTALL_RPATH "$ORIGIN/../lib:${TensorFlow_LIBRARY_PATH}") - endif() - - install(TARGETS ${LIB_DEEPMD_NATIVE} DESTINATION lib/) - install(TARGETS ${dp_mdnn_name} DESTINATION bin/) - if(MAKE_FF_AD) - install(TARGETS ${dp_mdad_name} DESTINATION bin/) - install(TARGETS ${dp_mdff_name} DESTINATION bin/) - endif() -endfunction() -_add_md_variant("${HIGH_PREC_VARIANT}" "${HIGH_PREC_DEF}") -# TODO: there is hard-code `DOUBLE` in the code -# _add_md_variant("${LOW_PREC_VARIANT}" "${LOW_PREC_DEF}") diff --git a/source/md/include/AdWeight.h b/source/md/include/AdWeight.h deleted file mode 100644 index 921185bf9f..0000000000 --- a/source/md/include/AdWeight.h +++ /dev/null @@ -1,63 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include - -#ifdef HIGH_PREC -typedef double VALUETYPE; -#else -typedef float VALUETYPE; -#endif - -class AdWeight { - public: - AdWeight(const VALUETYPE& pl); - virtual void zone_tag(std::vector& tag, - const std::vector& coord) const = 0; - virtual void atom_weight(std::vector& weight, - std::vector& weight_x, - const std::vector& coord) const = 0; - virtual std::vector get_center() const = 0; - void sel_nn_atom(std::vector& nn_coord, - std::vector& nn_type, - std::vector& nn_idx, - std::vector& nn_tag, - const std::vector& dcoord, - const std::vector& dtype) const; - void force_intpl(std::vector& of, - const std::vector& dcoord, - const std::vector& ff_force, - const std::vector& nn_force, - const std::vector& nn_idx) const; - void force_intpl(std::vector& of, - const std::vector& dcoord, - const std::vector& ff_bd_force, - const std::vector& ff_nb_force, - const std::vector& nn_force, - const std::vector& nn_idx) const; - - private: - VALUETYPE protect_level; -}; - -// slab model, axis x -class SlabWeight : public AdWeight { - public: - SlabWeight(const std::vector& box, - const VALUETYPE& rnn, - const VALUETYPE& rhy, - const VALUETYPE& rc, - const VALUETYPE& protect_level = 1e-3); - virtual void zone_tag(std::vector& tag, - const std::vector& coord) const; - virtual void atom_weight(std::vector& weight, - std::vector& weight_x, - const std::vector& coord) const; - virtual std::vector get_center() const { return center; } - - private: - std::vector center; - VALUETYPE rnn; - VALUETYPE rhy; - VALUETYPE rc; -}; diff --git a/source/md/include/Convert.h b/source/md/include/Convert.h deleted file mode 100644 index a0a11884e3..0000000000 --- a/source/md/include/Convert.h +++ /dev/null @@ -1,40 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include -#include -#include - -template -class Convert { - public: - Convert(const std::vector& atomname, - std::map& name_type_map, - std::map& name_mass_map, - std::map& name_charge_map, - const bool sort = true); - void gro2nnp(std::vector& coord, - std::vector& veloc, - std::vector& box, - const std::vector >& posi, - const std::vector >& velo, - const std::vector& box_size) const; - void nnp2gro(std::vector >& posi, - std::vector >& velo, - std::vector& box_size, - const std::vector& coord, - const std::vector& veloc, - const std::vector& box) const; - void idx_gro2nnp(std::vector& out, const std::vector& in) const; - void idx_nnp2gro(std::vector& out, const std::vector& in) const; - const std::vector& get_type() const { return atype; } - const std::vector& get_mass() const { return amass; } - const std::vector& get_charge() const { return acharge; } - - private: - std::vector idx_map_nnp2gro; - std::vector idx_map_gro2nnp; - std::vector atype; - std::vector amass; - std::vector acharge; -}; diff --git a/source/md/include/CosSwitch.h b/source/md/include/CosSwitch.h deleted file mode 100644 index 150cf1a8ba..0000000000 --- a/source/md/include/CosSwitch.h +++ /dev/null @@ -1,57 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once -#include - -#ifdef HIGH_PREC -typedef double VALUETYPE; -#else -typedef float VALUETYPE; -#endif - -class CosSwitch { - public: - CosSwitch(const VALUETYPE& rmin_ = 0, const VALUETYPE& rmax_ = 0) { - reinit(rmin_, rmax_); - } - void reinit(const VALUETYPE& rmin_, const VALUETYPE& rmax_); - - public: - void eval(VALUETYPE& vv, const VALUETYPE xx) const; - - private: - VALUETYPE rmin, rmax; -}; - -void CosSwitch::reinit(const VALUETYPE& rmin_, const VALUETYPE& rmax_) { - rmin = rmin_; - rmax = rmax_; -} - -void CosSwitch::eval(VALUETYPE& vv, const VALUETYPE xx) const { - VALUETYPE dd; - if (xx >= 0) { - if (xx < rmin) { - dd = 0; - vv = 1; - } else if (xx < rmax) { - VALUETYPE value = (xx - rmin) / (rmax - rmin) * M_PI; - dd = -0.5 * sin(value) * M_PI / (rmax - rmin); - vv = 0.5 * (cos(value) + 1); - } else { - dd = 0; - vv = 0; - } - } else { - if (xx > -rmin) { - dd = 0; - vv = 1; - } else if (xx > -rmax) { - VALUETYPE value = (-xx - rmin) / (rmax - rmin) * M_PI; - dd = 0.5 * sin(value) * M_PI / (rmax - rmin); - vv = 0.5 * (cos(value) + 1); - } else { - dd = 0; - vv = 0; - } - } -} diff --git a/source/md/include/Gaussian.h b/source/md/include/Gaussian.h deleted file mode 100644 index c10cc11c2a..0000000000 --- a/source/md/include/Gaussian.h +++ /dev/null @@ -1,13 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include -#include - -#include "RandomGenerator.h" - -class Gaussian { - public: - void set_seed(unsigned long seed); - void gen(double* vec, const int numb_gen); -}; diff --git a/source/md/include/GroFileManager.h b/source/md/include/GroFileManager.h deleted file mode 100644 index f80c3b57fe..0000000000 --- a/source/md/include/GroFileManager.h +++ /dev/null @@ -1,50 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#ifndef __GroFileManager_wanghan__ -#define __GroFileManager_wanghan__ - -#include -#include -#include - -namespace GroFileManager { -void read(const std::string& name, - std::vector& resdindex, - std::vector& resdname, - std::vector& atomname, - std::vector& atomindex, - std::vector >& posi, - std::vector >& velo, - std::vector& boxsize); -void write(const std::string& name, - const std::vector& resdindex, - const std::vector& resdname, - const std::vector& atomname, - const std::vector& atomindex, - const std::vector >& posi, - const std::vector >& velo, - const std::vector& boxsize); - -bool readTop(const std::string& filename, - std::vector& molnames, - std::vector& nmols); - -template -bool writePotenFile(const double& rmin, - const double& rcut, - const double& interval, - UnitaryFunction1& f, - UnitaryFunction2& fp, - UnitaryFunction3& g, - UnitaryFunction4& gp, - UnitaryFunction5& h, - UnitaryFunction6& hp, - const std::string& filename); - -}; // namespace GroFileManager - -#endif diff --git a/source/md/include/HarmonicAngle.h b/source/md/include/HarmonicAngle.h deleted file mode 100644 index b20430f599..0000000000 --- a/source/md/include/HarmonicAngle.h +++ /dev/null @@ -1,30 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include - -#include "SimulationRegion.h" - -#ifdef HIGH_PREC -typedef double VALUETYPE; -#else -typedef float VALUETYPE; -#endif - -class HarmonicAngle { - public: - HarmonicAngle(const VALUETYPE& kk, const VALUETYPE& tt); - - public: - void compute(VALUETYPE& ener, - std::vector& force, - std::vector& virial, - const std::vector& coord, - const std::vector& atype, - const SimulationRegion& region, - const std::vector& alist); - - private: - VALUETYPE ka, tt; - void hb_inner(VALUETYPE& ae, VALUETYPE& af, const VALUETYPE& r2); -}; diff --git a/source/md/include/HarmonicBond.h b/source/md/include/HarmonicBond.h deleted file mode 100644 index dd13491cc6..0000000000 --- a/source/md/include/HarmonicBond.h +++ /dev/null @@ -1,30 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include - -#include "SimulationRegion.h" - -#ifdef HIGH_PREC -typedef double VALUETYPE; -#else -typedef float VALUETYPE; -#endif - -class HarmonicBond { - public: - HarmonicBond(const VALUETYPE& kk, const VALUETYPE& bb); - - public: - void compute(VALUETYPE& ener, - std::vector& force, - std::vector& virial, - const std::vector& coord, - const std::vector& atype, - const SimulationRegion& region, - const std::vector& blist); - - private: - VALUETYPE kk, bb; - void hb_inner(VALUETYPE& ae, VALUETYPE& af, const VALUETYPE& r2); -}; diff --git a/source/md/include/Integrator.h b/source/md/include/Integrator.h deleted file mode 100644 index 73e84e7097..0000000000 --- a/source/md/include/Integrator.h +++ /dev/null @@ -1,50 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include - -#include "Gaussian.h" -#include "UnitManager.h" - -template -class Integrator { - public: - Integrator() : massConst(UnitManager::IntegratorMassConstant){}; - - public: - void stepVeloc(std::vector& vv, - const std::vector& ff, - const std::vector& mass, - const double& dt, - const std::vector& freez = std::vector()) const; - void stepCoord(std::vector& rr, - const std::vector& vv, - const double& dt) const; - - private: - VALUETYPE massConst; -}; - -template -class ThermostatLangevin { - public: - ThermostatLangevin(const VALUETYPE T = 300., - const VALUETYPE tau = 1., - const long long int seed = 0); - void reinit(const VALUETYPE T = 300., - const VALUETYPE tau = 1., - const long long int seed = 0); - void stepOU(std::vector& vv, - const std::vector& mass, - const double& dt, - const std::vector& freez = std::vector()) const; - - private: - mutable Gaussian gaussian; - std::string scheme; - VALUETYPE temperature; - VALUETYPE gamma; - VALUETYPE sigma; - VALUETYPE kT; - VALUETYPE sigmainvsqrt2gamma; -}; diff --git a/source/md/include/Interpolation.h b/source/md/include/Interpolation.h deleted file mode 100644 index ad64d6114e..0000000000 --- a/source/md/include/Interpolation.h +++ /dev/null @@ -1,70 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#ifndef __wanghan__Interpolation_h__ -#define __wanghan__Interpolation_h__ - -#include -#include -#include -#include - -#include "Poly.h" - -namespace Interpolation { -// linear interpolations -void pieceLinearInterpol(const double& a, - const double& b, - const double& va, - const double& vb, - Poly& p); -void piecewiseLinear(const std::vector& x, - const std::vector& y, - PiecewisePoly& ps); -// spline interpolations -void pieceHermiteInterpol(const double& a, - const double& b, - const double& va, - const double& vb, - const double& da, - const double& db, - Poly& p); -void pieceSecondDerivativeInterpol(const double& a, - const double& b, - const double& va, - const double& vb, - const double& dda, - const double& ddb, - Poly& p); -void piece6OrderInterpol(const double& a, - const double& b, - const double& va, - const double& vb, - const double& da, - const double& db, - const double& dda, - const double& ddb, - Poly& p); - -bool spline(const std::vector& x, - const std::vector& y, - PiecewisePoly& ps); -bool spline(const std::vector::const_iterator xbegin, - const std::vector::const_iterator xend, - const std::vector::const_iterator ybegin, - PiecewisePoly& ps); -bool splinePeriodic(const std::vector& x, - const std::vector& y, - PiecewisePoly& ps); -bool solverForSplinePeriodic(const std::vector::const_iterator& lbegin, - const std::vector::const_iterator& lend, - const std::vector::iterator& ubegin, - const std::vector::iterator& uend); -void secondDerivativeInterpol( - const std::vector::const_iterator& xbegin, - const std::vector::const_iterator& xend, - const std::vector::const_iterator& vbegin, - const std::vector::const_iterator& ddbegin, - PiecewisePoly& ps); - -} // namespace Interpolation - -#endif diff --git a/source/md/include/LJInter.h b/source/md/include/LJInter.h deleted file mode 100644 index fe507a0e90..0000000000 --- a/source/md/include/LJInter.h +++ /dev/null @@ -1,31 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include - -#include "SimulationRegion.h" - -#ifdef HIGH_PREC -typedef double VALUETYPE; -#else -typedef float VALUETYPE; -#endif - -class LJInter { - public: - LJInter(const VALUETYPE& c6, const VALUETYPE& c12, const VALUETYPE& rc); - - public: - void compute(VALUETYPE& ener, - std::vector& force, - std::vector& virial, - const std::vector& coord, - const std::vector& atype, - const SimulationRegion& region, - const std::vector >& nlist); - - private: - VALUETYPE c6, c12, rc, rc2, one_over_6, one_over_12, one_over_rc6, - one_over_rc12; - void lj_inner(VALUETYPE& ae, VALUETYPE& af, const VALUETYPE& r2); -}; diff --git a/source/md/include/LJTab.h b/source/md/include/LJTab.h deleted file mode 100644 index 73f46697ba..0000000000 --- a/source/md/include/LJTab.h +++ /dev/null @@ -1,31 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include - -#include "Tabulated.h" - -#ifdef HIGH_PREC -typedef double VALUETYPE; -#else -typedef float VALUETYPE; -#endif - -class LJTab { - public: - LJTab(const VALUETYPE& c6, const VALUETYPE& c12, const VALUETYPE& rc); - - public: - void compute(VALUETYPE& ener, - std::vector& force, - std::vector& virial, - const std::vector& coord, - const std::vector& atype, - const SimulationRegion& region, - const std::vector >& nlist) { - lj_tab.compute(ener, force, virial, coord, atype, region, nlist); - }; - - private: - Tabulated lj_tab; -}; diff --git a/source/md/include/MaxShift.h b/source/md/include/MaxShift.h deleted file mode 100644 index 21634df671..0000000000 --- a/source/md/include/MaxShift.h +++ /dev/null @@ -1,28 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include - -#include "SimulationRegion.h" - -#ifdef HIGH_PREC -typedef double VALUETYPE; -#else -typedef float VALUETYPE; -#endif - -class MaxShift { - public: - MaxShift(const std::vector& dcoord, const VALUETYPE& shell); - - bool rebuild(const std::vector& coord, - const SimulationRegion& region); - - private: - VALUETYPE - max_shift2(const std::vector& coord, - const SimulationRegion& region); - std::vector record; - VALUETYPE shell; - VALUETYPE max_allow2; -}; diff --git a/source/md/include/Poly.h b/source/md/include/Poly.h deleted file mode 100644 index 1b3c3d2e15..0000000000 --- a/source/md/include/Poly.h +++ /dev/null @@ -1,88 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#ifndef __wanghan_Poly_h__ -#define __wanghan_Poly_h__ - -#include -#include -#include - -class Poly { - std::vector a; - unsigned order; - - public: - Poly(); - Poly(const std::vector& out); - void reinit(const std::vector& out); - void zero() { - a.clear(); - a.resize(1, 0); - order = 0; - } - void one() { - a.clear(); - a.resize(1, 1); - order = 0; - } - - public: - Poly& operator=(const Poly& poly); - Poly& operator+=(const Poly& poly); - Poly& operator+=(const double& b); - Poly& operator*=(const Poly& poly); - Poly& operator*=(const double& scale); - Poly& derivative(); - - public: - unsigned& getOrder() { return order; } - const unsigned& getOrder() const { return order; } - std::vector& getCoeffs() { return a; } - const std::vector& getCoeffs() const { return a; } - - public: - void print(); - void print(const std::string& x); - void printCode(const std::string& x); - - public: - double value(const double& x) const; - - public: - // p = f(ax + b) - Poly& valueLinearPoly(const double& a, const double& b, Poly& p); -}; - -class PiecewisePoly { - public: - std::vector& get_x() { return x; } - std::vector& get_p() { return p; } - const std::vector& get_x() const { return x; } - const std::vector& get_p() const { return p; } - - public: - void clear() { - x.clear(); - p.clear(); - } - bool valid() const; - - public: - double value(const double& r) const; - void value(const std::vector& r, std::vector& y) const; - double value_periodic(const double& r) const; - void value_periodic(const std::vector& r, - std::vector& y) const; - - private: - std::vector x; - std::vector p; - void value(const unsigned& xbegin, - const unsigned& xend, - const std::vector& r, - const unsigned& rbegin, - const unsigned& rend, - std::vector& y) const; - double value(const double& xx, unsigned& begin, unsigned& end) const; -}; - -#endif diff --git a/source/md/include/RandomGenerator.h b/source/md/include/RandomGenerator.h deleted file mode 100644 index 0f14648e8d..0000000000 --- a/source/md/include/RandomGenerator.h +++ /dev/null @@ -1,13 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -namespace RandomGenerator_MT19937 { -void init_by_array(unsigned long init_key[], int key_length); -void init_genrand(unsigned long s); -unsigned long genrand_int32(void); -long genrand_int31(void); -double genrand_real1(void); // in [0,1] -double genrand_real2(void); // in [0,1) -double genrand_real3(void); // in (0,1) -double genrand_res53(void); -} // namespace RandomGenerator_MT19937 diff --git a/source/md/include/Statistics.h b/source/md/include/Statistics.h deleted file mode 100644 index 19e8896cbd..0000000000 --- a/source/md/include/Statistics.h +++ /dev/null @@ -1,40 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include - -#include "SimulationRegion.h" - -template -class Statistics { - public: - Statistics(const VALUETYPE e_corr = 0, const VALUETYPE p_corr = 0); - void record(const VALUETYPE& ener, - const std::vector& virial, - const std::vector& veloc, - const std::vector& mass, - const SimulationRegion& region); - - public: - double get_T() const; - double get_V() const; - double get_P() const; - double get_E() const { return get_ekin() + get_epot(); }; - double get_ekin() const { return r_kin_ener; } - double get_epot() const { return r_pot_ener + e_corr; } - - public: - void print(std::ostream& os, const int& step, const double time) const; - void print_head(std::ostream& os) const; - - private: - int natoms; - double r_ener; - double r_pot_ener; - double r_kin_ener; - // std::vector r_box; - SimulationRegion region; - std::vector r_vir; - double e_corr; - double p_corr; -}; diff --git a/source/md/include/StringSplit.h b/source/md/include/StringSplit.h deleted file mode 100644 index 9243e8ec6d..0000000000 --- a/source/md/include/StringSplit.h +++ /dev/null @@ -1,18 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#ifndef __StringSplit_h_wanghan__ -#define __StringSplit_h_wanghan__ - -#include -#include -#include -#include -#include - -namespace StringOperation { -void split(const std::string& in, std::vector& out); -void split(const std::string& in, - const std::string& delimiter, - std::vector& out); -} // namespace StringOperation - -#endif diff --git a/source/md/include/TF.h b/source/md/include/TF.h deleted file mode 100644 index be12a4f6ac..0000000000 --- a/source/md/include/TF.h +++ /dev/null @@ -1,30 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include -#include - -#include "AdWeight.h" -#include "common.h" - -#ifdef HIGH_PREC -typedef double VALUETYPE; -#else -typedef float VALUETYPE; -#endif - -class TF { - public: - TF(const std::string& filename); - - public: - void apply(std::vector& force, - const std::vector& coord, - const AdWeight& weight) const; - - private: - VALUETYPE meas(const VALUETYPE& xx) const; - std::vector data; - double hh; - double xup; -}; diff --git a/source/md/include/TableFileLoader.h b/source/md/include/TableFileLoader.h deleted file mode 100644 index 31461f5014..0000000000 --- a/source/md/include/TableFileLoader.h +++ /dev/null @@ -1,30 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#ifndef __TableFileLoader_h_wanghan__ -#define __TableFileLoader_h_wanghan__ - -#include -#include - -class TableFileLoader { - public: - unsigned getNumbColumns(); - - public: - TableFileLoader(const char* file); - void reinit(const char* file); - void setColumns(const std::vector& cols); - void setEvery(const unsigned every); - - public: - void loadAll(std::vector >& data); - bool loadLine(std::vector& data); - - private: - std::ifstream data; - std::string file; - unsigned count_read; - unsigned every; - std::vector inter_cols; -}; - -#endif diff --git a/source/md/include/Tabulated.h b/source/md/include/Tabulated.h deleted file mode 100644 index 5ab6e02bc3..0000000000 --- a/source/md/include/Tabulated.h +++ /dev/null @@ -1,46 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include - -#include "SimulationRegion.h" - -#ifdef HIGH_PREC -typedef double VALUETYPE; -#else -typedef float VALUETYPE; -#endif - -class Tabulated { - public: - Tabulated(){}; - Tabulated(const VALUETYPE rc, - const VALUETYPE hh, - const std::vector& tab); - void reinit(const VALUETYPE rc, - const VALUETYPE hh, - const std::vector& tab); - - public: - void compute(VALUETYPE& ener, - std::vector& force, - std::vector& virial, - const std::vector& coord, - const std::vector& atype, - const SimulationRegion& region, - const std::vector >& nlist); - void compute(VALUETYPE& ener, - std::vector& force, - std::vector& virial, - const std::vector& coord, - const std::vector& charge, - const std::vector& atype, - const SimulationRegion& region, - const std::vector >& nlist); - void tb_inner(VALUETYPE& ae, VALUETYPE& af, const VALUETYPE& r2); - - private: - VALUETYPE rc2, hi; - std::vector data; - void compute_posi(int& idx, VALUETYPE& eps, const VALUETYPE& rr); -}; diff --git a/source/md/include/Trajectory.h b/source/md/include/Trajectory.h deleted file mode 100644 index 862b393ea4..0000000000 --- a/source/md/include/Trajectory.h +++ /dev/null @@ -1,58 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#ifndef __MDFileManager_Trajectory_h_wanghan__ -#define __MDFileManager_Trajectory_h_wanghan__ - -// #include "Defines.h" -#include - -#include "xdrfile/xdrfile.h" -#include "xdrfile/xdrfile_trr.h" -#include "xdrfile/xdrfile_xtc.h" - -class XtcSaver { - public: - XtcSaver() : inited(false), prec(1000){}; - ~XtcSaver(); - XtcSaver(const char *filename, const int &natoms); - bool reinit(const char *filename, const int &natoms); - - public: - void save(const int &step, - const double &time, - const std::vector > &frame, - const std::vector &box); - - private: - XDRFILE *xd; - int natoms; - rvec *xx; - float prec; - bool inited; - void clear(); -}; - -class TrrSaver { - public: - TrrSaver() : inited(false), lambda(0){}; - ~TrrSaver(); - TrrSaver(const char *filename, const int &natoms); - bool reinit(const char *filename, const int &natoms); - - public: - void save(const int &step, - const double &time, - const std::vector > &ixx, - const std::vector > &ivv, - const std::vector > &iff, - const std::vector &box); - - private: - XDRFILE *xd; - int natoms; - rvec *xx, *vv, *ff; - float lambda; - bool inited; - void clear(); -}; - -#endif diff --git a/source/md/include/UnitManager.h b/source/md/include/UnitManager.h deleted file mode 100644 index 70393c406e..0000000000 --- a/source/md/include/UnitManager.h +++ /dev/null @@ -1,26 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include - -class UnitManager { - protected: - UnitManager(){}; - - public: - static double Degree2Radian; - static double Radian2Degree; - - static double IntegratorMassConstant; - static double PressureConstant; - static double BoltzmannConstant; - static double ElectrostaticConvertion; - - static double DefaultTableUpperLimit; - static double DefaultTableStep; - static double DefaultTableExtension; - static void set(const std::string& name_of_system); - - private: - static std::string unit_names[]; -}; diff --git a/source/md/include/XyzFileManager.h b/source/md/include/XyzFileManager.h deleted file mode 100644 index d557f0f771..0000000000 --- a/source/md/include/XyzFileManager.h +++ /dev/null @@ -1,20 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#ifndef __XyzFileManager_h_wanghan__ -#define __XyzFileManager_h_wanghan__ - -#include - -namespace XyzFileManager { - -void read(const std::string& file, - std::vector& atom_name, - std::vector >& posi, - std::vector >& velo, - std::vector >& forc, - std::vector& boxsize); - -void getBoxSize(const std::string& name, std::vector& boxsize); - -}; // namespace XyzFileManager - -#endif diff --git a/source/md/include/ZM.h b/source/md/include/ZM.h deleted file mode 100644 index 9255b0c17d..0000000000 --- a/source/md/include/ZM.h +++ /dev/null @@ -1,45 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include - -#include "SimulationRegion.h" -#include "Tabulated.h" -#include "ZMFunctions.h" - -#ifdef HIGH_PREC -typedef double VALUETYPE; -#else -typedef float VALUETYPE; -#endif - -class ZM { - public: - ZM(const int& order, const VALUETYPE& alpha, const VALUETYPE& rc); - - public: - void compute(VALUETYPE& ener, - std::vector& force, - std::vector& virial, - const std::vector& coord, - const std::vector& charge, - const std::vector& atype, - const SimulationRegion& region, - const std::vector >& nlist) { - zm_tab.compute(ener, force, virial, coord, charge, atype, region, nlist); - }; - void exclude(VALUETYPE& ener, - std::vector& force, - std::vector& virial, - const std::vector& coord, - const std::vector& charge, - const std::vector& atype, - const SimulationRegion& region, - const std::vector& elist); - VALUETYPE e_corr(const std::vector& charge) const; - - private: - Tabulated zm_tab; - void ex_inner(VALUETYPE& ae, VALUETYPE& af, const VALUETYPE& r2); - ZeroMultipole::Potential potzm; -}; diff --git a/source/md/include/ZMFunctions.h b/source/md/include/ZMFunctions.h deleted file mode 100644 index aba6ce34d9..0000000000 --- a/source/md/include/ZMFunctions.h +++ /dev/null @@ -1,38 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#ifndef __Functions_h_ZM_wanghan__ -#define __Functions_h_ZM_wanghan__ - -#include - -namespace ZeroMultipole { -double funcV(const double& alpha, const double& r); -double funcD1V(const double& alpha, const double& r); -double funcD2V(const double& alpha, const double& r); -double funcD3V(const double& alpha, const double& r); -double funcD4V(const double& alpha, const double& r); - -void calCoefficients(const int& ll, - const double& alpha, - const double& rc, - std::vector& coeff); - -class Potential { - double alpha, rc; - int ll; - std::vector coeff; - - public: - Potential(); - Potential(const int& ll, const double& alpha, const double& rc); - void reinit(const int& ll, const double& alpha, const double& rc); - double pot(const double& rr); - double ulpot(const double& rr); - double mpotp(const double& rr); - double mulpotp(const double& rr); - - public: - double energyCorr(const std::vector& charges) const; -}; -} // namespace ZeroMultipole - -#endif diff --git a/source/md/include/common.h b/source/md/include/common.h deleted file mode 100644 index f1662f1206..0000000000 --- a/source/md/include/common.h +++ /dev/null @@ -1,43 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -#include - -#include "SimulationRegion.h" - -const double b2m_l = 10; -const double b2m_e = 1.660539040e-21 / 1.602176621e-19; - -template -void clear(VALUETYPE& ener, - std::vector& force, - std::vector& virial) { - ener = 0; - std::fill(force.begin(), force.end(), 0.); - std::fill(virial.begin(), virial.end(), 0.); -} - -template -void normalize_coord(std::vector& coord, - const SimulationRegion& region) { - int natoms = coord.size() / 3; - for (int ii = 0; ii < natoms; ++ii) { - double phys[3]; - for (int dd = 0; dd < 3; ++dd) { - phys[dd] = coord[ii * 3 + dd]; - } - double inter[3]; - region.phys2Inter(inter, phys); - for (int dd = 0; dd < 3; ++dd) { - if (inter[dd] < 0) { - inter[dd] += 1.; - } else if (inter[dd] >= 1) { - inter[dd] -= 1.; - } - } - region.inter2Phys(phys, inter); - for (int dd = 0; dd < 3; ++dd) { - coord[ii * 3 + dd] = phys[dd]; - } - } -} diff --git a/source/md/include/mymath.h b/source/md/include/mymath.h deleted file mode 100644 index aaae92704e..0000000000 --- a/source/md/include/mymath.h +++ /dev/null @@ -1,73 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#pragma once - -template -inline VALUETYPE dot(const VALUETYPE* r0, const VALUETYPE* r1) { - return (r0[0] * r1[0] + r0[1] * r1[1] + r0[2] * r1[2]); -} - -template -inline TYPE dot(const TYPE& x0, - const TYPE& y0, - const TYPE& z0, - const TYPE& x1, - const TYPE& y1, - const TYPE& z1) { - return x0 * x1 + y0 * y1 + z0 * z1; -} - -template -inline VALUETYPE cos(const VALUETYPE* r0, const VALUETYPE* r1) { - double ip = dot(r0, r1); - double ip0 = dot(r0, r0); - double ip1 = dot(r1, r1); - double ip01 = ip0 * ip1; - - double cosval; - if (ip01 > 0) { - cosval = ip / sqrt(ip01); - } else { - cosval = 1.0; - } - if (cosval > 1.0) { - return 1.0; - } - if (cosval < -1.0) { - return -1.0; - } - return cosval; -} - -template -inline TYPE cos(const TYPE& x0, - const TYPE& y0, - const TYPE& z0, - const TYPE& x1, - const TYPE& y1, - const TYPE& z1) { - double dblx0 = (double)(x0); - double dblx1 = (double)(x1); - double dbly0 = (double)(y0); - double dbly1 = (double)(y1); - double dblz0 = (double)(z0); - double dblz1 = (double)(z1); - - double ip = dot(dblx0, dbly0, dblz0, dblx1, dbly1, dblz1); - double ip0 = dot(dblx0, dbly0, dblz0, dblx0, dbly0, dblz0); - double ip1 = dot(dblx1, dbly1, dblz1, dblx1, dbly1, dblz1); - double ip01 = ip0 * ip1; - - double cosval; - if (ip01 > 0) { - cosval = ip / sqrt(ip01); - } else { - cosval = 1.0; - } - if (cosval > 1.0) { - return 1.0; - } - if (cosval < -1.0) { - return -1.0; - } - return cosval; -} diff --git a/source/md/mdnn.cc b/source/md/mdnn.cc deleted file mode 100644 index 80287db891..0000000000 --- a/source/md/mdnn.cc +++ /dev/null @@ -1,219 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "Convert.h" -#include "DeepPot.h" -#include "GroFileManager.h" -#include "Integrator.h" -#include "Statistics.h" -#include "Trajectory.h" -#include "XyzFileManager.h" -#include "common.h" -#include "json.hpp" -using json = nlohmann::json; - -#include - -#ifdef HIGH_PREC -typedef double VALUETYPE; -#else -typedef float VALUETYPE; -#endif - -void print_vec(const vector& vec) { - int nloc = vec.size() / 3; - for (int kk = 0; kk < nloc; ++kk) { - for (int dd = 0; dd < 3; ++dd) { - cout << vec[kk * 3 + dd] << " \t "; - } - cout << endl; - } -} - -int main(int argc, char* argv[]) { - UnitManager::set("metal"); - - if (argc == 0) { - cerr << "usage " << endl; - cerr << argv[0] << " input_script " << endl; - return 1; - } - - ifstream fp(argv[1]); - json jdata; - fp >> jdata; - cout << "# using data base" << endl; - cout << setw(4) << jdata << endl; - - int nframes = 1; - - vector dcoord, dveloc, dbox, dmass; - vector dtype; - vector freez; - - // load_raw (dcoord, dtype, dbox); - // dveloc.resize(dcoord.size(), 0.); - string conf_format = jdata["conf_format"]; - string conf_file = jdata["conf_file"]; - vector resdindex, atomindex; - vector resdname, atomname; - vector > posi, velo, tmp_forc; - vector boxsize; - if (conf_format == "gro") { - GroFileManager::read(conf_file, resdindex, resdname, atomname, atomindex, - posi, velo, boxsize); - } else if (conf_format == "xyz") { - XyzFileManager::read(conf_file, atomname, posi, velo, tmp_forc, boxsize); - if (velo.size() == 0) { - for (unsigned ii = 0; ii < posi.size(); ++ii) { - velo.push_back(vector(3, 0.)); - } - } - // convert to nanometer - for (unsigned ii = 0; ii < posi.size(); ++ii) { - for (unsigned dd = 0; dd < 3; ++dd) { - posi[ii][dd] *= .1; - velo[ii][dd] *= .1; - } - } - for (unsigned dd = 0; dd < 9; ++dd) { - boxsize[dd] *= .1; - } - for (unsigned ii = 0; ii < posi.size(); ++ii) { - resdindex.push_back(ii + 1); - atomindex.push_back(ii + 1); - } - resdname = atomname; - } else { - cerr << "unknown conf file format: " << conf_format << endl; - return 1; - } - map name_type_map = jdata["atom_type"]; - map name_mass_map = jdata["atom_mass"]; - map name_charge_map; - if (jdata.find("atom_charge") == jdata.end()) { - for (map::iterator iter = name_mass_map.begin(); - iter != name_mass_map.end(); ++iter) { - name_charge_map[iter->first] = 0.; - } - } else { - map name_charge_map_tmp = jdata["atom_charge"]; - name_charge_map = name_charge_map_tmp; - } - if (jdata.find("freeze_atoms") != jdata.end()) { - freez = jdata["freeze_atoms"].get >(); - } - - // convert but do not sort - Convert cvt(atomname, name_type_map, name_mass_map, - name_charge_map, false); - cvt.gro2nnp(dcoord, dveloc, dbox, posi, velo, boxsize); - dtype = cvt.get_type(); - dmass = cvt.get_mass(); - - int nloc = dtype.size(); - SimulationRegion region; - region.reinitBox(&dbox[0]); - normalize_coord(dcoord, region); - - vector dforce(nloc * 3, 0.); - vector dae(nloc * 1, 0.); - vector dav(nloc * 9, 0.); - vector dvirial(9, 0.0); - VALUETYPE dener = 0; - - string graph_file = jdata["graph_file"]; - VALUETYPE dt = jdata["dt"]; - int nsteps = jdata["nsteps"]; - int nener = jdata["ener_freq"]; - int nxtc = jdata["xtc_freq"]; - int ntrr = jdata["trr_freq"]; - string ener_file = jdata["ener_file"]; - string xtc_file = jdata["xtc_file"]; - string trr_file = jdata["trr_file"]; - double temperature = jdata["T"]; - double tau_t = jdata["tau_T"]; - long long int seed = 0; - if (jdata.find("rand_seed") != jdata.end()) { - seed = jdata["rand_seed"]; - } - bool print_f = false; - if (jdata.find("print_force") != jdata.end()) { - print_f = jdata["print_force"]; - } - - Integrator inte; - ThermostatLangevin thm(temperature, tau_t, seed); - deepmd::DeepPot nnp(graph_file); - - Statistics st; - XtcSaver sxtc(xtc_file.c_str(), nloc); - TrrSaver strr(trr_file.c_str(), nloc); - - // compute force at step 0 - nnp.compute(dener, dforce, dvirial, dcoord, dtype, dbox); - // change virial to gromacs convention - for (int ii = 0; ii < 9; ++ii) { - dvirial[ii] *= -0.5; - } - st.record(dener, dvirial, dveloc, dmass, region); - ofstream efout(ener_file); - ofstream pforce; - if (print_f) { - pforce.open("force.out"); - } - st.print_head(efout); - st.print(efout, 0, 0); - - for (int ii = 0; ii < nsteps; ++ii) { - inte.stepVeloc(dveloc, dforce, dmass, 0.5 * dt, freez); - inte.stepCoord(dcoord, dveloc, 0.5 * dt); - thm.stepOU(dveloc, dmass, dt, freez); - inte.stepCoord(dcoord, dveloc, 0.5 * dt); - normalize_coord(dcoord, region); - nnp.compute(dener, dforce, dvirial, dae, dav, dcoord, dtype, dbox); - // change virial to gromacs convention - for (int ii = 0; ii < 9; ++ii) { - dvirial[ii] *= -0.5; - } - inte.stepVeloc(dveloc, dforce, dmass, 0.5 * dt, freez); - if ((ii + 1) % nener == 0) { - st.record(dener, dvirial, dveloc, dmass, region); - st.print(efout, ii + 1, (ii + 1) * dt); - efout.flush(); - } - if (nxtc > 0 && (ii + 1) % nxtc == 0) { - cvt.nnp2gro(posi, velo, boxsize, dcoord, dveloc, dbox); - sxtc.save(ii + 1, (ii + 1) * dt, posi, boxsize); - } - if (ntrr > 0 && (ii + 1) % ntrr == 0) { - cvt.nnp2gro(posi, velo, boxsize, dcoord, dveloc, dbox); - strr.save(ii + 1, (ii + 1) * dt, posi, velo, vector >(), - boxsize); - if (print_f) { - for (int jj = 0; jj < dforce.size(); ++jj) { - pforce << dforce[jj] << " "; - } - pforce << endl; - } - } - } - - cvt.nnp2gro(posi, velo, boxsize, dcoord, dveloc, dbox); - GroFileManager::write("out.gro", resdindex, resdname, atomname, atomindex, - posi, velo, boxsize); - // ofstream oxyz ("out.xyz"); - // oxyz << nloc << endl; - // oxyz << setprecision(12) ; - // for (int ii = 0; ii < dbox.size(); ++ii) { - // oxyz << dbox[ii] * 1 << " " ; - // } - // oxyz << endl; - // for (int ii = 0; ii < posi.size(); ++ii){ - // oxyz << atomname[ii] << " \t" ; - // for (int dd = 0; dd < 3; ++dd){ - // oxyz << posi[ii][dd] * 10 << " "; - // } - // oxyz << endl; - // } - - return 0; -} diff --git a/source/md/src/AdWeight.cc b/source/md/src/AdWeight.cc deleted file mode 100644 index f1b9922f94..0000000000 --- a/source/md/src/AdWeight.cc +++ /dev/null @@ -1,258 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "AdWeight.h" - -#include -#include -#include - -#include "CosSwitch.h" - -AdWeight::AdWeight(const VALUETYPE& pl) { protect_level = pl; } - -void AdWeight::sel_nn_atom(vector& nn_coord, - vector& nn_type, - vector& nn_idx, - vector& nn_tag, - const vector& dcoord, - const vector& dtype) const { - nn_coord.clear(); - nn_type.clear(); - nn_idx.clear(); - - vector& tag(nn_tag); - zone_tag(tag, dcoord); - for (int ii = 0; ii < tag.size(); ++ii) { - if (tag[ii] != 0) { - nn_coord.push_back(dcoord[3 * ii + 0]); - nn_coord.push_back(dcoord[3 * ii + 1]); - nn_coord.push_back(dcoord[3 * ii + 2]); - nn_type.push_back(dtype[ii]); - nn_idx.push_back(ii); - } - } -} - -void AdWeight::force_intpl(vector& of, - const vector& dcoord, - const vector& ff_force, - const vector& nn_force, - const vector& nn_idx) const { - int nall = dcoord.size() / 3; - - vector weight, weight_x; - atom_weight(weight, weight_x, dcoord); - assert(nall == weight.size()); - // for (unsigned ii = 0; ii < weight.size(); ++ii){ - // cout << ii << " " << weight[ii] << " " << dcoord[ii*3] << endl; - // } - - // cout << "of " << of.size() << endl; - // cout << "dcoord " << dcoord.size() << endl; - // cout << "ff_f " << ff_force.size() << endl; - // cout << "nn_f " << nn_force.size() << endl; - // cout << "nn_i " << nn_idx.size() << endl; - // cout << "w " << weight.size() << endl; - vector nn_sum(3, 0.); - vector ff_sum(3, 0.); - // for (int ii = 0; ii < ff_force.size() / 3; ++ii){ - // for (int dd = 0; dd < 3; ++dd){ - // ff_sum[dd] += ff_force[ii*3+dd]; - // } - // } - // for (int ii = 0; ii < nn_force.size() / 3; ++ii){ - // for (int dd = 0; dd < 3; ++dd){ - // nn_sum[dd] += nn_force[ii*3+dd]; - // } - // } - // cout << ff_sum[0] << " " << ff_sum[1] << " " << ff_sum[2] << " " - // <& of, - const vector& dcoord, - const vector& ff_bd_force, - const vector& ff_nb_force, - const vector& nn_force, - const vector& nn_idx) const { - int nall = dcoord.size() / 3; - - vector weight, weight_x; - atom_weight(weight, weight_x, dcoord); - assert(nall == weight.size()); - - vector nn_sum(3, 0.); - vector ff_sum(3, 0.); - - for (int ii = 0; ii < nn_idx.size(); ++ii) { - int idx = nn_idx[ii]; - for (int dd = 0; dd < 3; ++dd) { - // nn_sum[dd] += weight[idx] * nn_force[ii*3+dd]; - // nn_sum[dd] += 1 * nn_force[ii*3+dd]; - of[idx * 3 + dd] += weight[idx] * nn_force[ii * 3 + dd]; - // if (fabs(nn_force[ii*3+dd]) > 1e6) { - // cout << " ii " << ii - // << " dd " << dd - // << " coord " << dcoord[ii*3+dd] - // << " nn_f " << nn_force[ii*3+dd] - // << " ww " << weight[ii] - // << endl; - // } - } - // cout << "nn " << dcoord[idx*3] << " " << weight[idx] << endl; - } - - // double protect_level = 1e-3; - // cout << "with protect_level " << protect_level << endl; - for (int ii = 0; ii < nall; ++ii) { - for (int dd = 0; dd < 3; ++dd) { - double pref = (1 - weight[ii]); - if (fabs(pref) < protect_level) { - pref = protect_level; - } - of[ii * 3 + dd] += pref * ff_bd_force[ii * 3 + dd]; - // if (fabs(ff_bd_force[ii*3+dd]) > 1e6) { - // cout << " ii " << ii - // << " dd " << dd - // << " coord " << dcoord[ii*3+dd] - // << " ff_f " << ff_bd_force[ii*3+dd] - // << " ww " << 1 - weight[ii] - // << endl; - // } - } - // cout << "ff " << dcoord[ii*3] << " " << 1-weight[ii] << endl; - } - for (int ii = 0; ii < nall; ++ii) { - for (int dd = 0; dd < 3; ++dd) { - of[ii * 3 + dd] += (1 - weight[ii]) * ff_nb_force[ii * 3 + dd]; - // if (fabs(ff_nb_force[ii*3+dd]) > 1e6) { - // cout << " ii " << ii - // << " dd " << dd - // << " coord " << dcoord[ii*3+dd] - // << " ff_f " << ff_nb_force[ii*3+dd] - // << " ww " << 1 - weight[ii] - // << endl; - // } - } - // cout << "ff " << dcoord[ii*3] << " " << 1-weight[ii] << endl; - } - - for (int ii = 0; ii < of.size() / 3; ++ii) { - for (int dd = 0; dd < 3; ++dd) { - ff_sum[dd] += ff_bd_force[ii * 3 + dd]; - } - } - // cout << ff_sum[0] << " " << ff_sum[1] << " " << ff_sum[2] << " " - // <& box, - const VALUETYPE& rnn_, - const VALUETYPE& rhy_, - const VALUETYPE& rc_, - const VALUETYPE& protect_level_) - : AdWeight(protect_level_) { - assert(box.size() == 9); - center.resize(3); - for (int ii = 0; ii < 3; ++ii) { - center[ii] = 0.5 * box[3 * ii + ii]; - } - rnn = rnn_; - rhy = rhy_; - rc = rc_; -} - -void SlabWeight::zone_tag(vector& tag, - const vector& coord) const { - int natoms = coord.size() / 3; - tag.resize(natoms, 0); - - // slab axis x - VALUETYPE radius = rnn + rhy; - for (int ii = 0; ii < natoms; ++ii) { - VALUETYPE posi = fabs(coord[ii * 3] - center[0]); - if (posi < radius) { - tag[ii] = 3; - } else if (posi < radius + rc) { - tag[ii] = 2; - } else if (posi < radius + rc * 2) { - tag[ii] = 1; - } else { - tag[ii] = 0; - } - } -} - -// dirty hacking -void SlabWeight::atom_weight(vector& weight, - vector& weight_x, - const vector& coord) const { - CosSwitch cs(rnn, rnn + rhy); - - int natoms = coord.size() / 3; - weight.resize(natoms, 0); - weight_x.resize(natoms, 0); - // slab axis x - // for (int ii = 0; ii < natoms; ++ii){ - // VALUETYPE posi = fabs(coord[ii*3] - center[0]); - // cs.eval (weight[ii], posi); - // // if (posi < radius){ - // // weight[ii] = 1.; - // // } - // // else { - // // weight[ii] = 0.; - // // } - // } - // for (int ii = 0; ii < natoms/3; ++ii){ - // VALUETYPE posi = fabs(coord[ii*3] - center[0]); - // cs.eval (weight[ii], posi); - // weight[natoms/3 + ii*2 + 0] = weight[ii]; - // weight[natoms/3 + ii*2 + 1] = weight[ii]; - // // weight_x - // weight_x[ii] = posi; - // weight_x[natoms/3 + ii*2 + 0] = posi; - // weight_x[natoms/3 + ii*2 + 1] = posi; - // // if (posi < radius){ - // // weight[ii] = 1.; - // // } - // // else { - // // weight[ii] = 0.; - // // } - // } - for (int ii = 0; ii < natoms; ii += 3) { - VALUETYPE posi = fabs(coord[ii * 3] - center[0]); - cs.eval(weight[ii], posi); - weight[ii + 1] = weight[ii]; - weight[ii + 2] = weight[ii]; - // weight_x - weight_x[ii] = posi; - weight_x[ii + 1] = posi; - weight_x[ii + 2] = posi; - } -} diff --git a/source/md/src/Convert.cc b/source/md/src/Convert.cc deleted file mode 100644 index b8014bf974..0000000000 --- a/source/md/src/Convert.cc +++ /dev/null @@ -1,115 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "Convert.h" - -#include -#include -#include - -template -Convert::Convert(const vector& atomname, - map& name_type_map, - map& name_mass_map, - map& name_charge_map, - const bool sort_) { - int natoms = atomname.size(); - atype.resize(natoms); - amass.resize(natoms); - vector tmp_charge(natoms); - for (unsigned ii = 0; ii < atype.size(); ++ii) { - atype[ii] = name_type_map[atomname[ii]]; - amass[ii] = name_mass_map[atomname[ii]]; - tmp_charge[ii] = name_charge_map[atomname[ii]]; - } - vector > > sorting(natoms); - for (unsigned ii = 0; ii < sorting.size(); ++ii) { - sorting[ii] = pair >( - atype[ii], pair(ii, amass[ii])); - } - if (sort_) { - sort(sorting.begin(), sorting.end()); - } - idx_map_nnp2gro.resize(natoms); - idx_map_gro2nnp.resize(natoms); - for (unsigned ii = 0; ii < idx_map_nnp2gro.size(); ++ii) { - idx_map_nnp2gro[ii] = sorting[ii].second.first; - idx_map_gro2nnp[sorting[ii].second.first] = ii; - atype[ii] = sorting[ii].first; - amass[ii] = sorting[ii].second.second; - } - acharge.resize(natoms); - for (int ii = 0; ii < natoms; ++ii) { - int gro_i = idx_map_nnp2gro[ii]; - acharge[ii] = tmp_charge[gro_i]; - } -} - -template -void Convert::gro2nnp(vector& coord, - vector& veloc, - vector& box, - const vector >& posi, - const vector >& velo, - const vector& box_size) const { - assert(posi.size() == idx_map_nnp2gro.size()); - assert(velo.size() == idx_map_nnp2gro.size()); - int natoms = idx_map_nnp2gro.size(); - coord.resize(3 * static_cast(natoms)); - veloc.resize(3 * static_cast(natoms)); - for (unsigned ii = 0; ii < natoms; ++ii) { - int gro_i = idx_map_nnp2gro[ii]; - for (int dd = 0; dd < 3; ++dd) { - coord[ii * 3 + dd] = posi[gro_i][dd] * 10; - veloc[ii * 3 + dd] = velo[gro_i][dd] * 10; - } - } - box.resize(9); - for (int dd = 0; dd < 9; ++dd) { - box[dd] = box_size[dd] * 10; - } -} - -template -void Convert::nnp2gro(vector >& posi, - vector >& velo, - vector& box_size, - const vector& coord, - const vector& veloc, - const vector& box) const { - int natoms = idx_map_nnp2gro.size(); - posi.resize(natoms); - velo.resize(natoms); - for (unsigned ii = 0; ii < posi.size(); ++ii) { - posi[ii].resize(3); - velo[ii].resize(3); - } - for (unsigned ii = 0; ii < posi.size(); ++ii) { - int gro_i = idx_map_nnp2gro[ii]; - for (int dd = 0; dd < 3; ++dd) { - posi[gro_i][dd] = coord[ii * 3 + dd] * 0.1; - velo[gro_i][dd] = veloc[ii * 3 + dd] * 0.1; - } - } - box_size.resize(9); - for (int dd = 0; dd < 9; ++dd) { - box_size[dd] = box[dd] * 0.1; - } -} - -template -void Convert::idx_gro2nnp(vector& out, - const vector& in) const { - for (unsigned ii = 0; ii < in.size(); ++ii) { - out[ii] = idx_map_gro2nnp[in[ii]]; - } -} - -template -void Convert::idx_nnp2gro(vector& out, - const vector& in) const { - for (unsigned ii = 0; ii < in.size(); ++ii) { - out[ii] = idx_map_nnp2gro[in[ii]]; - } -} - -template class Convert; -template class Convert; diff --git a/source/md/src/Gaussian.cc b/source/md/src/Gaussian.cc deleted file mode 100644 index 80d78c9385..0000000000 --- a/source/md/src/Gaussian.cc +++ /dev/null @@ -1,20 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "Gaussian.h" - -void Gaussian::set_seed(unsigned long s) { - RandomGenerator_MT19937::init_genrand(s); -} - -void Gaussian::gen(double* vec, const int numb_gen) { - const double epsilon = std::numeric_limits::min(); - const double two_pi = 2.0 * M_PI; - - for (int ii = 0; ii < numb_gen; ++ii) { - double u0, u1; - do { - u0 = RandomGenerator_MT19937::genrand_real3(); - u1 = RandomGenerator_MT19937::genrand_real3(); - } while (u0 < epsilon); - vec[ii] = sqrt(-2.0 * log(u0)) * cos(two_pi * u1); - } -} diff --git a/source/md/src/GroFileManager.cc b/source/md/src/GroFileManager.cc deleted file mode 100644 index d61fbb7b97..0000000000 --- a/source/md/src/GroFileManager.cc +++ /dev/null @@ -1,283 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "GroFileManager.h" - -#include -#include -#include -#include - -using namespace std; - -class WrongFileFormat {}; - -bool GroFileManager::readTop(const std::string &filename, - std::vector &molnames, - std::vector &nmols) { - molnames.clear(); - nmols.clear(); - - std::ifstream in(filename.c_str()); - if (in.bad()) { - std::cerr << "cannot open file " << filename << std::endl; - return false; - } - char line[1024]; - std::string target("[ molecules ]"); - bool find = false; - while (!in.eof()) { - in.getline(line, 1024, '\n'); - if (target == std::string(line)) { - find = true; - break; - } - } - if (!find) { - std::cerr << "cannot find [ molecules ] in file " << filename - << ". please check there is no space after \"]\"\n"; - return false; - } - - // while (!(in.getline (line, 1024, '\n')).eof()){ - // // if (line[0] == '['){ - // // break; - // // } - // // char name[1024]; - // // int number; - // // sscanf (line, "%s%d", name, &number); - // // molnames.push_back (std::string(name)); - // // nmols.push_back (number); - // } - - std::string name; - int number; - while (!(in >> name >> number).eof()) { - if (name[0] == '[') { - break; - } - if (name.empty()) { - break; - } - // std::cout << name << std::endl; - molnames.push_back(name); - nmols.push_back(number); - } - - return true; -} - -template -bool GroFileManager::writePotenFile(const double &rmin, - const double &rcut, - const double &interval, - UnitaryFunction1 &f, - UnitaryFunction2 &fp, - UnitaryFunction3 &g, - UnitaryFunction4 &gp, - UnitaryFunction5 &h, - UnitaryFunction6 &hp, - const std::string &filename) { - FILE *filep = fopen(filename.c_str(), "w"); - if (filep == NULL) { - std::cerr << "cannot open file " << filename << std::endl; - return false; - } - - double upper = rcut + 1; - double nx; - if (int(upper / interval) != upper / interval) { - nx = int(upper / interval) + 1; - } else { - nx = int(upper / interval); - } - upper = interval * nx; - - int i = 0; - for (i = 0; i <= nx + 1; ++i) { - double x = i * interval; - if (x < rmin) { - fprintf(filep, "%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\n", x, 0., - 0., 0., 0., 0., 0.); - } else { - fprintf(filep, "%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\n", x, - f(x), -fp(x), g(x), -gp(x), h(x), -hp(x)); - } - } - - fclose(filep); - return true; -} - -void GroFileManager::read(const std::string &name, - std::vector &resdindex, - std::vector &resdname, - std::vector &atomname, - std::vector &atomindex, - std::vector > &posi, - std::vector > &velo, - std::vector &boxsize_) { - FILE *fp = fopen(name.c_str(), "r"); - if (fp == NULL) { - std::cerr << "cannot open file " << name << std::endl; - return; - } - while (fgetc(fp) != '\n'); - int npart; - fscanf(fp, "%d\n", &npart); - fclose(fp); - - resdindex.clear(); - resdname.clear(); - atomname.clear(); - atomindex.clear(); - posi.clear(); - velo.clear(); - vector boxsize; - boxsize.resize(3); - - fp = fopen(name.c_str(), "r"); - while (fgetc(fp) != '\n'); - while (fgetc(fp) != '\n'); - char line[1024]; - for (int i = 0; i < npart; ++i) { - fgets(line, 1024, fp); - char tmp[1024]; - int tmpd; - char tmps[1024]; - for (unsigned j = 0; j < 5; ++j) { - tmp[j] = line[j]; - } - tmp[5] = '\0'; - if (sscanf(tmp, "%d", &tmpd) != 1) { - throw WrongFileFormat(); - } - resdindex.push_back(tmpd); - - for (unsigned j = 0; j < 5; ++j) { - tmp[j] = line[j + 5]; - } - tmp[5] = '\0'; - if (sscanf(tmp, "%s", tmps) != 1) { - throw WrongFileFormat(); - } - resdname.push_back(tmps); - - for (unsigned j = 0; j < 5; ++j) { - tmp[j] = line[j + 10]; - } - tmp[5] = '\0'; - if (sscanf(tmp, "%s", tmps) != 1) { - throw WrongFileFormat(); - } - atomname.push_back(tmps); - - for (unsigned j = 0; j < 5; ++j) { - tmp[j] = line[j + 15]; - } - tmp[5] = '\0'; - if (sscanf(tmp, "%d", &tmpd) != 1) { - throw WrongFileFormat(); - } - atomindex.push_back(tmpd); - - double a, b, c; - double d, e, f; - std::vector tmpp(3); - std::vector tmpv(3); - - int tag = sscanf(&line[20], "%lf%lf%lf%lf%lf%lf", &a, &b, &c, &d, &e, &f); - tmpp[0] = a; - tmpp[1] = b; - tmpp[2] = c; - switch (tag) { - case 6: - tmpv[0] = d; - tmpv[1] = e; - tmpv[2] = f; - break; - case 3: - tmpv[0] = 0.f; - tmpv[1] = 0.f; - tmpv[2] = 0.f; - break; - default: - throw WrongFileFormat(); - } - - posi.push_back(tmpp); - velo.push_back(tmpv); - } - int tag = 0; - double rbox[9]; - tag = fscanf(fp, "%lf%lf%lf%lf%lf%lf%lf%lf%lf", rbox + 0, rbox + 1, rbox + 2, - rbox + 3, rbox + 4, rbox + 5, rbox + 6, rbox + 7, rbox + 8); - fclose(fp); - - boxsize_.resize(9, 0.); - fill(boxsize_.begin(), boxsize_.end(), 0.); - - if (tag == 9) { - boxsize_[0] = rbox[0]; - boxsize_[4] = rbox[1]; - boxsize_[8] = rbox[2]; - boxsize_[0 * 3 + 1] = rbox[3]; - boxsize_[0 * 3 + 2] = rbox[4]; - boxsize_[1 * 3 + 0] = rbox[5]; - boxsize_[1 * 3 + 2] = rbox[6]; - boxsize_[2 * 3 + 0] = rbox[7]; - boxsize_[2 * 3 + 1] = rbox[8]; - } else { - assert(tag == 3); - boxsize_[0] = rbox[0]; - boxsize_[4] = rbox[1]; - boxsize_[8] = rbox[2]; - } -} - -void GroFileManager::write(const std::string &name, - const std::vector &resdindex, - const std::vector &resdname, - const std::vector &atomname, - const std::vector &atomindex, - const std::vector > &posi, - const std::vector > &velo, - const std::vector &boxsize) { - FILE *fp = fopen(name.c_str(), "w"); - if (fp == NULL) { - std::cerr << "cannot open file " << name << std::endl; - return; - } - // std::copy (atomname.begin(), atomname.end(), - // std::ostream_iterator(std::cout, "\n")); - - fprintf(fp, "\n%d\n", int(resdindex.size())); - for (int i = 0; i < int(resdindex.size()); ++i) { - fprintf(fp, "%5d%5s%5s%5d%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n", - resdindex[i] % 100000, (char *)(resdname[i].c_str()), - (char *)(atomname[i].c_str()), atomindex[i] % 100000, posi[i][0], - posi[i][1], posi[i][2], velo[i][0], velo[i][1], velo[i][2]); - } - // vector box(3); - // for (int ii = 0; ii < 3; ++ii) box[ii] = boxsize[3*ii+ii]; - if (boxsize.size() == 3) { - fprintf(fp, "%f %f %f\n", boxsize[0], boxsize[1], boxsize[2]); - } else if (boxsize.size() == 9) { - fprintf(fp, "%f %f %f %f %f %f %f %f %f \n", boxsize[0 * 3 + 0], - boxsize[1 * 3 + 1], boxsize[2 * 3 + 2], boxsize[0 * 3 + 1], - boxsize[0 * 3 + 2], boxsize[1 * 3 + 0], boxsize[1 * 3 + 2], - boxsize[2 * 3 + 0], boxsize[2 * 3 + 1]); - } - - fclose(fp); -} - -struct F { - double operator()(double x) { return 1. / x; } -}; -struct Zero { - double operator()(double x) { return 0; } -}; diff --git a/source/md/src/HarmonicAngle.cc b/source/md/src/HarmonicAngle.cc deleted file mode 100644 index 71915683ab..0000000000 --- a/source/md/src/HarmonicAngle.cc +++ /dev/null @@ -1,99 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "HarmonicAngle.h" - -#include -#include - -#include "common.h" -#include "mymath.h" - -HarmonicAngle::HarmonicAngle(const VALUETYPE& ka_, const VALUETYPE& tt_) - : ka(ka_), tt(tt_) {} - -inline bool compute_variable(const VALUETYPE* rij, - const VALUETYPE* rkj, - VALUETYPE* var, - VALUETYPE* dvardcos, - VALUETYPE* cos_theta) { - *cos_theta = cos(rij[0], rij[1], rij[2], rkj[0], rkj[1], rkj[2]); - *var = acos(*cos_theta); - - VALUETYPE cos_theta2 = *cos_theta * *cos_theta; - if (cos_theta2 >= 1) { - *dvardcos = 1.; - return false; - } - *dvardcos = -1. / sqrt(1. - cos_theta2); - return true; -} - -void HarmonicAngle::compute(VALUETYPE& ener, - vector& force, - vector& virial, - const vector& coord, - const vector& atype, - const SimulationRegion& region, - const vector& alist) { - // all set zeros - for (unsigned _ = 0; _ < alist.size(); _ += 3) { - int ii = alist[_]; - int jj = alist[_ + 1]; - int kk = alist[_ + 2]; - - VALUETYPE rij[3], rkj[3]; - region.diffNearestNeighbor(&coord[ii * 3], &coord[jj * 3], rij); - region.diffNearestNeighbor(&coord[kk * 3], &coord[jj * 3], rkj); - - VALUETYPE var(0), dvardcos(0), cos_theta(0); - bool apply_force = compute_variable(rij, rkj, &var, &dvardcos, &cos_theta); - - VALUETYPE dudvar(0), angle_energy(0); - VALUETYPE diff = var - tt; - VALUETYPE pdiff = ka * diff; - dudvar = -pdiff; - angle_energy = VALUETYPE(0.5) * pdiff * diff; - - ener += angle_energy; - - // VALUETYPE fijx, fijy, fijz; - // VALUETYPE fkjx, fkjy, fkjz; - VALUETYPE fij[3]; - VALUETYPE fkj[3]; - - if (apply_force) { - VALUETYPE dudcos = dudvar * dvardcos; - VALUETYPE rij2 = dot(rij, rij); - VALUETYPE rkj2 = dot(rkj, rkj); - VALUETYPE invrij = 1. / sqrt(rij2); - VALUETYPE invrkj = 1. / sqrt(rkj2); - VALUETYPE invrij2 = invrij * invrij; - VALUETYPE invrkj2 = invrkj * invrkj; - VALUETYPE invrijrkj = invrij * invrkj; - // can be further optimized: - fij[0] = dudcos * (rkj[0] * invrijrkj - rij[0] * invrij2 * cos_theta); - fij[1] = dudcos * (rkj[1] * invrijrkj - rij[1] * invrij2 * cos_theta); - fij[2] = dudcos * (rkj[2] * invrijrkj - rij[2] * invrij2 * cos_theta); - fkj[0] = dudcos * (rij[0] * invrijrkj - rkj[0] * invrkj2 * cos_theta); - fkj[1] = dudcos * (rij[1] * invrijrkj - rkj[1] * invrkj2 * cos_theta); - fkj[2] = dudcos * (rij[2] * invrijrkj - rkj[2] * invrkj2 * cos_theta); - } else { - fij[0] = fij[1] = fij[2] = fkj[0] = fkj[1] = fkj[2] = VALUETYPE(0); - } - - force[3 * ii + 0] += fij[0]; - force[3 * ii + 1] += fij[1]; - force[3 * ii + 2] += fij[2]; - force[3 * kk + 0] += fkj[0]; - force[3 * kk + 1] += fkj[1]; - force[3 * kk + 2] += fkj[2]; - force[3 * jj + 0] -= fij[0] + fkj[0]; - force[3 * jj + 1] -= fij[1] + fkj[1]; - force[3 * jj + 2] -= fij[2] + fkj[2]; - for (int dd0 = 0; dd0 < 3; ++dd0) { - for (int dd1 = 0; dd1 < 3; ++dd1) { - virial[dd0 * 3 + dd1] -= 0.5 * fij[dd0] * rij[dd1]; - virial[dd0 * 3 + dd1] -= 0.5 * fkj[dd0] * rkj[dd1]; - } - } - } -} diff --git a/source/md/src/HarmonicBond.cc b/source/md/src/HarmonicBond.cc deleted file mode 100644 index cadbf46b53..0000000000 --- a/source/md/src/HarmonicBond.cc +++ /dev/null @@ -1,50 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "HarmonicBond.h" - -#include -#include - -#include "common.h" - -HarmonicBond::HarmonicBond(const VALUETYPE& kk_, const VALUETYPE& bb_) - : kk(kk_), bb(bb_) {} - -void HarmonicBond::hb_inner(VALUETYPE& ae, VALUETYPE& af, const VALUETYPE& r1) { - VALUETYPE diff = r1 - bb; - // cout << bb << " " << r1 << endl; - VALUETYPE pdiff = kk * diff; - af = -pdiff / r1; - ae = 0.5 * pdiff * diff; -} - -void HarmonicBond::compute(VALUETYPE& ener, - vector& force, - vector& virial, - const vector& coord, - const vector& atype, - const SimulationRegion& region, - const vector& blist) { - // all set zeros - for (unsigned _ = 0; _ < blist.size(); _ += 2) { - int ii = blist[_]; - int jj = blist[_ + 1]; - VALUETYPE diff[3]; - region.diffNearestNeighbor(&coord[ii * 3], &coord[jj * 3], diff); - VALUETYPE r2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]; - VALUETYPE r1 = sqrt(r2); - VALUETYPE ae, af; - hb_inner(ae, af, r1); - for (int dd = 0; dd < 3; ++dd) { - force[ii * 3 + dd] += af * diff[dd]; - } - for (int dd = 0; dd < 3; ++dd) { - force[jj * 3 + dd] -= af * diff[dd]; - } - ener += ae; - for (int dd0 = 0; dd0 < 3; ++dd0) { - for (int dd1 = 0; dd1 < 3; ++dd1) { - virial[dd0 * 3 + dd1] -= 0.5 * diff[dd0] * af * diff[dd1]; - } - } - } -} diff --git a/source/md/src/Integrator.cc b/source/md/src/Integrator.cc deleted file mode 100644 index 333776c7f2..0000000000 --- a/source/md/src/Integrator.cc +++ /dev/null @@ -1,94 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "Integrator.h" - -#include - -template -void Integrator::stepVeloc(vector &vv, - const vector &ff, - const vector &mass, - const double &dt, - const vector &freez) const { - int natoms = ff.size() / 3; - for (int kk = 0; kk < natoms; ++kk) { - VALUETYPE invmdt = dt / (mass[kk] * massConst); - vv[kk * 3 + 0] += ff[kk * 3 + 0] * invmdt; - vv[kk * 3 + 1] += ff[kk * 3 + 1] * invmdt; - vv[kk * 3 + 2] += ff[kk * 3 + 2] * invmdt; - } - for (unsigned ii = 0; ii < freez.size(); ++ii) { - int kk = freez[ii]; - vv[kk * 3 + 0] = 0; - vv[kk * 3 + 1] = 0; - vv[kk * 3 + 2] = 0; - } -} - -template -void Integrator::stepCoord(vector &rr, - const vector &vv, - const double &dt) const { - for (unsigned kk = 0; kk < vv.size(); ++kk) { - rr[kk] += dt * vv[kk]; - } -} - -template -ThermostatLangevin::ThermostatLangevin(const VALUETYPE T_, - const VALUETYPE tau_, - const long long int seed) { - reinit(T_, tau_, seed); -} - -template -void ThermostatLangevin::reinit(const VALUETYPE T_, - const VALUETYPE tau_, - const long long int seed) { - gaussian.set_seed(seed); - temperature = T_; - kT = UnitManager::BoltzmannConstant * T_; - gamma = 1. / tau_; - VALUETYPE twogammakT = 2. * gamma * kT; - sigma = 1. / sqrt(twogammakT) * twogammakT; - sigmainvsqrt2gamma = VALUETYPE(sigma / sqrt(2. * gamma)); -} - -template -void ThermostatLangevin::stepOU(vector &vv, - const vector &mass, - const double &dt, - const vector &freez) const { - VALUETYPE emgammat = exp(-gamma * dt); - VALUETYPE sqrt1memgammat2 = sqrt(1. - emgammat * emgammat); - VALUETYPE prefR = sigmainvsqrt2gamma * sqrt1memgammat2; - - int numb_part = mass.size(); - assert(int(vv.size()) == 3 * numb_part); - - double *all_rands = (double *)malloc(sizeof(double) * numb_part * 3); - gaussian.gen(all_rands, numb_part * 3); - - for (int kk = 0; kk < numb_part; ++kk) { - VALUETYPE sm = mass[kk] * UnitManager::IntegratorMassConstant; - VALUETYPE invsqrtm = 1. / sqrt(sm); - vv[kk * 3 + 0] = - emgammat * vv[kk * 3 + 0] + prefR * invsqrtm * all_rands[kk * 3 + 0]; - vv[kk * 3 + 1] = - emgammat * vv[kk * 3 + 1] + prefR * invsqrtm * all_rands[kk * 3 + 1]; - vv[kk * 3 + 2] = - emgammat * vv[kk * 3 + 2] + prefR * invsqrtm * all_rands[kk * 3 + 2]; - } - for (unsigned ii = 0; ii < freez.size(); ++ii) { - int kk = freez[ii]; - vv[kk * 3 + 0] = 0; - vv[kk * 3 + 1] = 0; - vv[kk * 3 + 2] = 0; - } - - free(all_rands); -} - -template class Integrator; -template class Integrator; -template class ThermostatLangevin; -template class ThermostatLangevin; diff --git a/source/md/src/Interpolation.cpp b/source/md/src/Interpolation.cpp deleted file mode 100644 index 0e5e7a9214..0000000000 --- a/source/md/src/Interpolation.cpp +++ /dev/null @@ -1,553 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "Interpolation.h" - -#include - -void Interpolation::piece6OrderInterpol(const double& a, - const double& b, - const double& va, - const double& vb, - const double& da, - const double& db, - const double& dda, - const double& ddb, - Poly& p) { - std::vector standardPolys(6); - for (unsigned i = 0; i < 6; ++i) { - standardPolys[i].getOrder() = 5; - standardPolys[i].getCoeffs().resize(6); - } - standardPolys[0].getCoeffs()[0] = 1; - standardPolys[0].getCoeffs()[1] = 0; - standardPolys[0].getCoeffs()[2] = 0; - standardPolys[0].getCoeffs()[3] = -10; - standardPolys[0].getCoeffs()[4] = 15; - standardPolys[0].getCoeffs()[5] = -6; - - standardPolys[1].getCoeffs()[0] = 0; - standardPolys[1].getCoeffs()[1] = 0; - standardPolys[1].getCoeffs()[2] = 0; - standardPolys[1].getCoeffs()[3] = 10; - standardPolys[1].getCoeffs()[4] = -15; - standardPolys[1].getCoeffs()[5] = 6; - - standardPolys[2].getCoeffs()[0] = 0; - standardPolys[2].getCoeffs()[1] = 1; - standardPolys[2].getCoeffs()[2] = 0; - standardPolys[2].getCoeffs()[3] = -6; - standardPolys[2].getCoeffs()[4] = 8; - standardPolys[2].getCoeffs()[5] = -3; - - standardPolys[3].getCoeffs()[0] = 0; - standardPolys[3].getCoeffs()[1] = 0; - standardPolys[3].getCoeffs()[2] = 0; - standardPolys[3].getCoeffs()[3] = -4; - standardPolys[3].getCoeffs()[4] = 7; - standardPolys[3].getCoeffs()[5] = -3; - - standardPolys[4].getCoeffs()[0] = 0; - standardPolys[4].getCoeffs()[1] = 0; - standardPolys[4].getCoeffs()[2] = 0.5; - standardPolys[4].getCoeffs()[3] = -1.5; - standardPolys[4].getCoeffs()[4] = 1.5; - standardPolys[4].getCoeffs()[5] = -0.5; - - standardPolys[5].getCoeffs()[0] = 0; - standardPolys[5].getCoeffs()[1] = 0; - standardPolys[5].getCoeffs()[2] = 0; - standardPolys[5].getCoeffs()[3] = 0.5; - standardPolys[5].getCoeffs()[4] = -1; - standardPolys[5].getCoeffs()[5] = 0.5; - - std::vector scaledPolys(6); - double tmpa(1. / (b - a)); - double tmpb(-a / (b - a)); - for (unsigned i = 0; i < 6; ++i) { - standardPolys[i].valueLinearPoly(tmpa, tmpb, scaledPolys[i]); - } - scaledPolys[2] *= 1. / tmpa; - scaledPolys[3] *= 1. / tmpa; - scaledPolys[4] *= 1. / tmpa / tmpa; - scaledPolys[5] *= 1. / tmpa / tmpa; - - p.zero(); - p += (scaledPolys[0] *= va); - p += (scaledPolys[1] *= vb); - p += (scaledPolys[2] *= da); - p += (scaledPolys[3] *= db); - p += (scaledPolys[4] *= dda); - p += (scaledPolys[5] *= ddb); - - return; -} - -void Interpolation::pieceLinearInterpol(const double& a, - const double& b, - const double& va, - const double& vb, - Poly& p) { - double k = (vb - va) / (b - a); - std::vector tmp(2, 0); - tmp[0] += va; - tmp[0] += k * (-a); - tmp[1] = k; - p.reinit(tmp); -} - -void Interpolation::piecewiseLinear(const std::vector& x, - const std::vector& y, - PiecewisePoly& ps) { - std::vector::const_iterator pxp1 = x.begin(); - std::vector::const_iterator px = (pxp1++); - std::vector::const_iterator pyp1 = y.begin(); - std::vector::const_iterator py = (pyp1++); - ps.clear(); - Poly tmpp; - for (; pxp1 != x.end(); ++pxp1, ++pyp1, ++px, ++py) { - pieceLinearInterpol(*px, *pxp1, *py, *pyp1, tmpp); - ps.get_x().push_back(*px); - ps.get_p().push_back(tmpp); - } - ps.get_x().push_back(*px); -} - -void Interpolation::pieceSecondDerivativeInterpol(const double& a, - const double& b, - const double& va, - const double& vb, - const double& dda, - const double& ddb, - Poly& p) { - std::vector tmp(2, 0); - double k = (vb - va) / (b - a); - tmp[0] += va; - tmp[0] += k * (-a); - tmp[1] = k; - p.reinit(tmp); - - tmp[1] = 1; - tmp[0] = -a; - Poly l1(tmp); - tmp[0] = -b; - Poly l2(tmp); - l1 *= l2; - - tmp[1] = 1. / 6. / (a - b); - tmp[0] = 1. / 6. * (a - 2 * b) / (a - b); - Poly p1(tmp); - p1 *= l1; - p1 *= dda; - - tmp[1] *= -1; - tmp[0] = 1. / 6. * (b - 2 * a) / (b - a); - Poly p2(tmp); - p2 *= l1; - p2 *= ddb; - - p += p1; - p += p2; -} - -void Interpolation::secondDerivativeInterpol( - const std::vector::const_iterator& xbegin, - const std::vector::const_iterator& xend, - const std::vector::const_iterator& vbegin, - const std::vector::const_iterator& ddbegin, - PiecewisePoly& ps) { - ps.clear(); - std::vector::const_iterator xb(xbegin), vb(vbegin), ddb(ddbegin); - std::vector::const_iterator xp(xbegin), vp(vbegin), ddp(ddbegin); - ++xp, ++vp, ++ddp; - while (xp != xend) { - ps.get_x().push_back(*xb); - Poly tmpp; - pieceSecondDerivativeInterpol(*(xb++), *(xp++), *(vb++), *(vp++), *(ddb++), - *(ddp++), tmpp); - ps.get_p().push_back(tmpp); - } - ps.get_x().push_back(*xb); -} - -void Interpolation::pieceHermiteInterpol(const double& a, - const double& b, - const double& va, - const double& vb, - const double& da, - const double& db, - Poly& p) { - std::vector tmp(2, 0); - Poly t; - tmp[0] = (-2 * a / (b - a) + 1); - tmp[1] = (2 / (b - a)); - Poly a0(tmp); - tmp[0] = -b / (a - b); - tmp[1] = 1 / (a - b); - t.reinit(tmp); - a0 *= t; - a0 *= t; - tmp[0] = -2 * b / (a - b) + 1; - tmp[1] = 2 / (a - b); - Poly a1(tmp); - tmp[0] = -a / (b - a); - tmp[1] = 1 / (b - a); - t.reinit(tmp); - a1 *= t; - a1 *= t; - - tmp[0] = -a; - tmp[1] = 1; - Poly b0(tmp); - tmp[0] = -b / (a - b); - tmp[1] = 1 / (a - b); - t.reinit(tmp); - b0 *= t; - b0 *= t; - tmp[0] = -b; - tmp[1] = 1; - Poly b1(tmp); - tmp[0] = -a / (b - a); - tmp[1] = 1 / (b - a); - t.reinit(tmp); - b1 *= t; - b1 *= t; - - p.zero(); - a0 *= va; - a1 *= vb; - b0 *= da; - b1 *= db; - p += a0; - p += a1; - p += b0; - p += b1; -} - -// lbegin--lend, stores lambda -// ubegin--uend, stores mu -bool Interpolation::solverForSplinePeriodic( - const std::vector::const_iterator& lbegin, - const std::vector::const_iterator& lend, - const std::vector::iterator& ubegin, - const std::vector::iterator& uend) { - std::vector la, lb, lc, ld; - for (std::vector::const_iterator i = lbegin; i != lend; ++i) { - la.push_back(1 - *i); - lb.push_back(2); - lc.push_back(*i); - ld.push_back(0); - } - // ld.front() = 1 - *lbegin; - ld[0] = 1 - lc[0]; - int num = ld.size(); - ld[num - 2] = lc[num - 2]; - ld[num - 1] = lb[num - 1]; - - std::vector::iterator pu = ubegin; - std::vector::iterator pu_1 = pu++; - for (int i = 1; i < num - 1; ++i, ++pu, ++pu_1) { - if (lb[i - 1] == 0) { - return false; - } - double ratio = -la[i] / lb[i - 1]; - lb[i] += lc[i - 1] * ratio; - ld[i] += ld[i - 1] * ratio; - *pu += *pu_1 * ratio; - } - int i = num - 1; - if (lb[i - 1] == 0) { - return false; - } - double ratio = -la[i] / lb[i - 1]; - lb[i] += ld[i - 1] * ratio; - ld[i] = lb[i]; - *pu += *pu_1 * ratio; - - // std::cout << lc.back() << std::endl; - // std::cout << lc.front() << std::endl; - ratio = -lb[0] / lc.back(); - ld[0] += ratio * ld[num - 1]; - *ubegin += ratio * *pu; - lb[0] = 0; - - // std::cout << ld.size() << std::endl; - ld.insert(ld.begin(), ld.back()); - // std::cout << ld.size() << std::endl; - ld.pop_back(); - // std::cout << ld.size() << std::endl; - double before = 0.; - // std::cout << "##############################" << std::endl; - // std::copy(ubegin, uend, std::ostream_iterator(std::cout, "\n")); - // std::cout << "##############################" << std::endl; - for (std::vector::iterator tmpu = ubegin; tmpu != uend; ++tmpu) { - if (tmpu == ubegin) { - before = *tmpu; - *tmpu = *pu; - } else { - double beforetmp = *tmpu; - *tmpu = before; - before = beforetmp; - } - } - // std::copy(ubegin, uend, std::ostream_iterator(std::cout, "\n")); - // std::cout << "##############################" << std::endl; - lc.insert(lc.begin(), *lbegin); - lc.pop_back(); - lc.back() = ld.back(); - lb.insert(lb.begin(), 0.); - lb.pop_back(); - - pu = ubegin; - pu++; - pu_1 = pu++; - for (int i = 2; i < num - 1; ++i, ++pu, ++pu_1) { - if (lc[i - 1] == 0) { - return false; - } - double ratio = -lb[i] / lc[i - 1]; - ld[i] += ld[i - 1] * ratio; - *pu += *pu_1 * ratio; - } - i = num - 1; - if (lc[i - 1] == 0) { - return false; - } - ratio = -lb[i] / lc[i - 1]; - lc[i] += ld[i - 1] * ratio; - ld[i] = lc[i]; - *pu += *pu_1 * ratio; - - *pu /= lc[num - 1]; - for (int i = num - 2; i >= 0; --i, --pu_1) { - *pu_1 = (*pu_1 - *pu * ld[i]) / lc[i]; - } - - return true; -} - -bool Interpolation::splinePeriodic(const std::vector& x, - const std::vector& y, - PiecewisePoly& ps) { - std::vector lambda(x.size() - 1); - std::vector mu(x.size() - 1); - std::vector dx; - - std::vector::const_iterator i = x.begin(); - std::vector::const_iterator j = i; - for (++j; j != x.end(); ++i, ++j) { - dx.push_back(*j - *i); - } - lambda[0] = dx.back() / (dx.back() + dx.front()); - mu[0] = 3 * ((1 - lambda.front()) / dx.back() * (y[0] - y[y.size() - 2]) + - lambda.front() / dx.front() * (y[1] - y[0])); - for (unsigned i = 1; i < lambda.size(); ++i) { - lambda[i] = dx[i - 1] / (dx[i - 1] + dx[i]); - mu[i] = 3 * ((1 - lambda[i]) / dx[i - 1] * (y[i] - y[i - 1]) + - lambda[i] / dx[i] * (y[i + 1] - y[i])); - } - - bool tag = solverForSplinePeriodic(lambda.begin(), lambda.end(), mu.begin(), - mu.end()); - if (!tag) { - return false; - } - - ps.get_x() = x; - ps.get_p().clear(); - for (unsigned i = 0; i < x.size() - 2; ++i) { - Poly tmpp; - pieceHermiteInterpol(x[i], x[i + 1], y[i], y[i + 1], mu[i], mu[i + 1], - tmpp); - ps.get_p().push_back(tmpp); - } - Poly tmpp; - pieceHermiteInterpol(x[x.size() - 2], x[x.size() - 2 + 1], y[x.size() - 2], - y[x.size() - 2 + 1], mu[x.size() - 2], mu[0], tmpp); - ps.get_p().push_back(tmpp); - return true; -} - -bool Interpolation::spline(const std::vector& x, - const std::vector& y, - PiecewisePoly& ps) { - std::vector lambda(x.size()); - std::vector mu(x.size()); - std::vector m(x.size()); - std::vector dx; - - std::vector::const_iterator i = x.begin(); - std::vector::const_iterator j = i; - for (++j; j != x.end(); ++i, ++j) { - dx.push_back(*j - *i); - } - - lambda.front() = 1; - lambda.back() = 0; - mu.front() = 3 * ((*(++(y.begin()))) - y.front()) / dx.front(); - mu.back() = 3 * (y.back() - (*(++(y.rbegin())))) / dx.back(); - std::vector::iterator pdx0 = dx.begin(); - std::vector::iterator pdx1 = pdx0; - ++pdx1; - std::vector::const_iterator py0 = y.begin(); - std::vector::const_iterator py1 = py0; - ++py1; - std::vector::const_iterator py2 = py1; - ++py2; - std::vector::iterator plambda = lambda.begin(); - ++plambda; - std::vector::iterator pmu = mu.begin(); - ++pmu; - for (; py2 != y.end(); - ++pdx0, ++pdx1, ++py0, ++py1, ++py2, ++plambda, ++pmu) { - *plambda = *pdx0 / (*pdx0 + *pdx1); - *pmu = 3 * ((1 - *plambda) / *pdx0 * (*py1 - *py0) + - *plambda / *pdx1 * (*py2 - *py1)); - } - - // for (unsigned i = 1; i < x.size()-1; ++i){ - // lambda[i] = dx[i-1] / (dx[i-1] + dx[i]); - // mu[i] = 3 * ((1-lambda[i]) / dx[i-1] * (y[i] - y[i-1]) + - // lambda[i] / dx[i] * (y[i+1] - y[i])); - // } - - double bet; - std::vector gam(x.size()); - m[0] = mu[0] / (bet = 2); - for (unsigned j = 1; j < x.size(); ++j) { - gam[j] = lambda[j - 1] / bet; - bet = 2 - (1 - lambda[j]) * gam[j]; - if (bet == 0) { - std::cerr << "a error in triangle solver\n"; - return false; - } - m[j] = (mu[j] - (1 - lambda[j]) * m[j - 1]) / bet; - } - for (int j = x.size() - 2; j >= 0; --j) { - m[j] -= gam[j + 1] * m[j + 1]; - } - - ps.clear(); - ps.get_x() = x; - std::vector::const_iterator px0 = x.begin(); - std::vector::const_iterator px1 = px0; - ++px1; - py0 = y.begin(); - py1 = py0; - ++py1; - std::vector::iterator pm0 = m.begin(); - std::vector::iterator pm1 = pm0; - ++pm1; - for (; px1 != x.end(); ++px0, ++px1, ++py0, ++py1, ++pm0, ++pm1) { - Poly tmpp; - pieceHermiteInterpol(*px0, *px1, *py0, *py1, *pm0, *pm1, tmpp); - ps.get_p().push_back(tmpp); - } - - return true; -} - -bool Interpolation::spline(const std::vector::const_iterator xbegin, - const std::vector::const_iterator xend, - const std::vector::const_iterator ybegin, - PiecewisePoly& ps) { - int xsize = 0; - std::vector::const_iterator itmp = xbegin; - while (itmp++ != xend) { - ++xsize; - } - - std::vector lambda(xsize); - std::vector mu(xsize); - std::vector m(xsize); - std::vector dx; - - // setup linear system - std::vector::const_iterator i = xbegin; - std::vector::const_iterator j = i; - for (++j; j != xend; ++i, ++j) { - dx.push_back(*j - *i); - } - lambda.front() = 1; - lambda.back() = 0; - mu.front() = 3 * ((*(++(itmp = ybegin))) - *ybegin) / dx.front(); - std::vector::iterator pdx0 = dx.begin(); - std::vector::iterator pdx1 = pdx0; - ++pdx1; - std::vector::const_iterator py0 = ybegin; - std::vector::const_iterator py1 = py0; - ++py1; - std::vector::const_iterator py2 = py1; - ++py2; - std::vector::iterator plambda = lambda.begin(); - ++plambda; - std::vector::iterator pmu = mu.begin(); - ++pmu; - for (; pdx1 != dx.end(); - ++pdx0, ++pdx1, ++py0, ++py1, ++py2, ++plambda, ++pmu) { - *plambda = *pdx0 / (*pdx0 + *pdx1); - *pmu = 3 * ((1 - *plambda) / *pdx0 * (*py1 - *py0) + - *plambda / *pdx1 * (*py2 - *py1)); - } - mu.back() = 3 * (*py1 - *py0) / dx.back(); - - // solve tridiangonal linear system - double bet; - std::vector gam(xsize); - m[0] = mu[0] / (bet = 2); - for (int j = 1; j < xsize; ++j) { - gam[j] = lambda[j - 1] / bet; - bet = 2 - (1 - lambda[j]) * gam[j]; - if (bet == 0) { - std::cerr << "a error in triangle solver\n"; - return false; - } - m[j] = (mu[j] - (1 - lambda[j]) * m[j - 1]) / bet; - } - for (int j = xsize - 2; j >= 0; --j) { - m[j] -= gam[j + 1] * m[j + 1]; - } - - // make piecewise polynominal - ps.get_p().clear(); - ps.get_x().resize(xsize); - std::copy(xbegin, xend, ps.get_x().begin()); - std::vector::const_iterator px0 = xbegin; - std::vector::const_iterator px1 = px0; - ++px1; - py0 = ybegin; - py1 = py0; - ++py1; - std::vector::iterator pm0 = m.begin(); - std::vector::iterator pm1 = pm0; - ++pm1; - for (; px1 != xend; ++px0, ++px1, ++py0, ++py1, ++pm0, ++pm1) { - Poly tmpp; - pieceHermiteInterpol(*px0, *px1, *py0, *py1, *pm0, *pm1, tmpp); - ps.get_p().push_back(tmpp); - } - - return true; -} - -// void tridag(float a[], float b[], float c[], float r[], float u[], -// unsigned long n) -// //Solves for a vector u[1..n] the tridiagonal linear set given by equation -// (2.4.1). a[1..n], -// // b[1..n], c[1..n], and r[1..n] are input vectors and are not modified. -// { -// unsigned long j; -// float bet,*gam; -// gam=vector(1,n); //One vector of workspace, gam is needed. -// //If this happens then you should rewrite your equations as a set of order -// N-1, w ith u2 -// //trivially eliminated. -// u[0]=r[0]/(bet=2); -// for (j=1;j<=n;j++) { //Decomposition and forward substitution. -// gam[j]=c[j-1]/bet; -// bet=2-a[j]*gam[j]; -// if (bet == 0.0) //nrerror("Error 2 in tridag"); //Algorithm fails; see be -// u[j]=(r[j]-a[j]*u[j-1])/bet; //low. -// } -// for (j=(n-1);j>=1;j--) -// u[j] -= gam[j+1]*u[j+1]; //Backsubstitution. -// free_vector(gam,1,n); -// } diff --git a/source/md/src/LJInter.cc b/source/md/src/LJInter.cc deleted file mode 100644 index e4c4f1097e..0000000000 --- a/source/md/src/LJInter.cc +++ /dev/null @@ -1,79 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "LJInter.h" - -#include - -#include "common.h" - -LJInter::LJInter(const VALUETYPE& c6_, - const VALUETYPE& c12_, - const VALUETYPE& rc_) - : c6(6. * c6_), c12(12. * c12_), rc(rc_), rc2(rc * rc) { - one_over_6 = 1. / 6.; - one_over_12 = 1. / 12.; - VALUETYPE rc6 = rc2 * rc2 * rc2; - one_over_rc6 = 1. / rc6; - one_over_rc12 = 1. / rc6 / rc6; -} - -void LJInter::lj_inner(VALUETYPE& ae, VALUETYPE& af, const VALUETYPE& r2) { - VALUETYPE rinv = 1. / sqrt(r2); - VALUETYPE rinv2 = rinv * rinv; - VALUETYPE rinv6 = rinv2 * rinv2 * rinv2; - VALUETYPE vvdw6 = c6 * rinv6; - VALUETYPE vvdw12 = c12 * rinv6 * rinv6; - ae = (vvdw12 - c12 * one_over_rc12) * one_over_12 - - (vvdw6 - c6 * one_over_rc6) * one_over_6; - af = (vvdw12 - vvdw6) * rinv2; -} - -void LJInter::compute(VALUETYPE& ener, - vector& force, - vector& virial, - const vector& coord, - const vector& atype, - const SimulationRegion& region, - const vector >& nlist) { - for (unsigned ii = 0; ii < nlist.size(); ++ii) { - for (unsigned _ = 0; _ < nlist[ii].size(); ++_) { - int jj = nlist[ii][_]; - if (jj < ii) { - continue; - } - VALUETYPE diff[3]; - region.diffNearestNeighbor(&coord[ii * 3], &coord[jj * 3], diff); - VALUETYPE r2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]; - if (r2 < rc2) { - VALUETYPE ae, af; - lj_inner(ae, af, r2); - for (int dd = 0; dd < 3; ++dd) { - force[ii * 3 + dd] += af * diff[dd]; - } - for (int dd = 0; dd < 3; ++dd) { - force[jj * 3 + dd] -= af * diff[dd]; - } - ener += ae; - for (int dd0 = 0; dd0 < 3; ++dd0) { - for (int dd1 = 0; dd1 < 3; ++dd1) { - virial[dd0 * 3 + dd1] -= 0.5 * diff[dd0] * af * diff[dd1]; - } - } - } - } - } - - // for (int ii = 0; ii < natoms; ++ii){ - // for (int jj = ii+1; jj < natoms; ++jj){ - // VALUETYPE diff[3]; - // for (int dd = 0; dd < 3; ++dd) diff[dd] = coord[ii*3+dd] - - // coord[jj*3+dd]; diff_pbc (diff, box); VALUETYPE r2 = diff[0] * diff[0] - // + diff[1] * diff[1] + diff[2] * diff[2]; if (r2 < rc2) { - // VALUETYPE ae, af; - // lj_inner (ae, af, r2); - // for (int dd = 0; dd < 3; ++dd) force[ii*3+dd] += af * diff[dd]; - // for (int dd = 0; dd < 3; ++dd) force[jj*3+dd] -= af * diff[dd]; - // ener += ae; - // } - // } - // } -} diff --git a/source/md/src/LJTab.cc b/source/md/src/LJTab.cc deleted file mode 100644 index a605c172fe..0000000000 --- a/source/md/src/LJTab.cc +++ /dev/null @@ -1,28 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "LJTab.h" - -LJTab::LJTab(const VALUETYPE& c6, const VALUETYPE& c12, const VALUETYPE& rc) { - VALUETYPE rcp = rc + 1; - VALUETYPE hh = 2e-3; - int nn = rcp / hh; - vector tab; - VALUETYPE rc6 = rc * rc * rc * rc * rc * rc; - VALUETYPE one_over_rc6 = 1. / rc6; - VALUETYPE one_over_rc12 = 1. / rc6 / rc6; - for (int ii = 0; ii < nn; ++ii) { - VALUETYPE xx = ii * hh; - VALUETYPE value, deriv; - if (xx <= rc) { - VALUETYPE xx3 = xx * xx * xx; - VALUETYPE xx6 = xx3 * xx3; - VALUETYPE xx12 = xx6 * xx6; - value = -c6 / xx6 + c12 / xx12 + c6 * one_over_rc6 - c12 * one_over_rc12; - deriv = -(6. * c6 / xx6 - 12. * c12 / xx12) / xx; - } else { - value = deriv = 0; - } - tab.push_back(value); - tab.push_back(deriv); - } - lj_tab.reinit(rcp, hh, tab); -} diff --git a/source/md/src/MaxShift.cc b/source/md/src/MaxShift.cc deleted file mode 100644 index 484078d51a..0000000000 --- a/source/md/src/MaxShift.cc +++ /dev/null @@ -1,43 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "MaxShift.h" - -#include - -#include "common.h" - -MaxShift::MaxShift(const vector& dcoord, const VALUETYPE& shell_) { - record = dcoord; - shell = shell_; - max_allow2 = shell * 0.5 * shell * 0.5; -} - -VALUETYPE -MaxShift::max_shift2(const vector& coord, - const SimulationRegion& region) { - assert(coord.size() == record.size()); - int natoms = coord.size() / 3; - - VALUETYPE maxv = 0; - - for (int ii = 0; ii < natoms; ++ii) { - VALUETYPE diff[3]; - region.diffNearestNeighbor(&coord[ii * 3], &record[ii * 3], diff); - VALUETYPE r2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]; - if (r2 > maxv) { - maxv = r2; - } - } - - return maxv; -} - -bool MaxShift::rebuild(const vector& coord, - const SimulationRegion& region) { - VALUETYPE maxv2 = max_shift2(coord, region); - if (maxv2 > max_allow2) { - record = coord; - return true; - } else { - return false; - } -} diff --git a/source/md/src/Poly.cpp b/source/md/src/Poly.cpp deleted file mode 100644 index 80db3a139f..0000000000 --- a/source/md/src/Poly.cpp +++ /dev/null @@ -1,284 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "Poly.h" - -bool PiecewisePoly::valid() const { - if (x.size() != p.size() + 1) { - return false; - } - std::vector::const_iterator i = x.begin(); - std::vector::const_iterator j = x.begin(); - for (++j; j != x.end(); ++i, ++j) { - if (*i > *j) { - return false; - } - } - return true; -} - -double PiecewisePoly::value(const double& xx) const { - unsigned begin = 0; - unsigned end = x.size() - 1; - unsigned mid = end / 2; - if (end == begin) { - return 0; - } - while (end - begin > 1) { - if (xx < x[mid]) { - end = mid; - mid = (begin + end) / 2; - } else { - begin = mid; - mid = (begin + end) / 2; - } - } - return p[begin].value(xx); -} - -double PiecewisePoly::value_periodic(const double& xx_) const { - double xx(xx_); - double T = x.back() - x.front(); - if (xx < x.front()) { - while ((xx += T) < x.front()); - } else if (xx >= x.back()) { - while ((xx -= T) >= x.back()); - } - unsigned begin = 0; - unsigned end = x.size() - 1; - unsigned mid = end / 2; - if (end == begin) { - return 0; - } - while (end - begin > 1) { - if (xx < x[mid]) { - end = mid; - mid = (begin + end) / 2; - } else { - begin = mid; - mid = (begin + end) / 2; - } - } - return p[begin].value(xx); -} - -double PiecewisePoly::value(const double& xx, - unsigned& begin, - unsigned& end) const { - if (end <= begin) { - return 0; - } - if (end - begin == 1) { - return p[begin].value(xx); - } - unsigned mid = (begin + end) / 2; - while (end - begin > 1) { - if (xx < x[mid]) { - end = mid; - mid = (begin + end) / 2; - } else { - begin = mid; - mid = (begin + end) / 2; - } - } - return p[begin].value(xx); -} - -void PiecewisePoly::value(const unsigned& xbegin, - const unsigned& xend, - const std::vector& r, - const unsigned& rbegin, - const unsigned& rend, - std::vector& y) const { - unsigned xbegin1 = xbegin; - unsigned xend1 = xend; - if (rend - rbegin <= 1) { - y[rbegin] = value(r[rbegin], xbegin1, xend1); - xbegin1 = xbegin; - xend1 = xend; - y[rend] = value(r[rend], xbegin1, xend1); - } else { - unsigned rmid = (rbegin + rend) / 2; - y[rmid] = value(r[rmid], xbegin1, xend1); - value(xbegin, xend1, r, rbegin, rmid - 1, y); - value(xbegin1, xend, r, rmid + 1, rend, y); - } -} - -// suppose that -void PiecewisePoly::value(const std::vector& r, - std::vector& y) const { - y.resize(r.size()); - value(0, x.size() - 1, r, 0, r.size() - 1, y); -} - -// suppose that -void PiecewisePoly::value_periodic(const std::vector& r, - std::vector& y) const { - std::vector tmpr; - std::vector tmpy; - std::vector > values; - unsigned presentEnd(0), presentStart(0); - double T = x.back() - x.front(); - - while (presentEnd < r.size()) { - tmpr.clear(); - presentStart = presentEnd; - double shift = 0; - if (r[presentStart] < x.front()) { - while (r[presentStart] + (shift += T) < x.front()); - } else if (r[presentStart] >= x.back()) { - while (r[presentStart] + (shift -= T) >= x.back()); - } - while (presentEnd < r.size() && r[presentEnd] + shift >= x.front() && - r[presentEnd] + shift < x.back()) { - tmpr.push_back(r[presentEnd++] + shift); - } - // while (presentEnd < r.size() && r[presentEnd] - r[presentStart] < T){ - // tmpr.push_back (r[presentEnd++]); - // } - // for (unsigned i = 0; i < tmpr.size(); ++i){ - // tmpr[i] += shift; - // } - value(tmpr, tmpy); - values.push_back(tmpy); - } - - y.clear(); - for (unsigned i = 0; i < values.size(); ++i) { - y.insert(y.end(), values[i].begin(), values[i].end()); - } -} - -Poly& Poly::valueLinearPoly(const double& a_, const double& b_, Poly& p) { - std::vector tmp(2, a_); - tmp[0] = b_; - Poly axb(tmp); - p.one(); - p *= a.back(); - for (int i = order - 1; i >= 0; i--) { - (p *= axb) += a[i]; - } - return p; -} - -double Poly::value(const double& x) const { - double value = a[a.size() - 1]; - for (int i = a.size() - 2; i >= 0; --i) { - value = value * x + a[i]; - } - return value; -} - -Poly::Poly() : a(1, 0.), order(0.) {} - -Poly::Poly(const std::vector& out) : a(out) { order = out.size() - 1; } - -void Poly::reinit(const std::vector& out) { - a = out; - order = out.size() - 1; -} - -Poly& Poly::operator=(const Poly& p) { - a = p.a; - order = p.order; - return *this; -} - -Poly& Poly::operator+=(const Poly& p) { - if (p.a.size() > a.size()) { - a.resize(p.a.size(), 0); - order = p.order; - for (unsigned i = 0; i <= order; i++) { - a[i] += p.a[i]; - } - } else { - for (unsigned i = 0; i <= p.order; i++) { - a[i] += p.a[i]; - } - } - return *this; -} - -Poly& Poly::operator+=(const double& b) { - a[0] += b; - return *this; -} - -Poly& Poly::derivative() { - if (order == 0) { - a[0] = 0; - return *this; - } - for (unsigned i = 0; i < order; i++) { - a[i] = a[i + 1] * (i + 1); - } - order--; - a.pop_back(); - return *this; -} - -Poly& Poly::operator*=(const double& scale) { - if (scale == 0) { - order = 0; - a.resize(1); - a[0] = 0; - } else { - for (std::vector::iterator i = a.begin(); i != a.end(); i++) { - *i *= scale; - } - } - return *this; -} - -Poly& Poly::operator*=(const Poly& p) { - std::vector a1(a); - unsigned order1(order); - - order += p.order; - a.resize(order + 1, 0); - - for (std::vector::iterator i = a.begin(); i != a.end(); i++) { - *i *= p.a[0]; - } - if (p.order >= 1) { - for (unsigned i = 1; i <= p.order; i++) { - for (unsigned j = 0; j <= order1; j++) { - a[i + j] += a1[j] * p.a[i]; - } - } - } - return *this; -} - -void Poly::print() { - for (unsigned i = 0; i <= order; i++) { - std::cout << a[i] << '\t'; - } - std::cout << std::endl; -} - -void Poly::print(const std::string& x) { - std::cout << a[0]; - for (unsigned i = 1; i <= order; i++) { - std::cout << " + " << a[i] << x << "^" << i; - } - std::cout << std::endl; -} - -void Poly::printCode(const std::string& x) { - std::cout.precision(16); - if (order == 0) { - std::cout << a[0] << std::endl; - return; - } - - for (unsigned i = 0; i < order - 1; i++) { - std::cout << "("; - } - std::vector::reverse_iterator p = a.rbegin(); - std::cout << *(p++) << " * " << x << " + "; - std::cout << *(p++); - for (; p != a.rend(); p++) { - std::cout << ") * " << x << " + " << *p; - } - std::cout << std::endl; -} diff --git a/source/md/src/RandomGenerator_MT19937.cc b/source/md/src/RandomGenerator_MT19937.cc deleted file mode 100644 index 12f445874a..0000000000 --- a/source/md/src/RandomGenerator_MT19937.cc +++ /dev/null @@ -1,180 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include - -#include "RandomGenerator.h" - -/* - A C-program for MT19937, with initialization improved 2002/1/26. - Coded by Takuji Nishimura and Makoto Matsumoto. - - Before using, initialize the state by using init_genrand(seed) - or init_by_array(init_key, key_length). - - Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, - All rights reserved. - - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions - are met: - - 1. Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. - - 2. Redistributions in binary form must reproduce the above copyright - notice, this list of conditions and the following disclaimer in the - documentation and/or other materials provided with the distribution. - - 3. The names of its contributors may not be used to endorse or promote - products derived from this software without specific prior written - permission. - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS - "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT - LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR - A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER - OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF - LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING - NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS - SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - - - Any feedback is very welcome. - http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html - email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) -*/ - -/* Period parameters */ -#define N 624 -#define M 397 -#define MATRIX_A 0x9908b0dfUL /* constant vector a */ -#define UPPER_MASK 0x80000000UL /* most significant w-r bits */ -#define LOWER_MASK 0x7fffffffUL /* least significant r bits */ - -static unsigned long mt[N]; /* the array for the state vector */ -static int mti = N + 1; /* mti==N+1 means mt[N] is not initialized */ - -// using namespace RandomGenerator_MT19937; - -/* initializes mt[N] with a seed */ -void RandomGenerator_MT19937::init_genrand(unsigned long s) { - mt[0] = s & 0xffffffffUL; - for (mti = 1; mti < N; mti++) { - mt[mti] = (1812433253UL * (mt[mti - 1] ^ (mt[mti - 1] >> 30)) + mti); - /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ - /* In the previous versions, MSBs of the seed affect */ - /* only MSBs of the array mt[]. */ - /* 2002/01/09 modified by Makoto Matsumoto */ - mt[mti] &= 0xffffffffUL; - /* for >32 bit machines */ - } -} - -/* initialize by an array with array-length */ -/* init_key is the array for initializing keys */ -/* key_length is its length */ -/* slight change for C++, 2004/2/26 */ -void RandomGenerator_MT19937::init_by_array(unsigned long init_key[], - int key_length) { - int i, j, k; - init_genrand(19650218UL); - i = 1; - j = 0; - k = (N > key_length ? N : key_length); - for (; k; k--) { - mt[i] = (mt[i] ^ ((mt[i - 1] ^ (mt[i - 1] >> 30)) * 1664525UL)) + - init_key[j] + j; /* non linear */ - mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ - i++; - j++; - if (i >= N) { - mt[0] = mt[N - 1]; - i = 1; - } - if (j >= key_length) { - j = 0; - } - } - for (k = N - 1; k; k--) { - mt[i] = (mt[i] ^ ((mt[i - 1] ^ (mt[i - 1] >> 30)) * 1566083941UL)) - - i; /* non linear */ - mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ - i++; - if (i >= N) { - mt[0] = mt[N - 1]; - i = 1; - } - } - - mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ -} - -/* generates a random number on [0,0xffffffff]-interval */ -unsigned long RandomGenerator_MT19937::genrand_int32(void) { - unsigned long y; - static unsigned long mag01[2] = {0x0UL, MATRIX_A}; - /* mag01[x] = x * MATRIX_A for x=0,1 */ - - if (mti >= N) { /* generate N words at one time */ - int kk; - - if (mti == N + 1) { /* if init_genrand() has not been called, */ - init_genrand(5489UL); /* a default initial seed is used */ - } - - for (kk = 0; kk < N - M; kk++) { - y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK); - mt[kk] = mt[kk + M] ^ (y >> 1) ^ mag01[y & 0x1UL]; - } - for (; kk < N - 1; kk++) { - y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK); - mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1UL]; - } - y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK); - mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ mag01[y & 0x1UL]; - - mti = 0; - } - - y = mt[mti++]; - - /* Tempering */ - y ^= (y >> 11); - y ^= (y << 7) & 0x9d2c5680UL; - y ^= (y << 15) & 0xefc60000UL; - y ^= (y >> 18); - - return y; -} - -/* generates a random number on [0,0x7fffffff]-interval */ -long RandomGenerator_MT19937::genrand_int31(void) { - return (long)(genrand_int32() >> 1); -} - -/* generates a random number on [0,1]-real-interval */ -double RandomGenerator_MT19937::genrand_real1(void) { - return genrand_int32() * (1.0 / 4294967295.0); - /* divided by 2^32-1 */ -} - -/* generates a random number on [0,1)-real-interval */ -double RandomGenerator_MT19937::genrand_real2(void) { - return genrand_int32() * (1.0 / 4294967296.0); - /* divided by 2^32 */ -} - -/* generates a random number on (0,1)-real-interval */ -double RandomGenerator_MT19937::genrand_real3(void) { - return (((double)genrand_int32()) + 0.5) * (1.0 / 4294967296.0); - /* divided by 2^32 */ -} - -/* generates a random number on [0,1) with 53-bit resolution*/ -double RandomGenerator_MT19937::genrand_res53(void) { - unsigned long a = genrand_int32() >> 5, b = genrand_int32() >> 6; - return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0); -} -/* These real versions are due to Isaku Wada, 2002/01/09 added */ diff --git a/source/md/src/Statistics.cc b/source/md/src/Statistics.cc deleted file mode 100644 index 6241d43d06..0000000000 --- a/source/md/src/Statistics.cc +++ /dev/null @@ -1,104 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include -#include -// #include -#include - -#include "Statistics.h" -#include "UnitManager.h" - -template -Statistics::Statistics(const VALUETYPE e_corr_, - const VALUETYPE p_corr_) - : e_corr(e_corr_), p_corr(p_corr_) {} - -template -void Statistics::record(const VALUETYPE& ener, - const vector& virial, - const vector& veloc, - const vector& mass, - const SimulationRegion& region_) { - r_pot_ener = ener; - r_vir.resize(9); - for (unsigned ii = 0; ii < 9; ++ii) { - r_vir[ii] = virial[ii]; - } - // r_box.resize(6); - // for (unsigned ii = 0; ii < 6; ++ii){ - // r_box[ii] = box[ii]; - // } - region.reinitBox(region_.getBoxTensor()); - natoms = mass.size(); - r_kin_ener = 0; - double pref = 0.5 * UnitManager::IntegratorMassConstant; - for (int ii = 0; ii < natoms; ++ii) { - r_kin_ener += pref * mass[ii] * veloc[3 * ii + 0] * veloc[3 * ii + 0]; - r_kin_ener += pref * mass[ii] * veloc[3 * ii + 1] * veloc[3 * ii + 1]; - r_kin_ener += pref * mass[ii] * veloc[3 * ii + 2] * veloc[3 * ii + 2]; - } -} - -template -double Statistics::get_T() const { - return get_ekin() / (natoms * 3. * UnitManager::BoltzmannConstant) * 2.; -} - -template -double Statistics::get_V() const { - // return (r_box[1] - r_box[0]) * (r_box[3] - r_box[2]) * (r_box[5] - - // r_box[4]); - return region.getVolume(); -} - -template -double Statistics::get_P() const { - return (get_ekin() - (r_vir[0] + r_vir[4] + r_vir[8])) * 2. / 3. / get_V() * - UnitManager::PressureConstant + - p_corr; -} - -template -void Statistics::print(ostream& os, - const int& step, - const double time) const { - char tmps[65536]; - sprintf(tmps, - "%13.4f %15.6f %15.6f %15.6f %15.6f %15.6f %15.6f %15.6f %15.6f\n", - time, get_ekin(), get_epot(), get_ekin() + get_epot(), get_T(), - get_P(), r_vir[0], r_vir[4], r_vir[8]); - os << tmps; - // os << setw(7) << setprecision(6) << time << setprecision (8) << setfill (' - // ') - // << setw(15) << get_ekin() << " " - // << setw(15) << get_epot() << " " - // << setw(15) << get_ekin() + get_epot() << " " - // << setw(15) << get_T() << " " - // << setw(15) << get_P() << " " - // << setw(15) << r_vir[0] << " " - // << setw(15) << r_vir[4] << " " - // << setw(15) << r_vir[8] << " " - // << endl; -} - -template -void Statistics::print_head(ostream& os) const { - char tmps[65536]; - sprintf(tmps, "#%12s %15s %15s %15s %15s %15s %15s %15s %15s\n", "time", - "Kinetic", "Potential", "E_tot", "Temperature", "Pressure", "Vxx", - "Vyy", "Vzz"); - os << tmps; - // os << "#"; - // os << setw(6) << "time" << setfill (' ') - // << setw(15) << "Kinetic" << " " - // << setw(15) << "Potential" << " " - // << setw(15) << "E_tot" << " " - // << setw(15) << "Temperature" << " " - // << setw(15) << "Pressure" << " " - // << setw(15) << "Vxx" << " " - // << setw(15) << "Vyy" << " " - // << setw(15) << "Vzz" << " " - // << endl; -} - -template class Statistics; -template class Statistics; diff --git a/source/md/src/StringSplit.cpp b/source/md/src/StringSplit.cpp deleted file mode 100644 index a34c6a43ab..0000000000 --- a/source/md/src/StringSplit.cpp +++ /dev/null @@ -1,40 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "StringSplit.h" - -void StringOperation::split(const std::string& in, - std::vector& out) { - std::istringstream iss(in); - out.clear(); - - do { - std::string sub; - iss >> sub; - out.push_back(sub); - // std::vector tokens; - // tokens.push_back (" "); - // tokens.push_back ("\t"); - // std::copy(std::istream_iterator(iss), - // std::istream_iterator(), - // std::back_inserter >(tokens)); - } while (iss); - - out.pop_back(); -} - -void StringOperation::split(const std::string& in, - const std::string& delimiter, - std::vector& out) { - size_t pos = 0; - size_t len = delimiter.length(); - std::string s(in); - std::string token; - - while ((pos = s.find(delimiter)) != std::string::npos) { - token = s.substr(0, pos); - out.push_back(token); - s.erase(0, pos + len); - } - if (!s.empty()) { - out.push_back(s); - } -} diff --git a/source/md/src/TF.cc b/source/md/src/TF.cc deleted file mode 100644 index dc37f34c26..0000000000 --- a/source/md/src/TF.cc +++ /dev/null @@ -1,58 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "TF.h" - -#include - -#include "Interpolation.h" -#include "TableFileLoader.h" - -TF::TF(const string& filename) { - vector > tmpdata; - TableFileLoader tfl(filename.c_str()); - tfl.setColumns({1, 3}); - tfl.loadAll(tmpdata); - data = tmpdata[1]; - hh = tmpdata[0][1] - tmpdata[0][0]; - xup = tmpdata[0].back(); - xup *= b2m_l; - hh *= b2m_l; - for (unsigned ii = 0; ii < data.size(); ++ii) { - data[ii] *= b2m_e / b2m_l; - } -} - -VALUETYPE -TF::meas(const VALUETYPE& xx) const { - VALUETYPE ff = 0; - if (xx >= xup) { - ff = 0; - } else { - int posi = int(xx / hh); - if (posi < 0) { - posi = 0; - } else if (posi >= data.size() - 1) { - posi = data.size() - 2; - } - Poly p; - Interpolation::pieceLinearInterpol(posi * hh, (posi + 1) * hh, data[posi], - data[posi + 1], p); - ff = p.value(xx); - } - return ff; -} - -void TF::apply(vector& dforce, - const vector& dcoord, - const AdWeight& adw) const { - vector weight, weight_x; - adw.atom_weight(weight, weight_x, dcoord); - vector center = adw.get_center(); - - for (unsigned ii = 0; ii < weight_x.size(); ++ii) { - VALUETYPE ff = meas(weight_x[ii]); - if (dcoord[ii * 3] < center[0]) { - ff = -ff; - } - dforce[ii * 3] += ff; - } -} diff --git a/source/md/src/TableFileLoader.cpp b/source/md/src/TableFileLoader.cpp deleted file mode 100644 index 99fec6f9cb..0000000000 --- a/source/md/src/TableFileLoader.cpp +++ /dev/null @@ -1,98 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "TableFileLoader.h" - -#include -#include - -#include "StringSplit.h" - -#define MaxLineLength 65536 - -using namespace std; - -TableFileLoader::TableFileLoader(const char* file) : every(1) { reinit(file); } - -unsigned TableFileLoader::getNumbColumns() { - char valueline[MaxLineLength]; - while (data.getline(valueline, MaxLineLength)) { - if (valueline[0] == '#' || valueline[0] == '@') { - continue; - } - break; - } - if (data.eof()) { - return 0; - } else if (!data.good()) { - cerr << "error file reading state!" << endl; - throw; - } - vector words; - StringOperation::split(string(valueline), words); - - data.close(); - reinit(file.c_str()); - return words.size(); -} - -void TableFileLoader::reinit(const char* file_) { - file = string(file_); - data.open(file.c_str()); - if (!data) { - cerr << "cannot open file \"" << file << "\"" << endl; - throw; - } - count_read = 0; - // inter_cols.push_back (0); -} - -void TableFileLoader::setColumns(const vector& cols) { - inter_cols = cols; - for (unsigned ii = 0; ii < inter_cols.size(); ++ii) { - if (inter_cols[ii] == 0) { - cerr << "invalid col index, should be larger than 0" << endl; - throw; - } - inter_cols[ii] -= 1; - } -} - -void TableFileLoader::setEvery(const unsigned every_) { every = every_; } - -bool TableFileLoader::loadLine(vector& odata) { - char valueline[MaxLineLength]; - - while (data.getline(valueline, MaxLineLength)) { - if (valueline[0] == '#' || valueline[0] == '@') { - continue; - } else if (count_read++ % every == 0) { - break; - } - } - - if (data.eof()) { - return false; - } else if (!data.good()) { - cerr << "error file reading state!" << endl; - throw; - } - - vector words; - StringOperation::split(string(valueline), words); - odata.resize(inter_cols.size()); - - for (unsigned ii = 0; ii < inter_cols.size(); ++ii) { - odata[ii] = atof(words[inter_cols[ii]].c_str()); - } - - return true; -} - -void TableFileLoader::loadAll(vector >& odata) { - odata.resize(inter_cols.size()); - vector line; - while (loadLine(line)) { - for (unsigned ii = 0; ii < inter_cols.size(); ++ii) { - odata[ii].push_back(line[ii]); - } - } -} diff --git a/source/md/src/Tabulated.cc b/source/md/src/Tabulated.cc deleted file mode 100644 index 6e9777ea29..0000000000 --- a/source/md/src/Tabulated.cc +++ /dev/null @@ -1,175 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "Tabulated.h" - -#include -#include - -#include "UnitManager.h" -#include "common.h" - -Tabulated::Tabulated(const VALUETYPE rc, - const VALUETYPE hh, - const vector &tab) { - reinit(rc, hh, tab); -} - -void Tabulated::reinit(const VALUETYPE rc, - const VALUETYPE hh, - const vector &tab) { - int numbFunc = 1; - int stride = numbFunc * 4; - int mystride = numbFunc * 2; - unsigned tableLength = tab.size() / mystride; - - hi = 1. / hh; - rc2 = rc * rc; - - data.resize(static_cast(tableLength) * stride); - - int ii; - for (ii = 0; ii < tableLength - 1; ++ii) { - const double &v0(tab[ii * mystride + 0]); - const double &f0(tab[ii * mystride + 1]); - const double &v1(tab[(ii + 1) * mystride + 0]); - const double &f1(tab[(ii + 1) * mystride + 1]); - VALUETYPE &dv(data[ii * stride + 0]); - VALUETYPE &df(data[ii * stride + 1]); - VALUETYPE &dg(data[ii * stride + 2]); - VALUETYPE &dh(data[ii * stride + 3]); - dv = v0; - df = -f0 * hh; - dg = 3 * (v1 - v0) + (f1 + 2 * f0) * hh; - dh = -2 * (v1 - v0) - (f1 + f0) * hh; - } - { - const double &v0(tab[ii * mystride + 0]); - const double &f0(tab[ii * mystride + 1]); - VALUETYPE &dv(data[ii * stride + 0]); - VALUETYPE &df(data[ii * stride + 1]); - VALUETYPE &dg(data[ii * stride + 2]); - VALUETYPE &dh(data[ii * stride + 3]); - dv = v0; - df = -f0 * hh; - dg = 0; - dh = 0; - } -} - -inline void Tabulated::compute_posi(int &idx, - VALUETYPE &eps, - const VALUETYPE &rr) { - VALUETYPE rt = rr * hi; - idx = int(rt); - eps = rt - idx; -} - -inline void Tabulated::tb_inner(VALUETYPE &ae, - VALUETYPE &af, - const VALUETYPE &r2) { - if (r2 > rc2) { - ae = af = 0; - return; - } - - VALUETYPE rr = sqrt(r2); - int idx; - VALUETYPE eps; - compute_posi(idx, eps, rr); - idx *= 4; - - VALUETYPE table_param[4]; - for (int ii = 0; ii < 4; ++ii) { - table_param[ii] = data[ii + idx]; - } - const VALUETYPE &Y(table_param[0]); - const VALUETYPE &F(table_param[1]); - const VALUETYPE &G(table_param[2]); - const VALUETYPE &H(table_param[3]); - - VALUETYPE Heps = eps * H; - VALUETYPE Fp = (F + eps * (G + Heps)); - VALUETYPE FF = (Fp + (eps * (G + (Heps + Heps)))); - - af = FF * hi; - af = -af / rr; - ae = (Y + (eps * Fp)); -} - -void Tabulated::compute(VALUETYPE &ener, - vector &force, - vector &virial, - const vector &coord, - const vector &atype, - const SimulationRegion ®ion, - const vector > &nlist) { - for (unsigned ii = 0; ii < nlist.size(); ++ii) { - for (unsigned _ = 0; _ < nlist[ii].size(); ++_) { - int jj = nlist[ii][_]; - if (jj < ii) { - continue; - } - VALUETYPE diff[3]; - region.diffNearestNeighbor(&coord[ii * 3], &coord[jj * 3], diff); - VALUETYPE r2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]; - if (r2 < rc2) { - VALUETYPE ae, af; - tb_inner(ae, af, r2); - for (int dd = 0; dd < 3; ++dd) { - force[ii * 3 + dd] += af * diff[dd]; - } - for (int dd = 0; dd < 3; ++dd) { - force[jj * 3 + dd] -= af * diff[dd]; - } - ener += ae; - for (int dd0 = 0; dd0 < 3; ++dd0) { - for (int dd1 = 0; dd1 < 3; ++dd1) { - virial[dd0 * 3 + dd1] -= 0.5 * diff[dd0] * af * diff[dd1]; - } - } - } - } - } -} - -void Tabulated::compute(VALUETYPE &ener, - vector &force, - vector &virial, - const vector &coord, - const vector &charge, - const vector &atype, - const SimulationRegion ®ion, - const vector > &nlist) { - for (unsigned ii = 0; ii < nlist.size(); ++ii) { - for (unsigned _ = 0; _ < nlist[ii].size(); ++_) { - int jj = nlist[ii][_]; - if (jj < ii) { - continue; - } - VALUETYPE diff[3]; - region.diffNearestNeighbor(&coord[ii * 3], &coord[jj * 3], diff); - VALUETYPE r2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]; - if (r2 < rc2) { - VALUETYPE ae, af; - tb_inner(ae, af, r2); - { - VALUETYPE qiqj = - charge[ii] * charge[jj] * UnitManager::ElectrostaticConvertion; - ae *= qiqj; - af *= qiqj; - } - for (int dd = 0; dd < 3; ++dd) { - force[ii * 3 + dd] += af * diff[dd]; - } - for (int dd = 0; dd < 3; ++dd) { - force[jj * 3 + dd] -= af * diff[dd]; - } - ener += ae; - for (int dd0 = 0; dd0 < 3; ++dd0) { - for (int dd1 = 0; dd1 < 3; ++dd1) { - virial[dd0 * 3 + dd1] -= 0.5 * diff[dd0] * af * diff[dd1]; - } - } - } - } - } -} diff --git a/source/md/src/Trajectory.cc b/source/md/src/Trajectory.cc deleted file mode 100644 index 8da1cff67e..0000000000 --- a/source/md/src/Trajectory.cc +++ /dev/null @@ -1,141 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "Trajectory.h" - -#include -#include - -#include -#include - -bool XtcSaver::reinit(const char *filename, const int &natoms_) { - char tmpname[2048]; - strncpy(tmpname, filename, 2047); - - xd = xdrfile_open(filename, "w"); - if (xd == NULL) { - std::cerr << "cannot open file " << filename << std::endl; - return false; - } - natoms = natoms_; - - xx = (rvec *)malloc(sizeof(rvec) * natoms); - inited = true; - return true; -} - -XtcSaver::~XtcSaver() { clear(); } - -XtcSaver::XtcSaver(const char *filename, const int &natoms_) - : inited(false), prec(1000) { - reinit(filename, natoms_); -} - -void XtcSaver::clear() { - if (inited) { - free(xx); - xdrfile_close(xd); - inited = false; - } -} - -void XtcSaver::save(const int &step, - const double &time, - const vector > &frame, - const vector &box) { - assert(box.size() == 9); - assert(inited); - matrix tmpBox; - for (int dd0 = 0; dd0 < 3; ++dd0) { - for (int dd1 = 0; dd1 < 3; ++dd1) { - tmpBox[dd0][dd1] = 0; - } - } - for (int dd = 0; dd < 3; ++dd) { - tmpBox[dd][dd] = box[3 * dd + dd]; - } - for (int ii = 0; ii < frame.size(); ++ii) { - for (int dd = 0; dd < 3; ++dd) { - xx[ii][dd] = frame[ii][dd]; - } - } - write_xtc(xd, natoms, step, time, tmpBox, xx, prec); -} - -bool TrrSaver::reinit(const char *filename, const int &natoms_) { - char tmpname[2048]; - strncpy(tmpname, filename, 2047); - - xd = xdrfile_open(filename, "w"); - if (xd == NULL) { - std::cerr << "cannot open file " << filename << std::endl; - return false; - } - natoms = natoms_; - - xx = (rvec *)malloc(sizeof(rvec) * natoms); - vv = (rvec *)malloc(sizeof(rvec) * natoms); - ff = (rvec *)malloc(sizeof(rvec) * natoms); - for (int ii = 0; ii < natoms; ++ii) { - for (int dd = 0; dd < 3; ++dd) { - vv[ii][dd] = 0; - ff[ii][dd] = 0; - } - } - inited = true; - return true; -} - -TrrSaver::~TrrSaver() { clear(); } - -TrrSaver::TrrSaver(const char *filename, const int &natoms_) - : inited(false), lambda(0) { - reinit(filename, natoms_); -} - -void TrrSaver::clear() { - if (inited) { - free(xx); - free(vv); - free(ff); - xdrfile_close(xd); - inited = false; - } -} - -void TrrSaver::save(const int &step, - const double &time, - const vector > &ixx, - const vector > &ivv, - const vector > &iff, - const vector &box) { - assert(box.size() == 9); - assert(inited); - matrix tmpBox; - for (int dd0 = 0; dd0 < 3; ++dd0) { - for (int dd1 = 0; dd1 < 3; ++dd1) { - tmpBox[dd0][dd1] = box[3 * dd0 + dd1]; - } - } - for (int ii = 0; ii < ixx.size(); ++ii) { - for (int dd = 0; dd < 3; ++dd) { - xx[ii][dd] = ixx[ii][dd]; - } - } - for (int ii = 0; ii < natoms; ++ii) { - for (int dd = 0; dd < 3; ++dd) { - vv[ii][dd] = 0; - ff[ii][dd] = 0; - } - } - for (int ii = 0; ii < ivv.size(); ++ii) { - for (int dd = 0; dd < 3; ++dd) { - vv[ii][dd] = ivv[ii][dd]; - } - } - for (int ii = 0; ii < iff.size(); ++ii) { - for (int dd = 0; dd < 3; ++dd) { - ff[ii][dd] = iff[ii][dd]; - } - } - write_trr(xd, natoms, step, time, lambda, tmpBox, xx, vv, ff); -} diff --git a/source/md/src/UnitManager.cc b/source/md/src/UnitManager.cc deleted file mode 100644 index 1982b67de3..0000000000 --- a/source/md/src/UnitManager.cc +++ /dev/null @@ -1,29 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "UnitManager.h" - -#include - -// unit independent constants -double UnitManager::Degree2Radian = M_PI / 180.; -double UnitManager::Radian2Degree = 180. / M_PI; -// unit dependent -double UnitManager::IntegratorMassConstant = 1.; -double UnitManager::PressureConstant = 16.60539040; -double UnitManager::BoltzmannConstant = 8.31445986144858164e-3; -double UnitManager::ElectrostaticConvertion = 138.93545756169981341199; - -string UnitManager::unit_names[] = {"biology", "metal", "unitless"}; - -void UnitManager::set(const string& unit) { - if (unit == "metal") { - IntegratorMassConstant = 1.03642695707516506071e-4; - PressureConstant = 1.602176621e6; - BoltzmannConstant = 8.6173303e-5; - ElectrostaticConvertion = 14.39964535475696995031; - } else if (unit == "unitless") { - IntegratorMassConstant = 1.; - PressureConstant = 1.; - BoltzmannConstant = 1.; - ElectrostaticConvertion = 1.; - } -} diff --git a/source/md/src/XyzFileManager.cc b/source/md/src/XyzFileManager.cc deleted file mode 100644 index c29a456b9f..0000000000 --- a/source/md/src/XyzFileManager.cc +++ /dev/null @@ -1,108 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "XyzFileManager.h" - -#include - -#include "StringSplit.h" -// #include -#include - -#include - -void XyzFileManager::read(const string& file, - vector& atom_name, - vector >& posi, - vector >& velo, - vector >& forc, - vector& boxsize) { - getBoxSize(file, boxsize); - - posi.clear(); - velo.clear(); - - ifstream data0(file.c_str()); - if (!data0.is_open()) { - cerr << "cannot open file " << file << endl; - exit(1); - } - - string valueline; - vector words; - words.reserve(10); - string tmpname; - vector tmpp(3); - vector tmpv(3); - vector tmpf(3); - std::getline(data0, valueline); - long long int numb_atom = atoll(valueline.c_str()); - std::getline(data0, valueline); - - for (long long int ii = 0; ii < numb_atom; ++ii) { - std::getline(data0, valueline); - StringOperation::split(string(valueline), words); - if (words.size() == 10) { - tmpp[0] = atof(words[1 + 0].c_str()); - tmpp[1] = atof(words[1 + 1].c_str()); - tmpp[2] = atof(words[1 + 2].c_str()); - tmpv[0] = atof(words[1 + 3].c_str()); - tmpv[1] = atof(words[1 + 4].c_str()); - tmpv[2] = atof(words[1 + 5].c_str()); - tmpf[0] = atof(words[1 + 6].c_str()); - tmpf[1] = atof(words[1 + 7].c_str()); - tmpf[2] = atof(words[1 + 8].c_str()); - posi.push_back(tmpp); - velo.push_back(tmpv); - forc.push_back(tmpf); - atom_name.push_back(words[0]); - } else if (words.size() == 7) { - tmpp[0] = atof(words[1 + 0].c_str()); - tmpp[1] = atof(words[1 + 1].c_str()); - tmpp[2] = atof(words[1 + 2].c_str()); - tmpv[0] = atof(words[1 + 3].c_str()); - tmpv[1] = atof(words[1 + 4].c_str()); - tmpv[2] = atof(words[1 + 5].c_str()); - posi.push_back(tmpp); - velo.push_back(tmpv); - atom_name.push_back(words[0]); - } else if (words.size() == 4) { - tmpp[0] = atof(words[1 + 0].c_str()); - tmpp[1] = atof(words[1 + 1].c_str()); - tmpp[2] = atof(words[1 + 2].c_str()); - posi.push_back(tmpp); - atom_name.push_back(words[0]); - } else { - cerr << "XyzFileManager::read: wrong format, line has " << words.size() - << " words" << endl; - exit(1); - } - } -} - -void XyzFileManager::getBoxSize(const string& file, vector& boxsize) { - ifstream data0(file.c_str()); - if (!data0.is_open()) { - cerr << "cannot open file " << file << endl; - } - string valueline; - vector words; - words.reserve(9); - std::getline(data0, valueline); - std::getline(data0, valueline); - StringOperation::split(valueline, words); - - boxsize.resize(9); - fill(boxsize.begin(), boxsize.end(), 0.); - if (words.size() == 3) { - for (int ii = 0; ii < 3; ++ii) { - boxsize[3 * ii + ii] = atof(words[ii].c_str()); - } - } else if (words.size() == 9) { - for (int ii = 0; ii < 9; ++ii) { - boxsize[ii] = atof(words[ii].c_str()); - } - } else { - cerr << "XyzFileManager::getBoxSize: wrong format, line has " - << words.size() << " words" << endl; - exit(1); - } -} diff --git a/source/md/src/ZM.cc b/source/md/src/ZM.cc deleted file mode 100644 index a1495ee650..0000000000 --- a/source/md/src/ZM.cc +++ /dev/null @@ -1,83 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "ZM.h" - -#include -#include - -#include "UnitManager.h" -#include "common.h" - -ZM::ZM(const int& order, const VALUETYPE& alpha, const VALUETYPE& rc) - : potzm(order, alpha, rc) { - VALUETYPE rcp = rc + 2; - VALUETYPE hh = 2e-3; - int nn = rcp / hh; - vector tab; - - for (int ii = 0; ii < nn; ++ii) { - VALUETYPE xx = ii * hh; - VALUETYPE value, deriv; - if (xx <= rc) { - value = potzm.pot(xx); - deriv = potzm.mpotp(xx); - } else { - value = deriv = 0; - } - tab.push_back(value); - tab.push_back(deriv); - } - zm_tab.reinit(rcp, hh, tab); -} - -VALUETYPE -ZM::e_corr(const vector& charge) const { - double sum = 0; - sum += potzm.energyCorr(charge); - return sum; -} - -inline void ZM::ex_inner(VALUETYPE& ae, VALUETYPE& af, const VALUETYPE& r2) { - VALUETYPE r1 = sqrt(r2); - ae = 1. / r1; - af = 1. / (r2 * r1); -} - -void ZM::exclude(VALUETYPE& ener, - vector& force, - vector& virial, - const vector& coord, - const vector& charge, - const vector& atype, - const SimulationRegion& region, - const vector& elist) { - for (unsigned _ = 0; _ < elist.size(); _ += 2) { - int ii = elist[_]; - int jj = elist[_ + 1]; - VALUETYPE diff[3]; - region.diffNearestNeighbor(&coord[ii * 3], &coord[jj * 3], diff); - VALUETYPE r2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]; - VALUETYPE ae, af; - ex_inner(ae, af, r2); - // VALUETYPE ae1, af1; - // zm_tab.tb_inner (ae1, af1, r2); - // cout << ae << " " << ae1 << endl; - { - VALUETYPE qiqj = - charge[ii] * charge[jj] * UnitManager::ElectrostaticConvertion; - ae *= qiqj; - af *= qiqj; - } - for (int dd = 0; dd < 3; ++dd) { - force[ii * 3 + dd] -= af * diff[dd]; - } - for (int dd = 0; dd < 3; ++dd) { - force[jj * 3 + dd] += af * diff[dd]; - } - ener -= ae; - for (int dd0 = 0; dd0 < 3; ++dd0) { - for (int dd1 = 0; dd1 < 3; ++dd1) { - virial[dd0 * 3 + dd1] += 0.5 * diff[dd0] * af * diff[dd1]; - } - } - } -} diff --git a/source/md/src/ZMFunctions.cpp b/source/md/src/ZMFunctions.cpp deleted file mode 100644 index fd85e04bdd..0000000000 --- a/source/md/src/ZMFunctions.cpp +++ /dev/null @@ -1,213 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -#include "ZMFunctions.h" - -#include -#include - -#include "UnitManager.h" - -#define M_inv2 (0.5) -#define M_inv4 (0.25) -#define M_inv8 (0.125) -#define M_inv16 (0.06250000000000000000) -#define M_inv48 (.02083333333333333333) - -static double f(const double& r) { return 1. / r; } - -static double D1f(const double& r) { return -1. / (r * r); } - -static double D2f(const double& r) { return 2. / (r * r * r); } - -static double D3f(const double& r) { return -6. / (r * r * r * r); } - -static double D4f(const double& r) { return 24. / (r * r * r * r * r); } - -static double g(const double& alpha, const double& r) { - return erfc(alpha * r); -} - -static double D1g(const double& alpha, const double& r) { - double tmp = alpha * r; - return -M_2_SQRTPI * alpha * exp(-tmp * tmp); -} - -static double D2g(const double& alpha, const double& r) { - double tmp = alpha * r; - return M_2_SQRTPI * 2 * alpha * alpha * alpha * r * exp(-tmp * tmp); -} - -static double D3g(const double& alpha, const double& r) { - double tmp = alpha * r; - return M_2_SQRTPI * 2 * alpha * alpha * alpha * (1. - 2. * tmp * tmp) * - exp(-tmp * tmp); -} - -static double D4g(const double& alpha, const double& r) { - double tmp = alpha * r; - double alpha5 = alpha * alpha; - alpha5 = alpha5 * alpha5 * alpha; - return M_2_SQRTPI * 4. * alpha5 * (-3. + 2. * tmp * tmp) * r * - exp(-tmp * tmp); -} - -double ZeroMultipole::funcV(const double& alpha, const double& r) { - return f(r) * g(alpha, r); -} - -double ZeroMultipole::funcD1V(const double& alpha, const double& r) { - return D1f(r) * g(alpha, r) + f(r) * D1g(alpha, r); -} - -double ZeroMultipole::funcD2V(const double& alpha, const double& r) { - return D2f(r) * g(alpha, r) + 2. * D1f(r) * D1g(alpha, r) + - f(r) * D2g(alpha, r); -} - -double ZeroMultipole::funcD3V(const double& alpha, const double& r) { - return D3f(r) * g(alpha, r) + 3. * D2f(r) * D1g(alpha, r) + - 3. * D1f(r) * D2g(alpha, r) + f(r) * D3g(alpha, r); -} - -double ZeroMultipole::funcD4V(const double& alpha, const double& r) { - return D4f(r) * g(alpha, r) + 4. * D3f(r) * D1g(alpha, r) + - 6. * D2f(r) * D2g(alpha, r) + 4. * D1f(r) * D3g(alpha, r) + - f(r) * D4g(alpha, r); -} - -void ZeroMultipole::calCoefficients(const int& ll, - const double& alpha, - const double& rc, - vector& coeff) { - coeff.clear(); - coeff.resize(ll + 1); - double b0, b1, b2, b3, b4; - double invrc, invrc2, invrc3, invrc4; - double rc2; - - switch (ll) { - case 0: - b0 = funcV(alpha, rc); - coeff[0] = b0; - break; - case 1: - b0 = funcV(alpha, rc); - b1 = funcD1V(alpha, rc); - coeff[0] = b0 - M_inv2 * b1 * rc; - coeff[1] = M_inv2 * b1 / rc; - break; - case 2: - b0 = funcV(alpha, rc); - b1 = funcD1V(alpha, rc); - b2 = funcD2V(alpha, rc); - invrc = 1. / rc; - coeff[0] = M_inv8 * b2 * rc * rc - 5. * M_inv8 * b1 * rc + b0; - coeff[1] = 3. * M_inv4 * b1 * invrc - M_inv4 * b2; - coeff[2] = - M_inv8 * b2 * invrc * invrc - M_inv8 * b1 * invrc * invrc * invrc; - break; - case 3: - b0 = funcV(alpha, rc); - b1 = funcD1V(alpha, rc); - b2 = funcD2V(alpha, rc); - b3 = funcD3V(alpha, rc); - invrc = 1. / rc; - invrc2 = invrc * invrc; - coeff[0] = -M_inv48 * b3 * rc * rc * rc + 3. * M_inv16 * b2 * rc * rc - - 11. * M_inv16 * b1 * rc + b0; - coeff[1] = - 15. * M_inv16 * b1 * invrc - 7. * M_inv16 * b2 + M_inv16 * b3 * rc; - coeff[2] = 5. * M_inv16 * b2 * invrc2 - - 5. * M_inv16 * b1 * invrc2 * invrc - M_inv16 * b3 * invrc; - coeff[3] = M_inv16 * b1 * invrc2 * invrc2 * invrc - - M_inv16 * b2 * invrc2 * invrc2 + M_inv48 * b3 * invrc2 * invrc; - break; - case 4: - b0 = funcV(alpha, rc); - b1 = funcD1V(alpha, rc); - b2 = funcD2V(alpha, rc); - b3 = funcD3V(alpha, rc); - b4 = funcD4V(alpha, rc); - rc2 = rc * rc; - invrc = 1. / rc; - invrc2 = invrc * invrc; - invrc3 = invrc2 * invrc; - invrc4 = invrc2 * invrc2; - coeff[0] = 1. / 384. * b4 * rc2 * rc2 - 7. / 192. * b3 * rc2 * rc + - 29. / 128. * b2 * rc2 - 93. / 128. * b1 * rc + b0; - coeff[1] = 35. / 32. * b1 * invrc - 19. / 32. * b2 - 1. / 96. * b4 * rc2 + - M_inv8 * b3 * rc; - coeff[2] = 1. / 64. * b4 - 35. / 64. * b1 * invrc3 + - 35. / 64. * b2 * invrc2 - 5. / 32. * b3 * invrc; - coeff[3] = 7. / 32. * b1 * invrc4 * invrc - 7. / 32. * b2 * invrc4 + - 1. / 12. * b3 * invrc3 - 1. / 96. * b4 * invrc2; - coeff[4] = 5. / 128. * b2 * invrc4 * invrc2 - - 5. / 128. * b1 * invrc4 * invrc3 - - 1. / 64. * b3 * invrc4 * invrc + 1. / 384 * b4 * invrc4; - break; - default: - cerr << "ll larger than 4 is not implemented" << endl; - break; - } -} - -ZeroMultipole::Potential::Potential() : alpha(0), rc(1.0), ll(0) { - calCoefficients(ll, alpha, rc, coeff); -} - -ZeroMultipole::Potential::Potential(const int& ll, - const double& alpha, - const double& rc) { - reinit(ll, alpha, rc); -} - -void ZeroMultipole::Potential::reinit(const int& ll_, - const double& alpha_, - const double& rc_) { - ll = ll_; - alpha = alpha_; - rc = rc_; - calCoefficients(ll, alpha, rc, coeff); -} - -double ZeroMultipole::Potential::pot(const double& rr) { - if (rr > rc) { - return 0.; - } - double tmp0 = funcV(alpha, rr); - // double tmp0 = 0.; - double tmp1 = coeff.back(); - for (int ii = ll - 1; ii >= 0; --ii) { - tmp1 = tmp1 * rr * rr + coeff[ii]; - } - return tmp0 - tmp1; -} - -double ZeroMultipole::Potential::ulpot(const double& rr) { - return pot(rr) + coeff[0]; -} - -double ZeroMultipole::Potential::mpotp(const double& rr) { - if (rr > rc) { - return 0.; - } - double tmp0 = -funcD1V(alpha, rr); - double tmp1 = 2 * ll * coeff[ll]; - for (int ii = ll - 1; ii >= 1; --ii) { - tmp1 = tmp1 * rr * rr + coeff[ii] * 2 * ii; - } - return tmp0 + tmp1 * rr; -} - -double ZeroMultipole::Potential::mulpotp(const double& rr) { return mpotp(rr); } - -double ZeroMultipole::Potential::energyCorr( - const vector& charges) const { - double sum = 0.; - double factor = UnitManager::ElectrostaticConvertion; - for (unsigned ii = 0; ii < charges.size(); ++ii) { - sum += charges[ii] * charges[ii]; - } - - // return - (coeff[0] * 0.5 + alpha / sqrt(M_PI)) * sum; - return -(coeff[0] * 0.5 + alpha / sqrt(M_PI)) * sum * factor; -}