openforcefield / smarty

Chemical perception tree automated exploration tool.
http://openforcefield.org
MIT License
19 stars 8 forks source link

Non-bonded potentials #221

Open ramess101 opened 7 years ago

ramess101 commented 7 years ago

We anticipate implementing the Mie m-n potential and the Exp-6 potential.

Here are some simple python codes for both potentials:

Note that the Exp-6 code does not account for the fact that at low distances (low values of r) the potential turns over. To correct for this it is common to assign a variable r_max to the minimum value of r at which dU_Exp6/dr = 0. U is assigned to infinity for all values of r less than r_max since these short distances have a very high repulsion.

from __future__ import division #In case using Python 2.7
import numpy as np

def U_Mie(r, epsilon, sigma, n = 12, m = 6):
    """
    The Mie potential calculated in the same units as epsilon. 
    Standard practice is for "epsilon" to actually be epsilon/kB which is [K].
    Note that r and sigma must be in the same units. 
    The exponents (n and m) are set to default values of 12 and 6  to return the LJ potential 
    """
    C = (n/(n-m))*(n/m)**(m/(n-m)) # The normalization constant for the Mie potential
    U = C*epsilon*((r/sigma)**-n - (r/sigma)**-m)
    return U

def U_Exp6(r, epsilon, r_min, alpha):
    """
    The Exp-6 potential calculated in the same units as epsilon. 
    Standard practice is for "epsilon" to actually be epsilon/kB which is [K].
    Note that r and r_min must be in the same units. 
    alpha determines the shape of the Exp-6 potential
    """
    U = epsilon/(1-6./alpha)*((6./alpha)*np.exp(alpha*(1-r/r_min))-(r_min/r)**6)
    return U
bannanc commented 7 years ago

@ramess101 This probably should be added to issue #161