Remove temperature from interface

CMakeLists.txt

@@ -39,14 +39,14 @@ if (NOT nanobind_FOUND)
   FetchContent_Declare(
     nanobind
     GIT_REPOSITORY https://github.com/wjakob/nanobind/
-    GIT_TAG        v2.6.0
+    GIT_TAG        v2.7.0
   )
   FetchContent_MakeAvailable(nanobind)
 endif()
 FetchContent_Declare(
   uammd
   GIT_REPOSITORY https://github.com/RaulPPelaez/uammd/
-  GIT_TAG        v2.8.0
+  GIT_TAG        v2.8.5
   EXCLUDE_FROM_ALL
 )
 FetchContent_MakeAvailable(uammd)
@@ -61,7 +61,6 @@ include_directories(${lanczos_SOURCE_DIR}/include)
 
 find_package(BLAS REQUIRED)
 find_package(LAPACK REQUIRED)
-
 include_directories(${CMAKE_SOURCE_DIR}/include)
 if(CMAKE_CUDA_COMPILER_ID STREQUAL "NVIDIA")
   set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} --expt-relaxed-constexpr -extended-lambda")

examples/cpp/example.cu

@@ -71,7 +71,6 @@ int main(){
   par.hydrodynamicRadius = {1};
   par.viscosity = 1;
   par.tolerance = 1e-4;
-  par.temperature = 1.0;
   par.includeAngular = false;
 
   //Create two different solvers

examples/cpp/example_gpu.cu

@@ -82,7 +82,6 @@ int main() {
     par.hydrodynamicRadius = {1};
     par.viscosity = 1;
     par.tolerance = 1e-4;
-    par.temperature = 1.0;
     par.includeAngular = false;
 
     // Create two different solvers

examples/python/example.py

@@ -38,7 +38,6 @@
 # solver.setParameters(psi=1, Lx=128, Ly=128, Lz=128,shearStrain=1)
 
 solver.initialize(
-    temperature=1.0,
     viscosity=1.0,
     hydrodynamicRadius=1.0,
     includeAngular=False,

include/MobilityInterface/MobilityInterface.h

@@ -25,7 +25,6 @@ enum class periodicity_mode {
 struct Parameters {
   std::vector<real> hydrodynamicRadius;
   real viscosity = 1;
-  real temperature = 0;
   real tolerance = 1e-4; // Tolerance for Lanczos fluctuations
   std::uint64_t seed = 0;
   bool includeAngular = false;
@@ -49,7 +48,6 @@ class Mobility {
   real lanczosTolerance;
   std::shared_ptr<LanczosStochasticVelocities> lanczos;
   std::vector<real> lanczosOutput;
-  real temperature;
   bool includeAngular = false;
   std::mt19937 rng;
 
@@ -98,7 +96,6 @@ class Mobility {
     this->initialized = true;
     this->lanczosSeed = this->rng();
     this->lanczosTolerance = par.tolerance;
-    this->temperature = par.temperature;
     this->includeAngular = par.includeAngular;
   }
 
@@ -120,8 +117,6 @@ class Mobility {
   // will be used automatically
   virtual void sqrtMdotW(device_span<real> ilinear, device_span<real> iangular,
                          real prefactor = 1) {
-    if (this->temperature == 0)
-      return;
     if (prefactor == 0)
       return;
     if (not this->initialized)
@@ -143,9 +138,6 @@ class Mobility {
           "[libMobility] The number of linear velocities does not match the "
           "number of particles");
     }
-    // if (this->needsTorque && linear.size() != angular.size())
-    // throw std::runtime_error("[libMobility] This solver requires angular "
-    //  "velocities when configured with torques");
     const auto numberElements =
         numberParticles + (this->includeAngular ? numberParticles : 0);
     if (not lanczos) {
@@ -178,25 +170,29 @@ class Mobility {
                      thrust::plus<real>());
   }
 
-  // Equivalent to calling Mdot, then stochasticDisplacements, and then thermal
-  // drift. Can be faster in some solvers
-  virtual void hydrodynamicVelocities(device_span<const real> forces,
-                                      device_span<const real> torques,
-                                      device_span<real> linear,
-                                      device_span<real> angular,
-                                      real prefactor = 1) {
+  // computes velocities according to the Langevin equation.
+  virtual void LangevinVelocities(real dt, real kbt,
+                                  device_span<const real> forces,
+                                  device_span<const real> torques,
+                                  device_span<real> linear,
+                                  device_span<real> angular) {
     if (!forces.empty() or !torques.empty()) {
       Mdot(forces, torques, linear, angular);
     }
-    sqrtMdotW(linear, angular, prefactor);
-    thermalDrift(linear, angular, prefactor);
+
+    // prefactors (last argument) are chosen for dimensional consistency
+    if (kbt != 0) {
+      const real sqrt_prefactor = std::sqrt(2 * kbt / dt);
+      sqrtMdotW(linear, angular, sqrt_prefactor);
+      divM(linear, angular, kbt);
+    }
   }
 
-  // Compute the thermal drift,
-  // :math:`k_BT\boldsymbol{\partial}_\boldsymbol{x}\cdot
+  // Compute the divergence of the mobility matrix,
+  // :math:`\boldsymbol{\partial}_\boldsymbol{x}\cdot
   // \boldsymbol{\mathcal{M}}`.
-  virtual void thermalDrift(device_span<real> ilinear,
-                            device_span<real> angular, real prefactor = 1) {
+  virtual void divM(device_span<real> ilinear, device_span<real> angular,
+                    real prefactor = 1) {
     // For open and periodic solvers the thermal drift is zero, so this function
     // does nothing by default.
   }

include/MobilityInterface/pythonify.h

@@ -171,8 +171,6 @@ Initialize the module with a given set of parameters.
 
 Parameters
 ----------
-temperature : float
-		Temperature of the system in energy units (i.e. kT).
 viscosity : float
 		Viscosity of the fluid.
 hydrodynamicRadius : float
@@ -201,7 +199,7 @@ template <class Solver> auto call_sqrtMdotW(Solver &myself, real prefactor) {
 }
 
 const char *sqrtMdotW_docstring = R"pbdoc(
-Computes the stochastic contribution, :math:`\text{prefactor}\sqrt{2T\boldsymbol{\mathcal{M}}}d\boldsymbol{W}`, where :math:`\boldsymbol{\mathcal{M}}` is the grand mobility matrix and :math:`d\boldsymbol{W}` is a Wiener process.
+Computes the stochastic contribution, :math:`{\mathcal{M}}^{1/2} \boldsymbol{W}`, where :math:`\boldsymbol{\mathcal{M}}` is the grand mobility matrix and :math:`\boldsymbol{W}` is a standard normal Gaussian process.
 
 It is required that :py:mod:`setPositions` has been called before calling this function.
 
@@ -255,10 +253,9 @@ array_like
 )pbdoc";
 
 template <class Solver>
-void call_initialize(Solver &myself, real T, real eta, real a,
-                     bool includeAngular, real tol) {
+void call_initialize(Solver &myself, real eta, real a, bool includeAngular,
+                     real tol) {
   libmobility::Parameters par;
-  par.temperature = T;
   par.viscosity = eta;
   par.hydrodynamicRadius = {a};
   par.tolerance = tol;
@@ -278,8 +275,8 @@ template <class Solver> void call_setPositions(Solver &myself, pyarray_c &pos) {
 }
 
 template <class Solver>
-auto call_hydrodynamicVelocities(Solver &myself, pyarray_c &forces,
-                                 pyarray_c &torques, real prefactor) {
+auto call_LangevinVelocities(Solver &myself, real dt, real kbt,
+                             pyarray_c &forces, pyarray_c &torques) {
   auto [f, t, mf, mt] = setup_arrays(myself, forces, torques);
 
   if (forces.size() == 0) // must check because this can be called without
@@ -295,30 +292,30 @@ auto call_hydrodynamicVelocities(Solver &myself, pyarray_c &forces,
 
   auto mf_ptr = cast_to_real(mf);
   auto mt_ptr = cast_to_real(mt);
-  myself.hydrodynamicVelocities(f, t, mf_ptr, mt_ptr, prefactor);
+  myself.LangevinVelocities(dt, kbt, f, t, mf_ptr, mt_ptr);
   int N = myself.getNumberParticles();
   return std::make_pair(mf, mt);
 }
 
-const char *hydrodynamicvelocities_docstring = R"pbdoc(
-Computes the hydrodynamic (deterministic and stochastic) velocities.
+const char *langevinvelocities_docstring = R"pbdoc(
+Computes the hydrodynamic (deterministic and stochastic) velocities according to the Langevin equation,
 
 .. math::
-        \boldsymbol{\mathcal{M}}\begin{bmatrix}\boldsymbol{F}\\\boldsymbol{T}\end{bmatrix} + \text{prefactor}\sqrt{2T\boldsymbol{\mathcal{M}}}d\boldsymbol{W}
+        \boldsymbol{\mathcal{M}}\begin{bmatrix}\boldsymbol{F}\\\boldsymbol{T}\end{bmatrix} + \sqrt{\frac{2k_BT}{\Delta t}}\boldsymbol{\mathcal{M}}^{1/2} \boldsymbol{W} + k_BT \partial_{\boldsymbol{q}}\cdot \boldsymbol{\mathcal{M}}.
 
-If the forces are omitted only the stochastic part is computed.
-If the temperature is zero the stochastic part is omitted.
-Calling this function is equivalent to calling :py:mod:`Mdot` and :py:mod:`sqrtMdotW` in sequence, but in some solvers this can be done more efficiently.
+If forces and torques are omitted then only the stochastic part is computed.
+Calling this function is equivalent to calling :py:mod:`Mdot`, :py:mod:`sqrtMdotW` and :py:mod:`divM` in sequence and applying their respective scalar coefficients, but this can be done more efficiently in a combined fashion in some solvers. By default, this is equvialent to an Euler-Maruyama scheme that uses a random finite difference to compute the divergence term, if necessary, which may not be accurate enough for some applications.
 
 Parameters
 ----------
+dt : float
+		Time step :math:`\Delta t` for the Langevin equation.
+kbt : float
+		Boltzmann constant times temperature, :math:`k_B T`, in units of energy.
 forces : array_like, optional
 		Forces acting on the particles.
 torques : array_like, optional
 		Torques acting on the particles. The solver must have been initialized with includeAngular=True.
-prefactor : float, optional
-		Prefactor to multiply the result by. Default is 1.0.
-
 Returns
 -------
 array_like
@@ -333,7 +330,7 @@ std::unique_ptr<Solver> call_construct(std::string perx, std::string pery,
   return std::make_unique<Solver>(createConfiguration(perx, pery, perz));
 }
 
-template <class Solver> auto call_thermalDrift(Solver &solver, real prefactor) {
+template <class Solver> auto call_divM(Solver &solver, real prefactor) {
   const size_t N = solver.getNumberParticles();
   if (N <= 0) {
     throw std::runtime_error(
@@ -347,12 +344,12 @@ template <class Solver> auto call_thermalDrift(Solver &solver, real prefactor) {
     angular = lp::create_with_framework<real>(last_shape, last_device,
                                               last_framework);
   }
-  solver.thermalDrift(cast_to_real(linear), cast_to_real(angular), prefactor);
+  solver.divM(cast_to_real(linear), cast_to_real(angular), prefactor);
   return std::make_pair(linear, angular);
 }
 
-const char *thermaldrift_docstring = R"pbdoc(
-Computes the thermal drift, :math:`k_BT\boldsymbol{\partial}_\boldsymbol{x}\cdot \boldsymbol{\mathcal{M}}`.
+const char *divM_docstring = R"pbdoc(
+Computes the divergence term, :math:`\boldsymbol{\partial}_\boldsymbol{x}\cdot \boldsymbol{\mathcal{M}}`.
 It is required that :py:mod:`setPositions` has been called before calling this function.
 
 Parameters
@@ -378,8 +375,8 @@ auto define_module_content(
       .def(nb::new_(&call_construct<MODULENAME>), constructor_docstring,
            "periodicityX"_a, "periodicityY"_a, "periodicityZ"_a)
       .def("initialize", call_initialize<MODULENAME>, initialize_docstring,
-           "temperature"_a, "viscosity"_a, "hydrodynamicRadius"_a,
-           "includeAngular"_a = false, "tolerance"_a = 1e-4)
+           "viscosity"_a, "hydrodynamicRadius"_a, "includeAngular"_a = false,
+           "tolerance"_a = 1e-4)
       .def("setPositions", call_setPositions<MODULENAME>,
            "The module will compute the mobility according to this set of "
            "positions.",
@@ -388,11 +385,10 @@ auto define_module_content(
            "forces"_a = pyarray(), "torques"_a = pyarray())
       .def("sqrtMdotW", call_sqrtMdotW<MODULENAME>, sqrtMdotW_docstring,
            "prefactor"_a = 1.0)
-      .def("hydrodynamicVelocities", call_hydrodynamicVelocities<MODULENAME>,
-           hydrodynamicvelocities_docstring, "forces"_a = pyarray_c(),
-           "torques"_a = pyarray_c(), "prefactor"_a = 1)
-      .def("thermalDrift", call_thermalDrift<MODULENAME>,
-           thermaldrift_docstring, "prefactor"_a = 1)
+      .def("LangevinVelocities", call_LangevinVelocities<MODULENAME>,
+           langevinvelocities_docstring, "dt"_a, "kbt"_a,
+           "forces"_a = pyarray_c(), "torques"_a = pyarray_c())
+      .def("divM", call_divM<MODULENAME>, divM_docstring, "prefactor"_a = 1)
       .def("clean", &MODULENAME::clean,
            "Frees any memory allocated by the module.")
       .def_prop_ro_static(

solvers/DPStokes/CMakeLists.txt

@@ -6,6 +6,7 @@ nanobind_add_module(
   STABLE_ABI
   python_wrapper.cu
 )
+uammd_setup_target(${NAME})
 target_link_libraries(${NAME} PRIVATE libMobility_${NAME})
 install(TARGETS libMobility_${NAME} DESTINATION ${CMAKE_INSTALL_PREFIX}/lib)
 install(TARGETS ${NAME} LIBRARY DESTINATION ${Python_SITEARCH}/libMobility)

solvers/DPStokes/mobility.h

@@ -28,7 +28,6 @@ class DPStokes : public libmobility::Mobility {
   Parameters par;
   std::shared_ptr<DPStokesUAMMD> dpstokes;
   DPStokesParameters dppar;
-  real temperature;
   real lanczosTolerance;
   std::string wallmode;
   std::mt19937 rng;
@@ -60,7 +59,6 @@ class DPStokes : public libmobility::Mobility {
 
   void initialize(Parameters ipar) override {
     this->dppar.viscosity = ipar.viscosity;
-    this->temperature = ipar.temperature;
     this->lanczosTolerance = ipar.tolerance;
     this->dppar.mode = this->wallmode;
     this->dppar.hydrodynamicRadius = ipar.hydrodynamicRadius[0];
@@ -146,9 +144,9 @@ class DPStokes : public libmobility::Mobility {
                    this->getNumberParticles(), this->getIncludeAngular());
   }
 
-  void thermalDrift(device_span<real> ilinear, device_span<real> iangular,
-                    real prefactor = 1) override {
-    if (temperature == 0 || prefactor == 0) {
+  void divM(device_span<real> ilinear, device_span<real> iangular,
+            real prefactor = 1) override {
+    if (prefactor == 0) {
       return;
     }
     if (ilinear.size() != 3 * this->numberParticles) {
@@ -182,8 +180,7 @@ class DPStokes : public libmobility::Mobility {
     device_vector thermal_drift_m(ilinear.size(), 0);
     device_vector thermal_drift_d(iangular.size(), 0);
     libmobility::random_finite_differences(mdot, original_pos, thermal_drift_m,
-                                           thermal_drift_d, seed,
-                                           prefactor * temperature);
+                                           thermal_drift_d, seed, prefactor);
     this->setPositions(original_pos);
     thrust::transform(thrust::cuda::par, thermal_drift_m.begin(),
                       thermal_drift_m.end(), linear.begin(), linear.begin(),

solvers/NBody/CMakeLists.txt

@@ -6,6 +6,7 @@ nanobind_add_module(
   STABLE_ABI
   python_wrapper.cu
 )
+uammd_setup_target(${NAME})
 target_link_libraries(${NAME} PRIVATE libMobility_${NAME})
 install(TARGETS libMobility_${NAME} DESTINATION ${CMAKE_INSTALL_PREFIX}/lib)
 install(TARGETS ${NAME} LIBRARY DESTINATION ${Python_SITEARCH}/libMobility)

solvers/NBody/mobility.h

@@ -33,7 +33,6 @@ class NBody : public libmobility::Mobility {
   int Nbatch;
   int NperBatch;
 
-  real temperature;
   std::mt19937 rng;
 
 public:
@@ -105,7 +104,6 @@ class NBody : public libmobility::Mobility {
                                     hydrodynamicRadius * hydrodynamicRadius);
     this->rotMobility = 1.0 / (8 * M_PI * ipar.viscosity * hydrodynamicRadius *
                                hydrodynamicRadius * hydrodynamicRadius);
-    this->temperature = ipar.temperature;
     if (ipar.seed == 0) {
       ipar.seed = std::random_device()();
     }
@@ -159,9 +157,9 @@ class NBody : public libmobility::Mobility {
                                 this->getIncludeAngular(), algorithm, kernel);
   }
 
-  void thermalDrift(device_span<real> ilinear, device_span<real> iangular,
-                    real prefactor = 1) override {
-    if (temperature == 0 || prefactor == 0) {
+  void divM(device_span<real> ilinear, device_span<real> iangular,
+            real prefactor = 1) override {
+    if (prefactor == 0) {
       return;
     }
     if (this->kernel == nbody_rpy::kernel_type::open_rpy)
@@ -198,8 +196,7 @@ class NBody : public libmobility::Mobility {
     device_vector thermal_drift_d(iangular.size(), 0);
 
     libmobility::random_finite_differences(mdot, original_pos, thermal_drift_m,
-                                           thermal_drift_d, seed,
-                                           prefactor * temperature);
+                                           thermal_drift_d, seed, prefactor);
     device_adapter<real> linear(ilinear, device::cuda);
     this->setPositions(original_pos);
     thrust::transform(thrust::cuda::par, thermal_drift_m.begin(),

solvers/PSE/CMakeLists.txt

@@ -6,6 +6,7 @@ nanobind_add_module(
   STABLE_ABI
   python_wrapper.cu
 )
+uammd_setup_target(${NAME})
 target_link_libraries(${NAME} PRIVATE libMobility_${NAME})
 install(TARGETS libMobility_${NAME} DESTINATION ${CMAKE_INSTALL_PREFIX}/lib)
 install(TARGETS ${NAME} LIBRARY DESTINATION ${Python_SITEARCH}/libMobility)