Evaporation dependent perturbed Lennard Jones potential

src/CMakeLists.txt

@@ -5,20 +5,24 @@ set(COMPONENT_NAME azplugins)
 set(_${COMPONENT_NAME}_sources
     ConstantFlow.cc
     export_ImagePotentialBondHarmonic.cc
+    PerturbedLennardJonesEvap.cc
+    export_VariantInterpolated.cc
     module.cc
     ParabolicFlow.cc
     VelocityCompute.cc
     )
 if (ENABLE_HIP)
     list(APPEND _${COMPONENT_NAME}_sources
     export_ImagePotentialBondHarmonicGPU.cc
+    PerturbedLennardJonesEvapGPU.cc
     VelocityComputeGPU.cc
     )
 endif()
 
 # TODO: List all GPU C++ source code files in _${COMPONENT_NAME}_cu_sources.
 set(_${COMPONENT_NAME}_cu_sources
     ImagePotentialBondGPUKernel.cu
+    PerturbedLennardJonesEvapGPU.cu
     VelocityComputeGPU.cu
     VelocityFieldComputeGPU.cu
     )
@@ -32,6 +36,7 @@ set(python_files
     external.py
     flow.py
     pair.py
+    variant.py
     wall.py
     )
 

src/LinearInterpolator1D.h

@@ -0,0 +1,93 @@
+// Copyright (c) 2018-2020, Michael P. Howard
+// Copyright (c) 2021-2026, Auburn University
+// Part of azplugins, released under the BSD 3-Clause License.
+
+#ifndef AZPLUGINS_LINEAR_INTERPOLATOR_1D_H_
+#define AZPLUGINS_LINEAR_INTERPOLATOR_1D_H_
+
+#include <cassert>
+#include <cmath>
+#include <cstdint>
+
+#include "hoomd/HOOMDMath.h"
+
+#if defined(__HIPCC__) || defined(__CUDACC__)
+#define AZPLUGINS_HOSTDEVICE __host__ __device__
+#define AZPLUGINS_FORCEINLINE __forceinline__
+#else
+#define AZPLUGINS_HOSTDEVICE
+#define AZPLUGINS_FORCEINLINE inline
+#endif
+
+namespace hoomd
+    {
+namespace azplugins
+    {
+
+template<typename T> class LinearInterpolator1D
+    {
+    public:
+    AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE LinearInterpolator1D()
+        : m_data(nullptr), m_lo(Scalar(0)), m_hi(Scalar(0)), m_dx(Scalar(0)), m_n(0)
+        {
+        }
+
+    //! Constructor
+    AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE
+    LinearInterpolator1D(const T* data, unsigned int n, Scalar lo, Scalar hi)
+        : m_data(data), m_lo(lo), m_hi(hi), m_n(n)
+        {
+        assert(n >= 2);
+        m_dx = (m_hi - m_lo) / Scalar(n - 1);
+        }
+
+    //! Implement piecewise linear interpolation on a uniform grid
+    AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar operator()(Scalar x) const
+        {
+        const Scalar f = (x - m_lo) / m_dx;
+        int b = static_cast<int>(std::floor(f));
+
+        // If exactly at the top boundary, shift into the last valid cell
+        if (f == Scalar(m_n - 1) && x == m_hi)
+            {
+            --b;
+            }
+
+        assert(b >= 0);
+        assert(b < static_cast<int>(m_n - 1));
+
+        const Scalar frac = f - Scalar(b);
+
+        const Scalar c0 = m_data[b];
+        const Scalar c1 = m_data[b + 1];
+
+        // Linear interpolation
+        return c0 * (Scalar(1) - frac) + c1 * frac;
+        }
+
+    //! Lower bound
+    AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getLo() const
+        {
+        return m_lo;
+        }
+    //! Upper bound
+    AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getHi() const
+        {
+        return m_hi;
+        }
+
+    private:
+    const T* m_data;
+    Scalar m_lo;
+    Scalar m_hi;
+    Scalar m_dx;
+    unsigned int m_n;
+    };
+
+    } // namespace azplugins
+    } // namespace hoomd
+
+#undef AZPLUGINS_HOSTDEVICE
+#undef AZPLUGINS_FORCEINLINE
+
+#endif // AZPLUGINS_LINEAR_INTERPOLATOR_1D_H_

src/LinearInterpolator2D.h

@@ -0,0 +1,131 @@
+// Copyright (c) 2018-2020, Michael P. Howard
+// Copyright (c) 2021-2026, Auburn University
+// Part of azplugins, released under the BSD 3-Clause License.
+
+#ifndef AZPLUGINS_LINEAR_INTERPOLATOR_2D_H_
+#define AZPLUGINS_LINEAR_INTERPOLATOR_2D_H_
+
+#include <cassert>
+#include <cmath>
+#include <cstdint>
+
+#include "hoomd/HOOMDMath.h"
+#include "hoomd/Index1D.h"
+
+#if defined(__HIPCC__) || defined(__CUDACC__)
+#define AZPLUGINS_HOSTDEVICE __host__ __device__
+#define AZPLUGINS_FORCEINLINE __forceinline__
+#else
+#define AZPLUGINS_HOSTDEVICE
+#define AZPLUGINS_FORCEINLINE inline
+#endif
+
+namespace hoomd
+    {
+namespace azplugins
+    {
+
+template<typename T> class LinearInterpolator2D
+    {
+    public:
+    AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE LinearInterpolator2D() : m_data(nullptr), m_indexer()
+        {
+        m_n[0] = 0;
+        m_n[1] = 0;
+        for (int d = 0; d < 2; ++d)
+            {
+            m_lo[d] = Scalar(0);
+            m_hi[d] = Scalar(0);
+            m_dx[d] = Scalar(0);
+            }
+        }
+
+    //! Constructor
+    AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE
+    LinearInterpolator2D(const T* data, const unsigned int* n, const Scalar* lo, const Scalar* hi)
+        : m_data(data), m_indexer(n[0], n[1])
+        {
+        for (int d = 0; d < 2; ++d)
+            {
+            const unsigned int nd = n[d];
+            assert(nd >= 2);
+
+            m_n[d] = nd;
+            m_lo[d] = lo[d];
+            m_hi[d] = hi[d];
+            m_dx[d] = (m_hi[d] - m_lo[d]) / Scalar(nd - 1);
+            }
+        }
+
+    //! Interpolate at (x0, x1)
+    AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar operator()(Scalar x0, Scalar x1) const
+        {
+        const Scalar x[2] = {x0, x1};
+
+        int bin[2];
+        Scalar frac[2];
+
+        for (int d = 0; d < 2; ++d)
+            {
+            const unsigned int nd = m_n[d];
+            const Scalar f = (x[d] - m_lo[d]) / m_dx[d];
+            int b = static_cast<int>(std::floor(f));
+
+            // If exactly at the top boundary, shift into the last valid cell
+            if (f == Scalar(nd - 1) && x[d] == m_hi[d])
+                {
+                --b;
+                }
+
+            assert(b >= 0);
+            assert(b < static_cast<int>(nd - 1));
+
+            bin[d] = b;
+            frac[d] = f - Scalar(b);
+            }
+
+        const unsigned int b0 = static_cast<unsigned int>(bin[0]);
+        const unsigned int b1 = static_cast<unsigned int>(bin[1]);
+        const Scalar c00 = m_data[m_indexer(b0, b1)];
+        const Scalar c01 = m_data[m_indexer(b0, b1 + 1)];
+        const Scalar c10 = m_data[m_indexer(b0 + 1, b1)];
+        const Scalar c11 = m_data[m_indexer(b0 + 1, b1 + 1)];
+
+        // Bilinear interpolation
+        const Scalar t0 = frac[0];
+        const Scalar omt0 = Scalar(1) - t0;
+        const Scalar c0 = c00 * omt0 + c10 * t0;
+        const Scalar c1 = c01 * omt0 + c11 * t0;
+
+        const Scalar t1 = frac[1];
+        return c0 * (Scalar(1) - t1) + c1 * t1;
+        }
+
+    //! Lower bound
+    AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getLo(int dim) const
+        {
+        return m_lo[dim];
+        }
+
+    //! Upper bound
+    AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getHi(int dim) const
+        {
+        return m_hi[dim];
+        }
+
+    private:
+    const T* m_data;
+    Scalar m_lo[2];
+    Scalar m_hi[2];
+    Scalar m_dx[2];
+    unsigned int m_n[2];
+    Index2D m_indexer;
+    };
+
+    } // namespace azplugins
+    } // namespace hoomd
+
+#undef AZPLUGINS_HOSTDEVICE
+#undef AZPLUGINS_FORCEINLINE
+
+#endif // AZPLUGINS_LINEAR_INTERPOLATOR_2D_H_