Toolkit for multiscale modeling with LAMMPS in Python 3.x. It integrates with LAMMPS-GUI and leverages Large Language Models for rapid prototyping.
- 📄 Read this file in
HTML
(includes 🎬 videos and enhanced visuals). - 📚 Online Documentation - 🐍Python *(including how to 🛠️ install,
▶️ run and ✏️ modify Pizza3) - 📚 Online Documentation - 📊Matlab/Octave (advanced 📈📉 post-treatment)
- 📐POST Examples and Analysis Workflow (examples of 🌊 fluid-structure interactions)
- 🔖 Releases (✨ new capabilities are detailed there, ⬇️ download page)
[TOC]
Pizza3 began as a fork and extension of the Pizza.py toolkit for LAMMPS, rewritten in Python 3.x to modernize the original Python 2.x framework. Over time, it has evolved into a flexible and loosely integrated toolkit designed for LAMMPS users, regardless of the simulation style. Pizza3 empowers users to facilitate multiscale modeling and seamlessly connect results across different scales and levels of coarse-graining.
All objects in Pizza3 have scripting and static compiling capabilities, and can be converted into LAMMPS instructions through mathematical operators (e.g., +
, |
, *
) or by calling methods like script() or do(). This design philosophy ensures that complex workflows remain accessible, flexible, and modular.
Work in Progress: The toolkit evolves regularly. Come back often for improved capabilities. Newer versions add more classes, dynamic scripting (
pizza.dscript()
,pizza.dforcefield()
), and better integration with large language models via intermediateDSCRIPT SAVE FILE
formats.
- 🛠️ Fully object-oriented architecture
- 🌐 Compatibility with modern Python standards
- 📦 Extensible classes for custom workflows
- 🎯 Focus on multiscale and multiphysics modeling
- 🚀 Advanced 2D/3D design capabilities
- 🌟 Design reusable templates with scripting classes like
pizza.script()
,pizza.dscript()
,pizza.forcefield()
,pizza.dforcefield()
,pizza.group()
. - 🚀 Dynamic design of complex 2D and 3D shapes with
pizza.raster()
,pizza.region()
. - 🧩 Combine simulation components via mathematical operators (
+
,|
,*
,&
). - 📜 Generate input/output scripts and dumps using classes like
pizza.data()
andpizza.dump()
. - 🔬 Perform advanced spatial/temporal analysis using the MATLAB/Octave-based library in the
post/
folder.
Our main goal is to achieve a set of tools that facilitate multiscale modeling in LAMMPS: connecting results across scales with different levels of detail. Below we highlight the current achievements and show how Pizza3 integrates with LAMMPS.
In the latest versions of Pizza3:
- 🚀 Almost all 3D design capabilities of LAMMPS have been implemented, particularly via new classes like
pizza.region()
andpizza.group()
. - 🧪 2D algorithms can be tested via
pizza.raster()
at lower computational cost without losing generality.
The class pizza.script() is central for coding LAMMPS instructions from Python "codelets" and handling
TEMPLATES
. Dynamic scripting concepts viapizza.dscript()
andpizza.dforcefield()
make it easy to transform existing LAMMPS code into reusable templates.
I/O Management:
pizza.data()
and pizza.dump()
objects streamline frame manipulation, output-to-input conversions (like restart files), and advanced spatial/temporal analysis.
Forcefields:
Forcefields can be managed through inheritance, and a user library can be built by overriding existing definitions. It's recommended to implement user rules with pizza.generic() or pizza.dforcefield().
Visualization and Post-Processing:
Some examples of visualization are proposed in 2D and 3D in Python. The post/
folder includes advanced analysis tools integrated with MATLAB/Octave. Users are encouraged to combine LAMMPS-GUI with OVITO for rapid prototyping.
Integration with Large Language Models:
The DSCRIPT SAVE FILE
format has been designed for AI text generators and facilitating the maintenance of a large library of codelets. The use of this new text format or language is not mandatory, but it enables code generation with almost zero Python.
Workshops (workable demos) |
Main classes and subclasses |
Low-level classes |
---|---|---|
class: ==workshop1== | 2D drawing class: pizza.raster |
generic struct class à la Matlab: pizza.private.mstruct() |
class: ==workshop2== | as above (with advanced features: time-steps, region of interest, salting) | as above |
script: geometry | class to read/write input files: pizza.data3.data |
self-evaluable struct with scripting/alias: pizza.private.param() |
script: assembly | manage dump files: pizza.dump3.dump |
|
advanced scripting: pizza.script.script , pizza.script.datascript , pizza.script.scriptobject , pizza.script.scriptobjectgroup , pizza.script.pipescript |
||
forcefields: pizza.forcefield.forcefield , pizza.forcefield.smd , pizza.forcefield.tlsph , pizza.forcefield.ulsph , pizza.forcefield.none , pizza.forcefield.water , pizza.forcefield.solidfood , pizza.forcefield.rigidwall |
Workshop1 illustrates the main steps to design, run, and analyze a LAMMPS project from Python 3.x.
graph TD;
G[geometry] --> R
R[raster] --> F[data.write]
W[workshop1] --> A
F --> A
A[assembly] --> D[dump]
Workshop1 involves 7 main steps and codes, as shown below:
# initizalization of the scheme
bead_kernel_radius = 0.0015
init = initialization(neighbor =[bead_kernel_radius,"bin"])
# scriptobject handles bead interactions
FLUID = scriptdata(
rho = 1000,
c0 = 100.0,
q1 = 1.0,
contact_stiffness = 10000000
)
SOLID = scriptdata(
rho = 2000,
c0 = 200.0,
sigma_yield = '0.1*${E}',
contact_stiffness = 10000000
)
WALL = scriptdata(
rho = 3000,
c0 = 200.0,
contact_stiffness = 10000000,
contact_scale = 1.5
)
b1 = scriptobject(name="bead 1",
group = ["rigid", "solid"],
filename='./raster_2_types.lmp',
forcefield=rigidwall(USER=WALL))
b2 = scriptobject(name="bead 2",
group = ["fluid", "ulsph"],
filename = './raster_2_types.lmp',
forcefield=water(USER=FLUID))
b3 = scriptobject(name="bead 3",
group = ["oscillating", "solid","tlsph"],
filename = './raster_4_types.lmp',
forcefield=solidfood(USER=SOLID))
b4 = scriptobject(name="bead 4",
group = ["solid", "tlsph"],
filename = './raster_4_types.lmp',
forcefield=solidfood(USER=SOLID))
inte = integration()
thermo = thermo_print()
equilsteps = equilibration(it=15)
dmp = smddump(outstep=2000,outputfile=["dump.workshop1"],)
moves = translation(vx = ["0.1*exp(-step/100)"],
vy = ["0"],vz = ["0"]) & \
run() & \
translation() & \
force() & \
run()
collection = b1+b2+b3+b4
fullscript = init + collection.script + inte + thermo + equilsteps + dmp + moves
fullscript.write("./tmp/in.swimmingpool")
Equivalent flowchart:
classDiagram
class initialization{
<<globalsection>>
...
}
class integration{
<<integrationsection>>
g=9.81
...
}
class thermo_print{
<<integrationsection>>
...
}
class equilibration{
<<integrationsection>>
it=50,
re=0.9
}
class smddump{
<<dumpsection>>
outstep=1000,
...
}
class translation{
<<runsection>>
eqvx=0
...
}
class force{
<<runsection>>
eqfx=0
...
}
class run{
<<runsection>>
runs=50000
}
initialization --o collection
collection --o integration
integration --o thermo_print
thermo_print --o equilibration
equilibration --o smddump
translation --o moves
run --o moves
translation --o moves
force --o moves
smddump --o moves
classDiagram
class raster{
...
}
class data{
...
}
class dump{
...
}
class mstruct{
<<private>>
var=value
generator()
scan()
}
class param{
<<private>>
var=value
eval()
formateval()
}
class scriptdata{
var1=value1
}
class forcefield{
...
}
class script{
...
}
class scriptobject{
...
}
class scriptobjectgroup{
...
}
class pipescript{
...
}
data --* raster
mstruct --|> param : extended
mstruct --* dump
param --* paramforcefield
param --|> scriptdata
param --* scriptobject
paramforcefield --* forcefield
forcefield --* scriptobject
scriptobject --|> scriptobjectgroup
scriptdata --* script
script --o pipescript
scriptobject --o pipescript
scriptobjectgroup --o pipescript
data --o dump
# Automatic LAMMPS script (version 0.4)
# olivi@LX-Olivier2021:/home/olivi/billy/python
# ven. 18 mars 2022 19:33:50
# INRAE\SAYFOOD\[email protected]
# INRAE\SAYFOOD\[email protected]
# Scheme initialization
units si
dimension 2
boundary p f p # f = fixed, p = periodic
comm_modify vel yes
comm_style tiled
atom_modify map array
newton off
neighbor 0.0015 bin # search radius for neighbor list (=kernel radius)
neigh_modify every 5 delay 0 check yes
atom_style smd
# ===== [ BEGIN INPUT FILES SECTION ] ===============================================================================
# load files for objects: bead 1, bead 2
read_data ./raster_2_types.lmp
# load files for objects: bead 3, bead 4
read_data ./raster_4_types.lmp add append
# ===== [ END INPUT FILES SECTION ] =================================================================================
# ===== [ BEGIN GROUP SECTION ] =====================================================================================
# Definition of group 1:rigid
# =bead 1
# Similar groups: rigid, solid
group rigid type 1
# Definition of group 2:solid
# =bead 1, bead 3, bead 4
# Similar groups: rigid, solid, oscillating, tlsph
group solid type 1 3 4
# Definition of group 3:fluid
# =bead 2
# Similar groups: fluid, ulsph
group fluid type 2
# Definition of group 4:ulsph
# =bead 2
# Similar groups: fluid, ulsph
group ulsph type 2
# Definition of group 5:oscillating
# =bead 3
# Similar groups: oscillating, solid, tlsph
group oscillating type 3
# Definition of group 6:tlsph
# =bead 3, bead 4
# Similar groups: oscillating, solid, tlsph
group tlsph type 3 4
# ===== [ END GROUP SECTION ] =======================================================================================
# ===== [ BEGIN FORCEFIELD SECTION ] ================================================================================
# [1:bead 1] PAIR STYLE SMD
pair_style hybrid/overlay smd/ulsph *DENSITY_CONTINUITY *VELOCITY_GRADIENT *NO_GRADIENT_CORRECTION &
smd/tlsph smd/hertz 1.5
# [1:bead 1 x 1:bead 1] Diagonal pair coefficient tlsph
pair_coeff 1 1 none
# [2:bead 2 x 2:bead 2] Pair diagonal coefficient ulsph
pair_coeff 2 2 smd/ulsph *COMMON 1000 100.0 1.0 1.0 0 &
*EOS_TAIT 7 &
*END
# [3:bead 3 x 3:bead 3] Diagonal pair coefficient tlsph
pair_coeff 3 3 smd/tlsph *COMMON 2000 400000000.0 0.3 1.0 0.0 10 1.0 &
*STRENGTH_LINEAR_PLASTIC 40000000.0 0 &
*EOS_LINEAR &
*END
# [4:bead 4 x 4:bead 4] Diagonal pair coefficient tlsph
pair_coeff 4 4 smd/tlsph *COMMON 2000 400000000.0 0.3 1.0 0.0 10 1.0 &
*STRENGTH_LINEAR_PLASTIC 40000000.0 0 &
*EOS_LINEAR &
*END
# [1:bead 1 x 2:bead 2] Off-diagonal pair coefficient (generic)
pair_coeff 1 2 smd/hertz 10000000
# [1:bead 1 x 3:bead 3] Off-diagonal pair coefficient (generic)
pair_coeff 1 3 smd/hertz 10000000
# [2:bead 2 x 3:bead 3] Off-diagonal pair coefficient (generic)
pair_coeff 2 3 smd/hertz 10000000
# [1:bead 1 x 4:bead 4] Off-diagonal pair coefficient (generic)
pair_coeff 1 4 smd/hertz 10000000
# [2:bead 2 x 4:bead 4] Off-diagonal pair coefficient (generic)
pair_coeff 2 4 smd/hertz 10000000
# [3:bead 3 x 4:bead 4] Off-diagonal pair coefficient (generic)
pair_coeff 3 4 smd/hertz 10000000
# ===== [ END FORCEFIELD SECTION ] ==================================================================================
# Gravity
fix gfix all gravity 9.81 vector 0 1 0
# Time integration conditions
fix dtfix fluid smd/adjust_dt 0.1 # dynamically adjust time increment every step
fix integration_fix_ulsph ulsph smd/integrate_ulsph adjust_radius 1.01 10 15
fix integration_fix_tlsph tlsph smd/integrate_tlsph
# balancing
fix balance_fix all balance 500 0.9 rcb
# thermodynamic information
thermo 100
thermo_modify lost ignore
thermo_style custom step dt
# Equilibration
dump dump_equilibrate all custom 100 dump.equilibrate id type x y z &
fx fy fz vx vy vz radius
dump_modify dump_equilibrate first yes
variable vmag atom sqrt(vx^2+vy^2+vz^2)
compute maxvel all reduce max v_vmag
variable maxvelre equal 0.6*c_maxvel
run 10000 every 1000 &
"print ${maxvelre}" &
"fix ulsph_equilibration ulsph smd/integrate_ulsph limit_velocity ${maxvelre}" &
"fix tlsph_equilibration tlsph smd/integrate_tlsph limit_velocity ${maxvelre}"
dump_modify dump_equilibrate every 1000000
fix ulsph_equilibration ulsph smd/integrate_ulsph
fix tlsph_equilibration tlsph smd/integrate_tlsph
# Additional computed values
# Compute supp particle data
compute contact_radius all smd/contact/radius
compute vol all smd/vol
compute rho all smd/rho
# dump command
dump dump_id all custom 100 dump.workshop1 &
id type x y z mol mass c_rho c_vol radius c_contact_radius &
vx vy vz fx fy fz
dump_modify dump_id first yes #every/time 0.001
# Translation
variable vx equal 0
variable vy equal 1
variable vz equal 0
fix setvelocities oscillating smd/setvel v_vx v_vy v_vz
# run section
run 50000
# Force
variable fx equal 0/count(all)
variable fy equal 0/count(all)
variable fz equal 0/count(all)
fix setforces all smd/setvel v_fx v_fy v_fz
# run section
run 50000
The videos are generated with ovito directly from dump files.
Equilibration step
Run step
- Inputs from images
- Time-step and salting
- Squeeze and shear examples
With workshop 2, more complex systems should be achievable to the user including the use of "sensing" surfaces and better control of system parameters such as viscosity.
A few classes have been updated/improved to be more usable and funcitonal, and two new classes have been created; file and group. Both these classes give additional methods of creating groups in lammps to suppliment the use of the intereactions class.
The file class can be used to read the files, and it is a new method of creating a group. If you want all the atoms in your file to be part of a group, use the arguments group=True and specify an ID for them with group_name = "your_group_name"
Example:
files = workshop2.file(file_name=["./file1"],group=True,group_name=["sensor"]) & \
workshop2.file(file_name=["./file2"],group=True,group_name=["substance"], append = True)
N.B. for subsequent files, if you want to preserve the files already uploaded use append=True
The group class can create new groups from an input of arguments to specify a sub-domain of the simulation. You have four arguments; x = [x1 x2], x = [x1 x2], x = [x1 x2] and groupID="sensor".
Example:
groups = workshop2.group(groupID = ['bottom'], y=['EDGE',0.2])
By default, all spatial arguments are 'EDGE' which means the edge of the domain box.
- by default, tlsph object are frozen but using static='yourgroup', an group of atoms can be frozen
- equilibrate is simplified, periodically all velocities in the simulation are set to 0 in order to remove the kinetic energy from the system, it= can be used to specify the number of such cycles
- bugs in the translation and force classes resolved
- names can be specified case-by-case and multiple actions can be performed
A block of fluid is squeezed between two surfaces, the lower surface is "sensing" and detects stresses.
pizza.raster()
is a 2D space-filling model for coarse-grained systems.
from pizza.raster import raster, emulsion, coreshell
Use
pizza.region()
to generate 3D complex objects and set their interactions. They can be generated and visualized with LAMMPS-GUI.pizza.group()
applies indifferently to 2D and 3D objects.
The example creates simple objects with different shapes and bead types.
# drawing area
R = raster()
# rectangle shapes
R.rectangle(1,24,2,20,name='rect1')
R.rectangle(60,80,50,81,
name='rect2',beadtype=2,angle=40)
R.rectangle(50,50,10,10,
mode="center",angle=45,beadtype=1)
# circules and ellipses
R.circle(45,20,5,name='C1',beadtype=3)
R.circle(35,10,5,name='C2',beadtype=3) R.circle(15,30,10,
name='p1',beadtype=4,shaperatio=0.2,angle=-30)
R.circle(12,40,8,
name='p2',beadtype=4,shaperatio=0.2,angle=20) R.circle(12,80,22,name='p3',beadtype=4,shaperatio=1.3,angle=20)
# other polygons
R.triangle(85,20,10,name='T1',beadtype=5,angle=20)
R.diamond(85,35,5,name='D1',beadtype=5,angle=20)
R.pentagon(50,35,5,name='P1',beadtype=5,angle=90)
R.hexagon(47,85,12,name='H1',beadtype=5,angle=90)
# labels and first plot
R.label("rect003")
R.plot()
R.list()
R.show()
# refresh with all labels except one
R.clear()
R.plot()
R.show(extra="label")
R.label("rect003")
R.unlabel('rect1')
# generate and write the corresponding pizza.data() object
X=R.data()
X.write("tmp/example.lmp") # scaling is available
Geometry objects can be copied between pizza.raster()
objects and copied along arbitrary paths or coordinates.
# source
draft = raster()
draft.rectangle(1,24,2,20,name='rect1'),
draft.circle(10,10,2,name="X",beadtype=4)
# copy an entire raster to the collection "C1"
A = raster()
A.collection(draft,name="C1",beadtype=1,translate=[10,30])
# copy along a path "PX"
B = raster()
B.copyalongpath(draft.X,name="PX",beadtype=2,
path=arc,
xmin=10,
ymin=10,
xmax=90,
ymax=50,
n=12)
B.plot()
B.show(extra="label")
# generate the suspension/emulsion object with multiple random insertions by indicating the radius of the object
e = emulsion(xmin=10, ymin=10, xmax=390, ymax=390)
e.insertion([60,50,40,30,20,15,15,10,8,20,12,8,6,4,11,13],beadtype=1)
e.insertion([30,10,20,2,4,5,5,10,12,20,25,12,14,16,17],beadtype=2)
e.insertion([40,2,8,6,6,5,5,2,3,4,4,4,4,4,10,16,12,14,13],beadtype=3)
# generate the corresponing raster
C = raster(width=400,height=400)
C.scatter(e,name="emulsion")
C.plot()
C.show()
# generate the core-shell model
cs = coreshell(xmin=10, ymin=10, xmax=390, ymax=390)
cs.insertion([60,50,40,30,20,15,15,10,8,20,12,8,11,13],beadtype=(1,2),thickness = 4)
# generate the corresponding raster
D = raster(width=400,height=400)
D.scatter(cs,name="core-shell")
D.plot()
D.show()
Forcefield customization with inheritance. Example outputs shown previously.
from pizza.forcefield import *
# basic examples from user library
w = water(beadtype=1, userid="fluid")
w.parameters.Cp = 20
print("\n"*2,w)
f = solidfood(beadtype=2, userid="elastic")
print("\n"*2,f)
r = rigidwall(beadtype=3, userid="wall")
print("\n"*2,r)
Example of output for r
============================ [ wall | version=0.1 ] ============================
Bead of type 3 = [LAMMPS:SMD:none:walls]
forcefield (FF object) with 4 parameters
............................... [ description ] ................................
# LAMMPS:SMD - solid, liquid, rigid forcefields (continuum mechanics)
# no interactions
# rigid walls
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [ methods ] ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
replace FFi,FFj by your variable names <<<
To assign a type, use: FFi.beadtype = integer value
Use the methods FFi.pair_style() and FFi.pair_coeff(FFj)
Note for pairs: the caller object is i (FFi), the argument is j (FFj or j)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [ template ] ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# [3:wall] PAIR STYLE SMD
pair_style hybrid/overlay smd/ulsph *DENSITY_CONTINUITY *VELOCITY_GRADIENT *NO_GRADIENT_CORRECTION &
smd/tlsph smd/hertz 1.5
# [3:wall x 3:wall] Diagonal pair coefficient tlsph
pair_coeff 3 3 none
# [3:wall x 2:none] Off-diagonal pair coefficient (generic)
pair_coeff 2 3 smd/hertz 750000.0
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Key scripting classes enabling dynamic variable space, template parsing, and operator overloading.
from pizza.script import *
classes | description | overloaded operators and attributes |
---|---|---|
scriptdata() |
variable space: |
|
scriptobject() |
"bead" object: |
name="bead name", group=["group1","group2"], filename="input_filename.lmp", forcefield=forcefield_class |
scriptobjectgroup() |
collection of beads: |
select(),group(), interactions(), forcefield(), script $P=C_1 |
script() |
script object: |
do(),write() Immediate execution: Differed execution: |
pipescript() |
pipescript object: |
do(),script,clear,rename() $P=P_0 |
Easily handle input/output files, including restarts:
# input data objects and methods
from pizza.data3 import data
# dump objects and methods
from pizza.dump3 import dump
# basic example, see documentation for details
datafile = "../data/play_data/data.play.lmp"
X = data(datafile)
Y = dump("../data/play_data/dump.play.restartme")
t = Y.time()
step = 2000
R = data(Y,step)
R.write("../tmp/data.myfirstrestart.lmp")
🖋️ Contributors:
- INRAE\Olivier Vitrac (main contact: [email protected])
- INRAE\William Jenkinson ([email protected])
- INRAE\Pei Sun ([email protected])
- INRAE\Han Chen ([email protected])
- community contributors
$ 2024-12-21 $
For further details, consult the online documentation and the release page for new capabilities.