mosdef-hub / gmso

Flexible storage of chemical topology for molecular simulation
https://gmso.mosdef.org
MIT License
52 stars 33 forks source link

Energy evaluations #152

Open mattwthompson opened 5 years ago

mattwthompson commented 5 years ago

I was initially resistant to bringing this up, because it's a big task to compute the potential energy of the entire system. But I wonder if it's tractable to do it for individual potentials.

My first idea is that we implement a Potential.evaluate function in a way that it works for all the potentials that inherit from this class. Usage would be something like this

import numpy as np
from topology.core.potential import Potential

mypotential = MyFavoritePotential(
    expression='ax+b'
    parameters={'a': 0.5: 'b': 1.0}
    independent_variables={'x'}
)

x = np.linspace(0, 1, 101)
f = mypotential.evaluate(substitution_variables={'x': x})

from which you could plot f(x). This seems like something sympy should be able to do but it would need to be implemented in a way that is completely agnostic to the particulars of an individual potential.

ahy3nz commented 5 years ago

https://docs.sympy.org/latest/modules/numeric-computation.html#so-which-should-i-use

ahy3nz commented 5 years ago

I haven't done much thinking on this particular matter - for hacking together to look at I-J lennard jones potentials and doing comparisons I've just been making top.AtomType objects, parsing their expression and parameters, and then my own sympy calls to evaluate them (i've been using evalf or subs), but it looks like sympy has advice for how to efficiently conduct large, vectorized calculations on sympy expressions

mattwthompson commented 5 years ago

Can't get it to work with unyt at the moment but here's the basics:

>> expr = sympy.sympify('sigma*epsilon')
>>> params = {sympy.symbols('sigma'): 1, sympy.symbols('epsilon'): 2}
>>> expr.evalf(subs=x) #@{**x, **y})
2.00000000000000

>>> LJ = Potential(expression='sigma*epsilon', parameters={'sigma': 1.0*u.m, 'epsilon': 2.0*u.m}, independent_variables={'x'})
>>> LJ.expression.evalf(subs=params) #{**x, **y}) # Note: using `params`, not `LJ.parameters`