From d876f84e8c9df7f5c2d1e2a0ea8a6ff7512f5bb4 Mon Sep 17 00:00:00 2001 From: hechtprojects Date: Wed, 11 Dec 2024 17:25:55 +0100 Subject: [PATCH 1/3] WCA potential added --- amep/utils.py | 97 ++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 96 insertions(+), 1 deletion(-) diff --git a/amep/utils.py b/amep/utils.py index 4415eca..7034fd6 100644 --- a/amep/utils.py +++ b/amep/utils.py @@ -2103,4 +2103,99 @@ def quaternion_multiply(a : np.ndarray, b : np.ndarray): ab[1:]=(a[0] * b[1:] + a[1:] * b[0] + np.cross(a[1:], b[1:])) - return ab \ No newline at end of file + return ab + + + +# ============================================================================= +# INTERACTION POTENTIALS +# ============================================================================= +def wca(r: float | np.ndarray, eps: float = 10.0, sig: float = 1.0): + """ + Weeks-Chandler-Anderson (WCA) potential.[1]_ + + References + ---------- + + .. [1] J. D. Weeks, D. Chandler, and H. C. Andersen, Role of Repulsive + Forces in Determining the Equilibrium Structure of Simple Liquids, + J. Chem. Phys. 54, 5237 (1971). https://doi.org/10.1063/1.1674820 + + Parameters + ---------- + r : float | np.ndarray + Distances. + eps : float, optional + Strength of the potential. The default is 10.0. + sig : float, optional + Effective particle diameter. The default is 1.0. + + Returns + ------- + epot : float | np.ndarray + Potential energy. + + """ + rcut = 2**(1/6)*sig + epot = np.where(r<=rcut, 4*eps*((sig/r)**12-(sig/r)**6)+eps, 0) + return epot + +def dr_wca(r: float | np.ndarray, eps: float = 10.0, sig: float = 1.0): + """ + First derivative of the Weeks-Chandler-Anderson (WCA) potential.[1]_ + + References + ---------- + + .. [1] J. D. Weeks, D. Chandler, and H. C. Andersen, Role of Repulsive + Forces in Determining the Equilibrium Structure of Simple Liquids, + J. Chem. Phys. 54, 5237 (1971). https://doi.org/10.1063/1.1674820 + + Parameters + ---------- + r : float | np.ndarray + Distances. + eps : float, optional + Strength of the potential. The default is 10.0. + sig : float, optional + Effective particle diameter. The default is 1.0. + + Returns + ------- + dr : float | np.ndarray + First derivative. + + """ + rcut = 2**(1/6)*sig + dr = np.where(r<=rcut, 4*eps*(6*sig**6/r**7 - 12*sig**12/r**13), 0) + return dr + +def dr2_wca(r: float | np.ndarray, eps: float = 10.0, sig: float = 1.0): + """ + Second derivative of the Weeks-Chandler-Anderson (WCA) potential.[1]_ + + References + ---------- + + .. [1] J. D. Weeks, D. Chandler, and H. C. Andersen, Role of Repulsive + Forces in Determining the Equilibrium Structure of Simple Liquids, + J. Chem. Phys. 54, 5237 (1971). https://doi.org/10.1063/1.1674820 + + Parameters + ---------- + r : float | np.ndarray + Distances. + eps : float, optional + Strength of the potential. The default is 10.0. + sig : float, optional + Effective particle diameter. The default is 1.0. + + Returns + ------- + dr2 : float | np.ndarray + Second derivative. + + """ + rcut = 2**(1/6)*sig + dr2 = np.where(r<=rcut, 4*eps*(156*sig**12/r**14 - 42*sig**6/r**8), 0) + return dr2 \ No newline at end of file From 06bae455264bf132649cc172af7b0ba21fb9fc63 Mon Sep 17 00:00:00 2001 From: hechtprojects Date: Wed, 11 Dec 2024 18:01:39 +0100 Subject: [PATCH 2/3] single frame temperatures added --- amep/thermo.py | 217 ++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 204 insertions(+), 13 deletions(-) diff --git a/amep/thermo.py b/amep/thermo.py index ee8dfa1..b26da4b 100644 --- a/amep/thermo.py +++ b/amep/thermo.py @@ -31,24 +31,31 @@ # ============================================================================= import numpy as np +from .pbc import distance_matrix, pbc_diff # ============================================================================= -# TEMPERATURE +# KINETIC TEMPERATURE # ============================================================================= def kintemp( - v: np.ndarray, mass: float | None | np.ndarray = None) -> np.ndarray: + v: np.ndarray, m: float | None | np.ndarray = None, + d: int = 2) -> np.ndarray: r''' - Calculates the kinetic temperature per particle. + Calculates the kinetic temperature per particle based on the 2nd + moment of the velocity distribution as described in Ref. [1]_. Notes ----- The kinetic temperature is defined as the kinetic energy. In the overdamped regime (i.e. mass=0), - it is simply :math:`\langle v^2\rangle` (see Ref. [1]_). + it is simply :math:`\langle v^2\rangle` (see Ref. [2]_). References ---------- - .. [1] Caprini, L., & Marini Bettolo Marconi, U. (2020). Active matter at + .. [1] L. Hecht, L. Caprini, H. Löwen, and B. Liebchen, + How to Define Temperature in Active Systems?, J. Chem. Phys. 161, + 224904 (2024). https://doi.org/10.1063/5.0234370 + + .. [2] Caprini, L., & Marini Bettolo Marconi, U. (2020). Active matter at high density: Velocity distribution and kinetic temperature. The Journal of Chemical Physics, 153(18), 184901. https://doi.org/10.1063/5.0029710 @@ -57,8 +64,10 @@ def kintemp( ---------- v : np.ndarray Velocity array of shape (N,3). - mass : float or np.ndarray or None, optional + m : float or np.ndarray or None, optional Particle mass(es). The default is None. + d : int, optional + Spatial dimension. The default is 2. Returns ------- @@ -71,7 +80,7 @@ def kintemp( >>> traj = amep.load.traj("../examples/data/lammps.h5amep") >>> frame = traj[-1] >>> tkin = amep.thermo.kintemp( - ... frame.velocities(), mass=frame.mass() + ... frame.velocities(), m=frame.mass() ... ) >>> fig, axs = amep.plot.new(figsize=(3.6,3)) >>> mp = amep.plot.particles( @@ -91,23 +100,205 @@ def kintemp( :align: center ''' - if mass is None: - mass = 2.0 - + if m is None: + # use version for overdamped particles + m = 2.0 if len(v.shape)==2: # multiple particles vmean = np.mean(v, axis=0) - temp = 0.5 * mass * np.sum((v-vmean)**2, axis=1) + temp = m * np.sum((v-vmean)**2, axis=1) / d elif len(v.shape)==1: # single particle vmean = 0 - temp = 0.5 * mass * np.sum((v-vmean)**2) + temp = m * np.sum((v-vmean)**2) / d else: raise RuntimeError('kintemp: v has the wrong shape. Only shape (N,3) or (3,) is allowed.') - return temp +def Tkin(v: np.ndarray, m: float | np.ndarray, d: int = 2): + ''' + Kinetic temperature based on the second moment of the velocity distribution + and averaged over all particles.[1]_ + + References + ---------- + + .. [1] L. Hecht, L. Caprini, H. Löwen, and B. Liebchen, + How to Define Temperature in Active Systems?, J. Chem. Phys. 161, + 224904 (2024). https://doi.org/10.1063/5.0234370 + + Parameters + ---------- + v : np.ndarray of shape (N,3) + Velocities. + m : float | np.ndarray + Mass(es). + d : int, optional + Spatial dimension. The default is 2. + + Returns + ------- + float + Mean kinetic temperature. + ''' + vmean = np.mean(v, axis=0) + temp = m * np.sum((v-vmean)**2, axis=1) / d + return temp.mean() + +def Tkin4(v: np.ndarray, m: float | np.ndarray, d: int = 2): + ''' + Kinetic temperature based on the 4th moment of the velocity distribution + and averaged over all particles.[1]_ + + References + ---------- + + .. [1] L. Hecht, L. Caprini, H. Löwen, and B. Liebchen, + How to Define Temperature in Active Systems?, J. Chem. Phys. 161, + 224904 (2024). https://doi.org/10.1063/5.0234370 + + Parameters + ---------- + v : np.ndarray of shape (N,3) + Velocities. + m : float | np.ndarray + Mass(es). + d : int, optional + Spatial dimension. The default is 2. + + Returns + ------- + float + Mean kinetic temperature. + ''' + vmean = np.mean(v, axis=0) + v4 = np.sum((v-vmean[None,:])**2, axis=1)**2 + return 0.5 * m * np.sqrt(4 * v4.mean() / (d*(d+2))) + +# ============================================================================= +# CONFIGURATIONAL TEMPERATURE +# ============================================================================= +def Tconf( + coords: np.ndarray, box_boundary: np.ndarray, drU, dr2U, d = 2, + rcut: float = 1.122, pbc: bool = True): + """ + Calculates the configurational temperature averaged over all particles. + + References + ---------- + + .. [1] L. Hecht, L. Caprini, H. Löwen, and B. Liebchen, + How to Define Temperature in Active Systems?, J. Chem. Phys. 161, + 224904 (2024). https://doi.org/10.1063/5.0234370 + + .. [2] S. Saw, L. Costigliola, and J. C. Dyre, Configurational Temperature + in Active Matter. I. Lines of Invariant Physics in the Phase Diagram of + the Ornstein-Uhlenbeck Model, Phys. Rev. E 107, 024609 (2023). + https://doi.org/10.1103/PhysRevE.107.024609 + + .. [3] S. Saw, L. Costigliola, and J. C. Dyre, Configurational Temperature + in Active Matter. II. Quantifying the Deviation from Thermal + Equilibrium, Phys. Rev. E 107, 024610 (2023). + https://doi.org/10.1103/PhysRevE.107.024610 + + Parameters + ---------- + coords : np.ndarray of shape (N,3) + Particle coordinates. + box_boundary : np.ndarray of shape (3,2) + Boundary of the simulation box in the form of + `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`. + drU : function + First derivative of the potential energy function of one particle. + dr2U : function + Second derivative of the potential energy function of one particle. + d : int, optional + Spatial dimension. The current version only works for d=2 and d=3. + The default is 2. + rcut : float, optional + Cutoff radius of the interaction potential. The default is 1.122. + pbc : bool, optional + If True, periodic boundary conditions are applied. The default is True. + + Returns + ------- + temp : float + Average configurational temperature. + nom : float + Nominator. + denom : float + Denominator. + """ + # calculate distances + dmatrix = distance_matrix( + coords, box_boundary, + maxdist = rcut, pbc = pbc + ) + + # get only relevant distances within the cutoff + mask = np.where(dmatrix!=0) + distances = dmatrix[mask] + #print(distances.shape) + + # pairwise difference vectors for relevant pairs + diff = pbc_diff( + coords[mask[0]], coords[mask[1]], box_boundary, pbc = pbc + ) + + # calculate force contribution for each particle + nom = 0.0 + idx_old = 0 + force = 0.0 + for i, idx in enumerate(mask[0]): + + if idx != idx_old: + nom += np.sum(force**2) + force = 0.0 + + force += (drU(distances[i])/distances[i])*diff[i] + idx_old = idx + + # calculate denominator + if d==2: + denom = np.sum(dr2U(distances) + drU(distances)/distances) + elif d==3: + denom = np.sum(dr2U(distances) + 2*drU(distances)/distances) + else: + raise ValueError("Invalid value for d. Use 2 or 3.") + + return nom/denom, nom, denom + + +# ===================================================================== +# OSCILLATOR TEMPERATURE +# ===================================================================== +def Tosc(coords: np.ndarray, k: float): + """ + Oscillator temperature averaged over all particles.[1]_ + + References + ---------- + + .. [1] L. Hecht, L. Caprini, H. Löwen, and B. Liebchen, + How to Define Temperature in Active Systems?, J. Chem. Phys. 161, + 224904 (2024). https://doi.org/10.1063/5.0234370 + + Parameters + ---------- + coords : np.ndarray of shape (N,3) + Particle coordinates. + k : float + Strength of the harmonic potential. + + Returns + ------- + float + Mean oscillator temperature. + """ + return 0.5 * k * np.sum(coords**2, axis=1).mean() + + # ============================================================================= # Kinetic energy # ============================================================================= From ec36307d913689eb3bd928dd6ad447d2bfc2b93c Mon Sep 17 00:00:00 2001 From: hechtprojects Date: Wed, 11 Dec 2024 18:28:50 +0100 Subject: [PATCH 3/3] temperatures added to evaluate module --- amep/evaluate.py | 559 ++++++++++++++++++++++++++++++++++++++++++ test/test_evaluate.py | 4 +- 2 files changed, 560 insertions(+), 3 deletions(-) diff --git a/amep/evaluate.py b/amep/evaluate.py index 12c97f6..70cf752 100644 --- a/amep/evaluate.py +++ b/amep/evaluate.py @@ -5180,3 +5180,562 @@ def indices(self): ''' return self.__indices + + +# ============================================================================= +# TEMPERATURES +# ============================================================================= +class Tkin(BaseEvaluation): + """Kinetic temperature based on the second moment of the velocity + distribution. + """ + def __init__( + self, traj: ParticleTrajectory, skip: float = 0.0, nav: int = 10, + ptype: int | None = None, **kwargs + ) -> None: + r''' + Calculates the kinetic temperature based on the second moment of the + velocity distribution.[1]_ + + References + ---------- + + .. [1] L. Hecht, L. Caprini, H. Löwen, and B. Liebchen, + How to Define Temperature in Active Systems?, J. Chem. Phys. 161, + 224904 (2024). https://doi.org/10.1063/5.0234370 + + Parameters + ---------- + traj : ParticleTrajectory + Trajectory object with simulation data. + skip : float, optional + Skip this fraction at the beginning of the trajectory. + The default is 0.0. + nav : int, optional + Number of frames to consider for the time average. + The default is 10. + ptype : float, optional + Particle type. The default is None. + **kwargs + Other keyword arguments are forwarded to + `amep.thermo.Tkin`. + + ''' + super(Tkin, self).__init__() + + self.name = 'Tkin' + + self.__traj = traj + self.__skip = skip + self.__nav = nav + self.__ptype = ptype + self.__kwargs = kwargs + + self.__frames, self.__avg, self.__indices = average_func( + self.__compute, + self.__traj, + skip = self.__skip, + nr = self.__nav, + indices = True + ) + + self.__times = self.__traj.times[self.__indices] + + + def __compute(self, frame): + r''' + Calculation for a single frame. + + Parameters + ---------- + frame : BaseFrame + A single frame of particle-based simulation data. + + Returns + ------- + temp: float + Mean temperature. + ''' + temp = thermo.Tkin( + frame.velocities(ptype=self.__ptype), + frame.mass(ptype=self.__ptype), + **self.__kwargs + ) + return temp + + + @property + def frames(self): + r''' + Spatial velocity correlation function for each frame. + + Returns + ------- + np.ndarray + Function value for each frame. + + ''' + return self.__frames + + @property + def times(self): + r''' + Times at which the spatial velocity correlation + function is evaluated. + + Returns + ------- + np.ndarray + Times at which the function is evaluated. + + ''' + return self.__times + + @property + def avg(self): + r''' + Time-averaged spatial velocity correlation function + (averaged over the given number of frames). + + Returns + ------- + np.ndarray + Time-averaged spatial velocity correlation function. + + ''' + return self.__avg + + @property + def indices(self): + r''' + Indices of all frames for which the spatial velocity + correlation function has been evaluated. + + Returns + ------- + np.ndarray + Frame indices. + + ''' + return self.__indices + +class Tkin4(BaseEvaluation): + """Kinetic temperature based on the 4th moment of the velocity + distribution. + """ + def __init__( + self, traj: ParticleTrajectory, + skip: float = 0.0, nav: int = 10, + ptype: int | None = None, **kwargs + ) -> None: + r''' + Calculates the kinetic temperature based on the 4th moment of the + velocity distribution.[1]_ + + References + ---------- + + .. [1] L. Hecht, L. Caprini, H. Löwen, and B. Liebchen, + How to Define Temperature in Active Systems?, J. Chem. Phys. 161, + 224904 (2024). https://doi.org/10.1063/5.0234370 + + Parameters + ---------- + traj : ParticleTrajectory + Trajectory object with simulation data. + skip : float, optional + Skip this fraction at the beginning of the trajectory. + The default is 0.0. + nav : int, optional + Number of frames to consider for the time average. + The default is 10. + ptype : float, optional + Particle type. The default is None. + **kwargs + Other keyword arguments are forwarded to + `amep.thermo.Tkin`. + + ''' + super(Tkin4, self).__init__() + + self.name = 'Tkin4' + + self.__traj = traj + self.__skip = skip + self.__nav = nav + self.__ptype = ptype + self.__kwargs = kwargs + + self.__frames, self.__avg, self.__indices = average_func( + self.__compute, + self.__traj, + skip = self.__skip, + nr = self.__nav, + indices = True + ) + + self.__times = self.__traj.times[self.__indices] + + + def __compute(self, frame): + r''' + Calculation for a single frame. + + Parameters + ---------- + frame : BaseFrame + A single frame of particle-based simulation data. + + Returns + ------- + temp: float + Mean temperature. + ''' + temp = thermo.Tkin4( + frame.velocities(ptype=self.__ptype), + frame.mass(ptype=self.__ptype), + **self.__kwargs + ) + return temp + + + @property + def frames(self): + r''' + Spatial velocity correlation function for each frame. + + Returns + ------- + np.ndarray + Function value for each frame. + + ''' + return self.__frames + + @property + def times(self): + r''' + Times at which the spatial velocity correlation + function is evaluated. + + Returns + ------- + np.ndarray + Times at which the function is evaluated. + + ''' + return self.__times + + @property + def avg(self): + r''' + Time-averaged spatial velocity correlation function + (averaged over the given number of frames). + + Returns + ------- + np.ndarray + Time-averaged spatial velocity correlation function. + + ''' + return self.__avg + + @property + def indices(self): + r''' + Indices of all frames for which the spatial velocity + correlation function has been evaluated. + + Returns + ------- + np.ndarray + Frame indices. + + ''' + return self.__indices + +class Tosc(): + """Oscillator temperature. + """ + def __init__( + self, traj: ParticleTrajectory, k: float, + skip: float = 0.0, nav: int = 10, + ptype: int | None = None, + ) -> None: + r''' + Calculates the oscillator temperature.[1]_ + + References + ---------- + + .. [1] L. Hecht, L. Caprini, H. Löwen, and B. Liebchen, + How to Define Temperature in Active Systems?, J. Chem. Phys. 161, + 224904 (2024). https://doi.org/10.1063/5.0234370 + + Parameters + ---------- + traj : ParticleTrajectory + Trajectory object with simulation data. + k : float + Strength of the harmonic confinement. + skip : float, optional + Skip this fraction at the beginning of the trajectory. + The default is 0.0. + nav : int, optional + Number of frames to consider for the time average. + The default is 10. + ptype : float, optional + Particle type. The default is None. + **kwargs + Other keyword arguments are forwarded to + `amep.thermo.Tkin`. + + ''' + super(Tkin4, self).__init__() + + self.name = 'Tosc' + + self.__traj = traj + self.__skip = skip + self.__nav = nav + self.__ptype = ptype + self.__k = k + + self.__frames, self.__avg, self.__indices = average_func( + self.__compute, + self.__traj, + skip = self.__skip, + nr = self.__nav, + indices = True + ) + + self.__times = self.__traj.times[self.__indices] + + + def __compute(self, frame): + r''' + Calculation for a single frame. + + Parameters + ---------- + frame : BaseFrame + A single frame of particle-based simulation data. + + Returns + ------- + temp: float + Mean temperature. + ''' + temp = thermo.Tosc( + frame.coords(ptype=self.__ptype), + self.__k + ) + return temp + + @property + def frames(self): + r''' + Spatial velocity correlation function for each frame. + + Returns + ------- + np.ndarray + Function value for each frame. + + ''' + return self.__frames + + @property + def times(self): + r''' + Times at which the spatial velocity correlation + function is evaluated. + + Returns + ------- + np.ndarray + Times at which the function is evaluated. + + ''' + return self.__times + + @property + def avg(self): + r''' + Time-averaged spatial velocity correlation function + (averaged over the given number of frames). + + Returns + ------- + np.ndarray + Time-averaged spatial velocity correlation function. + + ''' + return self.__avg + + @property + def indices(self): + r''' + Indices of all frames for which the spatial velocity + correlation function has been evaluated. + + Returns + ------- + np.ndarray + Frame indices. + + ''' + return self.__indices + +class Tconf(BaseEvaluation): + """Configurational temperature. + """ + def __init__( + self, traj: ParticleTrajectory, drU, dr2U, + skip: float = 0.0, nav: int = 10, + ptype: int | None = None, **kwargs + ) -> None: + r''' + Calculates the configurational temperature. + + References + ---------- + + .. [1] L. Hecht, L. Caprini, H. Löwen, and B. Liebchen, + How to Define Temperature in Active Systems?, J. Chem. Phys. 161, + 224904 (2024). https://doi.org/10.1063/5.0234370 + + .. [2] S. Saw, L. Costigliola, and J. C. Dyre, Configurational Temperature + in Active Matter. I. Lines of Invariant Physics in the Phase Diagram of + the Ornstein-Uhlenbeck Model, Phys. Rev. E 107, 024609 (2023). + https://doi.org/10.1103/PhysRevE.107.024609 + + .. [3] S. Saw, L. Costigliola, and J. C. Dyre, Configurational Temperature + in Active Matter. II. Quantifying the Deviation from Thermal + Equilibrium, Phys. Rev. E 107, 024610 (2023). + https://doi.org/10.1103/PhysRevE.107.024610 + + Parameters + ---------- + traj : ParticleTrajectory + Trajectory object with simulation data. + drU : function + First derivative of the potential energy function of one particle. + For example, one can use amep.utils.dr_wca. + dr2U : function + Second derivative of the potential energy function of one particle. + For example, one can use amep.utils.dr2_wca. + skip : float, optional + Skip this fraction at the beginning of the trajectory. + The default is 0.0. + nav : int, optional + Number of frames to consider for the time average. + The default is 10. + ptype : float, optional + Particle type. The default is None. + **kwargs + Other keyword arguments are forwarded to + `amep.thermo.Tkin`. + + ''' + super(Tconf, self).__init__() + + self.name = 'Tconf' + + self.__traj = traj + self.__skip = skip + self.__nav = nav + self.__ptype = ptype + self.__kwargs = kwargs + self.__drU = drU + self.__dr2U = dr2U + + self.__frames, self.__avg, self.__indices = average_func( + self.__compute, + self.__traj, + skip = self.__skip, + nr = self.__nav, + indices = True + ) + + self.__times = self.__traj.times[self.__indices] + + + def __compute(self, frame): + r''' + Calculation for a single frame. + + Parameters + ---------- + frame : BaseFrame + A single frame of particle-based simulation data. + + Returns + ------- + temp: float + Mean temperature. + ''' + temp = thermo.Tconf( + frame.coords(ptype=self.__ptype), + frame.box, + self.__drU, + self.__dr2U, + **self.__kwargs + ) + return temp + + + @property + def frames(self): + r''' + Spatial velocity correlation function for each frame. + + Returns + ------- + np.ndarray + Function value for each frame. + + ''' + return self.__frames + + @property + def times(self): + r''' + Times at which the spatial velocity correlation + function is evaluated. + + Returns + ------- + np.ndarray + Times at which the function is evaluated. + + ''' + return self.__times + + @property + def avg(self): + r''' + Time-averaged spatial velocity correlation function + (averaged over the given number of frames). + + Returns + ------- + np.ndarray + Time-averaged spatial velocity correlation function. + + ''' + return self.__avg + + @property + def indices(self): + r''' + Indices of all frames for which the spatial velocity + correlation function has been evaluated. + + Returns + ------- + np.ndarray + Frame indices. + + ''' + return self.__indices diff --git a/test/test_evaluate.py b/test/test_evaluate.py index 409d9d9..66c400f 100644 --- a/test/test_evaluate.py +++ b/test/test_evaluate.py @@ -17,9 +17,7 @@ # # Contact: Lukas Hecht (lukas.hecht@pkm.tu-darmstadt.de) # ============================================================================= -"""Test units for the continuum class. - -Including it's main class, readers and methods.""" +"""Test units for the evaluate module.""" import unittest from pathlib import Path from matplotlib import use