Merge pull request #482 from akohlmey/add-pair-python

examples of Python code you can write for both pure Python operations

and callbacks to LAMMPS. See "fix python"_fix_python.html to learn about

possibilities to execute Python code during each time step.

Through the "python pair style"_pair_python.html it is also possible

to define potential functions as python code.

To run pure Python code from LAMMPS, you only need to build LAMMPS

with the PYTHON package installed:

fix_press_berendsen.html

fix_print.html

fix_property_atom.html

fix_python.html

fix_qbmsst.html

fix_qeq.html

fix_qeq_comb.html

pair_oxdna2.html

pair_peri.html

pair_polymorphic.html

pair_python.html

pair_quip.html

pair_reax.html

pair_reaxc.html

"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c

:link(lws,http://lammps.sandia.gov)

:link(ld,Manual.html)

:link(lc,Section_commands.html#comm)

:line

pair_style python command :h3

[Syntax:]

pair_style python cutoff :pre

cutoff = global cutoff for interactions in python potential classes

[Examples:]

pair_style python 2.5

pair_coeff * * py_pot.LJCutMelt lj :pre

pair_style hybrid/overlay coul/long 12.0 python 12.0

pair_coeff * * coul/long

pair_coeff * * python py_pot.LJCutSPCE OW NULL :pre

[Description:]

The {python} pair style provides a way to define pairwise additive

potential functions as python script code that is loaded into LAMMPS

from a python file which must contain specific python class definitions.

This allows to rapidly evaluate different potential functions without

having to modify and recompile LAMMPS. Due to python being an

interpreted language, however, the performance of this pair style is

going to be significantly slower (often between 20x and 100x) than

corresponding compiled code. This penalty can be significantly reduced

through generating tabulations from the python code through the

"pair_write"_pair_write.html command, which is supported by this style.

Only a single pair_coeff command is used with the {python} pair style

which specifies a python class inside a python module or file that

LAMMPS will look up in the current directory, the folder pointed to by

the LAMMPS_POTENTIALS environment variable or somewhere in your python

path. A single python module can hold multiple python pair class

definitions. The class definitions itself have to follow specific rules

that are explained below.

Atom types in the python class are specified through symbolic constants,

typically strings. These are mapped to LAMMPS atom types by specifying

N additional arguments after the class name in the pair_coeff command,

where N must be the number of currently defined atom types:

As an example, imagine a file {py_pot.py} has a python potential class

names {LJCutMelt} with parameters and potential functions for a two

Lennard-Jones atom types labeled as 'LJ1' and 'LJ2'. In your LAMMPS

input and you would have defined 3 atom types, out of which the first

two are supposed to be using the 'LJ1' parameters and the third

the 'LJ2' parameters, then you would use the following pair_coeff

command:

pair_coeff * * py_pot.LJCutMelt LJ1 LJ1 LJ2 :pre

The first two arguments [must] be * * so as to span all LAMMPS atom types.

The first two LJ1 arguments map LAMMPS atom types 1 and 2 to the LJ1

atom type in the LJCutMelt class of the py_pot.py file.  The final LJ2

argument maps LAMMPS atom type 3 to the LJ2 atom type the python file.

If a mapping value is specified as NULL, the mapping is not performed,

any pair interaction with this atom type will be skipped. This can be

used when a {python} potential is used as part of the {hybrid} or

{hybrid/overlay} pair style. The NULL values are then placeholders for

atom types that will be used with other potentials.

:line

The python potential file has to start with the following code:

from __future__ import print_function

class LAMMPSPairPotential(object):

    def __init__(self):

        self.pmap=dict()

        self.units='lj'

    def map_coeff(self,name,ltype):

        self.pmap\[ltype\]=name

    def check_units(self,units):

        if (units != self.units):

           raise Exception("Conflicting units: %s vs. %s" % (self.units,units))

:pre

Any classes with definitions of specific potentials have to be derived

from this class and should be initialize in a similar fashion to the

example given below. NOTE: The class constructor has to set up a data

structure containing the potential parameters supported by this class.

It should also define a variable {self.units} containing a string

matching one of the options of LAMMPS' "units"_units.html command, which

is used to verify, that the potential definition in the python class and

in the LAMMPS input match. Example for a single type Lennard-Jones

potential class {LJCutMelt} in reducted units, which defines an atom

type {lj} for which the parameters epsilon and sigma are both 1.0:

class LJCutMelt(LAMMPSPairPotential):

    def __init__(self):

        super(LJCutMelt,self).__init__()

        # set coeffs: 48*eps*sig**12, 24*eps*sig**6,

        #              4*eps*sig**12,  4*eps*sig**6

        self.units = 'lj'

        self.coeff = {'lj'  : {'lj'  : (48.0,24.0,4.0,4.0)}}

:pre

The class also has to provide two methods for the computation of the

potential energy and forces, which have be named {compute_force},

and {compute_energy}, which both take 3 numerical arguments:

  rsq   = the square of the distance between a pair of atoms (float) :li

  itype = the (numerical) type of the first atom :li

  jtype = the (numerical) type of the second atom :ul

This functions need to compute the force and the energy, respectively,

and use the result as return value. The functions need to use the

{pmap} dictionary to convert the LAMMPS atom type number to the symbolic

value of the internal potential parameter data structure. Following

the {LJCutMelt} example, here are the two functions:

   def compute_force(self,rsq,itype,jtype):

        coeff = self.coeff\[self.pmap\[itype\]\]\[self.pmap\[jtype\]\]

        r2inv  = 1.0/rsq

        r6inv  = r2inv*r2inv*r2inv

        lj1 = coeff\[0\]

        lj2 = coeff\[1\]

        return (r6inv * (lj1*r6inv - lj2))*r2inv :pre

    def compute_energy(self,rsq,itype,jtype):

        coeff = self.coeff\[self.pmap\[itype\]\]\[self.pmap\[jtype\]\]

        r2inv  = 1.0/rsq

        r6inv  = r2inv*r2inv*r2inv

        lj3 = coeff\[2\]

        lj4 = coeff\[3\]

        return (r6inv * (lj3*r6inv - lj4)) :pre

IMPORTANT NOTE: for consistency with the C++ pair styles in LAMMPS,

the {compute_force} function follows the conventions of the Pair::single()

methods and does not return the full force, but the force scaled by

the distance between the two atoms, so this value only needs to be

multiplied by delta x, delta y, and delta z to conveniently obtain

the three components of the force vector between these two atoms.

:line

IMPORTANT NOTE: The evaluation of scripted python code will slow down

the computation pair-wise interactions quite significantly. However,

this can be largely worked around through using the python pair style

not for the actual simulation, but to generate tabulated potentials

on the fly using the "pair_write"_pair_write.html command. Please

see below for an example LAMMPS input of how to build a table file:

pair_style python 2.5

pair_coeff * * py_pot.LJCutMelt lj

shell rm -f melt.table

pair_write  1 1 2000 rsq 0.01 2.5 lj1_lj2.table lj :pre

Note, that it is strong recommended to try to [delete] the potential

table file before generating it. Since the {pair_write} command will

always append to a table file, which pair style table will use the first

match. Thus when changing the potential function in the python class,

the table pair style will still read the old variant.

After switching the pair style to {table}, the potential tables need

to be assigned to the LAMMPS atom types like this:

pair_style      table linear 2000

pair_coeff      1  1  melt.table lj :pre

This can also be done for more complex systems. Please see the

{examples/python} folders for a few more examples.

:line

[Mixing, shift, table, tail correction, restart, rRESPA info]:

Mixing of potential parameters has to be handled inside the provided

python module. The python pair style simply assumes that force and

energy computation can be correctly performed for all pairs of atom

types as they are mapped to the atom type labels inside the python

potential class.

This pair style does not support the "pair_modify"_pair_modify.html

shift, table, and tail options.

This pair style does not write its information to "binary restart

files"_restart.html, since it is stored in potential files.  Thus, you

need to re-specify the pair_style and pair_coeff commands in an input

script that reads a restart file.

This pair style can only be used via the {pair} keyword of the

"run_style respa"_run_style.html command.  It does not support the

{inner}, {middle}, {outer} keywords.

:line

[Restrictions:]

This pair style is part of the PYTHON package.  It is only enabled

if LAMMPS was built with that package.  See

the "Making LAMMPS"_Section_start.html#start_3 section for more info.

[Related commands:]

"pair_coeff"_pair_coeff.html, "pair_write"_pair_write.html,

"pair style table"_pair_table.html

[Default:] none

   pair_oxdna2

   pair_peri

   pair_polymorphic

   pair_python

   pair_quip

   pair_reax

   pair_reaxc

This folder contains several LAMMPS input scripts and a python module

file py_pot.py to demonstrate the use of the pair style python.

in.pair_python_melt:

This is a version of the melt example using the python pair style.  The first

part of the output should have identical energies, temperature and pressure

than the melt example.  The following two sections then demonstrate how to

restart with pair style python from a restart file and a data file.

in.pair_python_hybrid:

This versions shows how to mix regular pair styles with a python pair style.

However, in this case both potentials are the same, so the energies and

pressure in the output should be identical to that of the previous example.

in.pair_python_spce:

This input shows a simulation of small bulk water system with the SPC/E

water potential. Since the python pair style does not support computing

coulomb contributions, pair style hybrid/overload is used to combine

the python style containing the Lennard-Jones part with the long-range coulomb.

Same as for the previous example, it also showcases restarting.

in.pair_python_table:

This input demonstrates the use of using the python pair style to build

a table file for use with pair style table. This will run much faster

than the python pair style. This example tabulates the melt example from

above. Note that tabulation is approximative, so the output will only

agree with the melt result to some degree. 

in.pair_python_coulomb:

This is another tabulation example, this time for the SPC/E water example

with cutoff coulomb interactions.

Please note, that tabulating long-range coulomb has a systematic error in

forces and energies for all systems with bonds, angle and dihedrals.

In this case, this will only affect the energies, since the water molecules

are held rigid with fix shake. To enable long-range coulomb the coul/cut

style needs to be replaced with coul/long, a suitable kspace style added

and the pppm keyword added to the table pair style definition.

in.pair_python_long:

The final example shows how to combine long-range coulomb with tabulation

for only the short range interactions via pair style hybrid/overlay.