Log Fermi wall potential

log_fermi_constraint.py

import numpy as np
from ase import units
from numba import njit


@njit(fastmath=True)
def log_fermi(positions, radius, temperature, beta):
    eps = 1e-9 # small number to avoid instability
    dists = np.sqrt(np.sum(positions * positions, axis=1))
    exp_term = np.exp(beta * (dists - radius))
    kT = units.kB * temperature 
    E_i = kT * np.log(1 + exp_term)
    E = E_i.sum()
    grad_multiplier = kT * beta * exp_term / (dists * (1 + exp_term) + eps)
    E_grad = grad_multiplier.reshape(-1, 1) * positions
    return E, E_grad


# TODO: Do not compute E and F twice
class LogFermiWallPotential:
    """ Apply logfermi potential for confined molecular dynamics.
    Confines the system to be inside a sphere by applying wall potential.
    
    Method referenced from https://xtb-docs.readthedocs.io/en/latest/xcontrol.html#confining-in-a-cavity
    """
    def __init__(self, radius=5.0, temperature=300, beta=6):
        self.radius = radius
        self.temperature = temperature
        self.beta = beta
        
    def _get_wall_energy_and_force(self, pos):
        E, E_grad = log_fermi(pos, self.radius, self.temperature, self.beta)
        return E, -E_grad
        
    def adjust_forces(self, atoms, forces):
        E_wall, F_wall = self._get_wall_energy_and_force(atoms.positions)
        forces += F_wall
    
    def adjust_potential_energy(self, atoms):
        E_wall, F_wall = self._get_wall_energy_and_force(atoms.positions)
        return E_wall
        
    def adjust_positions(self, atoms, new):
        pass
        
    def get_removed_dof(self, atoms):
        return 0

External spherical wall potential for molecular dynamics in confined area (no PBC), as described in xTB docs.

Formulation

$$
V_\mathrm{wall} = \sum_{i=1}^N k_BT \log \left[1 + e^{\beta ||\mathbf{R}i - \mathbf{O}|| - R\mathrm{sph}}\right]
$$

$\beta$: steepness of wall potential, $\mathbf{O}$: Origin, $R_\mathrm{sph}$: Radius of confining sphere

Usage

Before running dynamics simulation, set the constraint:

atoms = ...
wall_constraint = LogFermiWallPotential()
atoms.set_constraint(wall_constraint)

dyn = ...
dyn.run(1000)

Results

Example system: close O₂ molecules (toy system for giving high repulsion)