Tutorial 05: Simlammps wrapper

In this tutorial we will go through a simple example of how to use the wrapper for the LAMMPS simulation engine. You can find the wrapper here.

Background

A wrapper for LAMMPS has been present in SimPhoNy since its initial version, and it is the first simmulation engine we supported in version 3.

This wrapper is a good example of the 3-layered-design where the Syntactic layer is a third party library. In this case we use PyLammps, a Python binding for LAMMPS created and maintained by the LAMMPS developers.

Let’s get hands on

Installation

We will start by quickly going through the installation of this tool. Like we explain in the wrapper development section, there are two options:

  • Using Docker: run ./docker_install.sh.
  • Local installation: remember that you must have a compatible version of osp-core installed.

Install the ontology via pico install ontology.simlammps.yml.

Run ./install_engine.sh.

  • Note that you will be asked for a superuser password to install required libraries for the installation (make, libjpeg, libpng…)
  • Currently we support Ubuntu and centos.

Install the wrapper by running python setup.py install.

That should be all needed to use simlammps!

Simple example

This is an adaptation of simlammps/examples/simple.py. As usual, we start importing the necessary components:

[1]:
from osp.core import simlammps_ontology
from osp.wrappers.simlammps import SimlammpsSession

We create the wrapper instance. All wrappers are created by defining their own session class.

There is no need to specify a syntactic layer (PyLammps). The session will generate one.

[2]:
simlammps_session = SimlammpsSession()
simlammps = simlammps_ontology.SimlammpsWrapper(session=simlammps_session)

LAMMPS output is captured by PyLammps wrapper

Next, we can define some necessary settings for the run:

[3]:
# Define the simulation box
box = simlammps_ontology.SimulationBox()
face_x = simlammps_ontology.FaceX(vector=(10, 0, 0))
face_x.add(simlammps_ontology.Periodic())
face_y = simlammps_ontology.FaceY(vector=(0, 10, 0))
face_y.add(simlammps_ontology.Periodic())
face_z = simlammps_ontology.FaceZ(vector=(0, 0, 10))
face_z.add(simlammps_ontology.Periodic())
box.add(face_x, face_y, face_z)
simlammps.add(box)

# molecular dynamics model
md_nve = simlammps_ontology.MolecularDynamics()
simlammps.add(md_nve)

# solver component:
sp = simlammps_ontology.SolverParameter()

# integration time:
steps = 100
itime = simlammps_ontology.IntegrationTime(steps=steps)

sp.add(itime)
verlet = simlammps_ontology.Verlet()

sp.add(verlet)
simlammps.add(sp)

# Mass and material for the atoms
mass = simlammps_ontology.Mass(value=0.2)
material = simlammps_ontology.Material()

material.add(mass)
simlammps.add(material)

# Interatomic force as material relation
lj = simlammps_ontology.LennardJones_6_12(cutoff_distance=2.5,
                                          energy_well_depth=1.0,
                                          van_der_waals_radius=1.0)
lj.add(material)
simlammps.add(lj)
[3]:
<SIMLAMMPS_ONTOLOGY.LENNARD_JONES_6_12: 9c9c1672-a9f2-4eba-85a4-060b56addf2a,  SimlammpsSession: @0x7fac60dc7a90>

Now we add some atoms:

[4]:
particle = simlammps_ontology.Atom()
particle.add(material,
             simlammps_ontology.Position(vector=(1, 6, 3)),
             simlammps_ontology.Velocity(vector=(1, 0, 0)))
simlammps.add(particle)

particle = simlammps_ontology.Atom()
particle.add(material,
             simlammps_ontology.Position(vector=(2, 1, 4)),
             simlammps_ontology.Velocity(vector=(2, 0, 2)))
simlammps.add(particle)

particle = simlammps_ontology.Atom()
# The velocity is not required (the position is)
particle.add(material,
             simlammps_ontology.Position(vector=(7, 3, 0)))
simlammps.add(particle)
[4]:
<SIMLAMMPS_ONTOLOGY.ATOM: db96a76f-4d70-4e19-b460-46ee286f831e,  SimlammpsSession: @0x7fac60dc7a90>

To run the simulation, we call the run() method of the session. The run method sends the information to the engine, and tells it to run the number of steps defined in the Integration Time entity (100):

[5]:
simlammps.session.run()

Since we will run the simulation a couple of times, we can define a simple function for showing the position and velocities of the atoms:

[6]:
def print_info():
    for atom in simlammps.iter(oclass=simlammps_ontology.Atom):
        # Remember that Cuds.get(oclass) returns a list
        # We now all atoms have one (and only one) position
        position = atom.get(oclass=simlammps_ontology.Position)[0]
        # But the atoms might not have a velocity
        velocity = atom.get(oclass=simlammps_ontology.Velocity)
        print("Atom " + str(atom.uid) + ":")
        print(" - Position: " + str(position.vector))
        if velocity:
            print(" - Velocity: " + str(velocity[0].vector))

Now we can easily print the results of the run:

[7]:
print_info()
Atom fd4199d4-4d1a-425c-8010-60efca65bd1c:
 - Position: [1.5 6.  3. ]
 - Velocity: [1. 0. 0.]
Atom f9a32d14-b638-4796-9407-4b1ae6be43cb:
 - Position: [3. 1. 5.]
 - Velocity: [2. 0. 2.]
Atom db96a76f-4d70-4e19-b460-46ee286f831e:
 - Position: [7. 3. 0.]

Finally, let’s change the velocities and run again, but now for 200 steps:

[8]:
from random import randint

for atom in simlammps.iter(oclass=simlammps_ontology.Atom):
    # But the atoms might not have a velocity
    velocity = atom.get(oclass=simlammps_ontology.Velocity)
    if velocity:
        velocity[0].vector = (randint(-3, 3), randint(-3, 3), randint(-3, 3))
    else:
        atom.add(simlammps_ontology.Velocity(vector = (randint(-3, 3), randint(-3, 3), randint(-3, 3))))

solver_parameter = simlammps.get(oclass=simlammps_ontology.SolverParameter)[0]
integration_time = solver_parameter.get(oclass=simlammps_ontology.IntegrationTime)[0]
integration_time.steps = 200

simlammps.session.run()
print_info()
Atom fd4199d4-4d1a-425c-8010-60efca65bd1c:
 - Position: [0.5 4.  6. ]
 - Velocity: [-1. -2.  3.]
Atom f9a32d14-b638-4796-9407-4b1ae6be43cb:
 - Position: [6. 8. 8.]
 - Velocity: [ 3. -3.  3.]
Atom db96a76f-4d70-4e19-b460-46ee286f831e:
 - Position: [9. 4. 0.]
 - Velocity: [2. 1. 0.]
[ ]: