Source code for pyretis.forcefield.potentials.pairpotentials.pairpotential

# -*- coding: utf-8 -*-
# Copyright (c) 2021, PyRETIS Development Team.
# Distributed under the LGPLv2.1+ License. See LICENSE for more info.
"""This module defines some helper functions for pair potentials.

Important methods defined here
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

mixing_parameters (:py:func:`.mixing_parameters`)
    Definition of mixing rules which are used in some force fields
    for generating cross interactions.

generate_pair_interactions (:py:func:`generate_pair_interactions`)
    Function to generate pair parameters from atom parameters.
"""
import itertools
import logging
import numpy as np
logger = logging.getLogger(__name__)  # pylint: disable=invalid-name
logger.addHandler(logging.NullHandler())


__all__ = ['mixing_parameters', 'generate_pair_interactions']


[docs]def _check_pair_parameters(parameters): """Check that the required pair parameters are given. If the parameters are not given, we will just set them to default values. Parameters ---------- parameters : dict The settings for the potential. Note ---- This function will not return anything, but it will update the input `parameters` with default values. """ for pair in parameters: for key in ('epsilon', 'sigma'): if key not in parameters[pair]: msg = '{} for {} not given. Set to 0.0'.format(key, pair) logger.warning(msg) parameters[pair][key] = 0.0 if 'rcut' not in parameters[pair]: factor = parameters[pair].get('factor', 0.0) rcut = factor * parameters[pair]['sigma'] parameters[pair]['rcut'] = rcut msg = ('rcut for {} not given. Using factor to set it to: ' '{} * {} = {}') msg = msg.format(pair, parameters[pair]['sigma'], factor, rcut) logger.info(msg)
[docs]def generate_pair_interactions(parameters, mixing): """Generate pair parameters from atom parameters. The parameters are given as a dictionary where the keys are either just integers -- which defines atom parameters -- or tuples which define pair interactions. Parameters ---------- parameters : dict This dict contain the atom parameters. mixing : string Determines how we should mix pair interactions. """ _check_pair_parameters(parameters) atoms = [] pair_param = {} for key in parameters: if key == 'mixing' or isinstance(key, tuple): continue atoms.append(key) for (atmi, atmj) in itertools.product(atoms, atoms): pari = parameters[atmi] parj = parameters[atmj] if (atmi, atmj) in pair_param: continue if (atmi, atmj) in parameters: pair_param[atmi, atmj] = dict(parameters[atmi, atmj]) pair_param[atmj, atmi] = pair_param[atmi, atmj] continue elif (atmj, atmi) in parameters: pair_param[atmj, atmi] = dict(parameters[atmj, atmi]) pair_param[atmi, atmj] = pair_param[atmj, atmi] continue if atmi == atmj: eps_ij = pari['epsilon'] sig_ij = pari['sigma'] rcut_ij = pari['rcut'] pair_param[atmi, atmi] = {'epsilon': eps_ij, 'sigma': sig_ij, 'rcut': rcut_ij} else: eps_ij, sig_ij, rcut_ij = mixing_parameters(pari['epsilon'], pari['sigma'], pari['rcut'], parj['epsilon'], parj['sigma'], parj['rcut'], mixing) pair_param[atmi, atmj] = {'epsilon': eps_ij, 'sigma': sig_ij, 'rcut': rcut_ij} pair_param[atmj, atmi] = pair_param[atmi, atmj] # if (atmj, atmi) not in pair_param: # pair_param[atmj, atmi] = pair_param[atmi, atmj] return pair_param
[docs]def mixing_parameters(epsilon_i, sigma_i, rcut_i, epsilon_j, sigma_j, rcut_j, mixing='geometric'): r"""Define the so-called mixing rules. These mixing rules are used for some force fields when generating cross interactions. The known mixing rules are: 1. Geometric: * .. math:: \epsilon_{ij} = \sqrt{\epsilon_{i} \times \epsilon_{j}} * .. math:: \sigma_{ij} = \sqrt{\sigma_{i} \times \sigma_{j}} * .. math:: r_{\text{c},ij} = \sqrt{r_{\text{c},i} \times r_{\text{c},j}} 2. Arithmetic: * .. math:: \epsilon_{ij} = \sqrt{\epsilon_{i} \times \epsilon_{j}} * .. math:: \sigma_{ij} = \frac{\sigma_{i} \times \sigma_{j}}{2} * .. math:: r_{\text{c},ij} = \frac{r_{\text{c},i} \times r_{\text{c},j}}{2} 3. Sixthpower * .. math:: \epsilon_{ij} = 2 \sqrt{\epsilon_{i} \times \epsilon_{j}} \frac{\sigma_i^3 \times \sigma_j^3}{\sigma_i^6 + \sigma_j^6} * .. math:: \sigma_{ij} = \left( \frac{\sigma_{i}^6 \times \sigma_{j}^6}{2} \right)^{1/6} * .. math:: r_{\text{c},ij} = \left(\frac{r_{\text{c},i}^6 \times r_{\text{c},j}^6}{2}\right)^{1/6} Parameters ---------- epsilon_i : float Lennard-Jones epsilon parameter for a particle of type `i`. sigma_i : float Lennard-Jones sigma parameter for a particle of type `i`. rcut_i : float Lennard-Jones cut-off value for a particle of type `i`. epsilon_j : float Lennard-Jones epsilon parameter for a particle of type `j`. sigma_j : float Lennard-Jones sigma parameter for a particle of type `j`. rcut_j : float Lennard-Jones cut-off value for a particle of type `j`. mixing : string, optional Represents what kind of mixing that should be done. Returns ------- out[0] : float The mixed ``epsilon_ij`` parameter. out[1] : float The mixed ``sigma_ij`` parameter. out[2] : float The mixed ``rcut_ij`` parameter. """ if mixing == 'geometric': epsilon_ij = np.sqrt(epsilon_i * epsilon_j) sigma_ij = np.sqrt(sigma_i * sigma_j) rcut_ij = np.sqrt(rcut_i * rcut_j) elif mixing == 'arithmetic': epsilon_ij = np.sqrt(epsilon_i * epsilon_j) sigma_ij = 0.5 * (sigma_i + sigma_j) rcut_ij = 0.5 * (rcut_i + rcut_j) elif mixing == 'sixthpower': si3 = sigma_i**3 si6 = si3**2 sj3 = sigma_j**3 sj6 = sj3**2 avgs6 = 0.5 * (si6 + sj6) epsilon_ij = np.sqrt(epsilon_i * epsilon_j) * si3 * sj3 / avgs6 sigma_ij = avgs6**(1.0 / 6.0) rcut_ij = (0.5 * (rcut_i**6 + rcut_j**6))**(1.0 / 6.0) else: epsilon_ij = 0.0 sigma_ij = 0.0 rcut_ij = 0.0 logger.warning('No mixing requested. Cross interactions set to 0') return epsilon_ij, sigma_ij, rcut_ij