Skip to content

Commit 57f9265

Browse files
sirmarcelclaude
andcommitted
Add i-PI driven dynamics example with BEC support
Wire up the Calculator → i-PI BEC data path: extract "apt" from LoremBEC model output, reshape to (3*natoms, 3), and expose as "BEC" for i-PI's driven_dynamics module. Implement LOREMCalculator as a proper Calculator subclass and clean up the i-PI driver (remove buggy compute_structure override that operated on JSON strings). Add md-ipi example with input.xml (eda-nve, bec mode="driver"), run.sh, start.xyz, and README. Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
1 parent 869b070 commit 57f9265

File tree

8 files changed

+204
-13
lines changed

8 files changed

+204
-13
lines changed

examples/md-ipi/README.md

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
# Driven dynamics with LOREM via i-PI
2+
3+
This example shows how to run an applied electric field simulation using a LOREM BEC model and i-PI's driven dynamics (Electric Dipole Approximation).
4+
5+
Unlike the standard i-PI efield example which reads fixed Born Effective Charges from a file, LOREM computes BECs on-the-fly from its learned model (`<bec mode="driver"/>`).
6+
7+
## Prerequisites
8+
9+
- Trained `LoremBEC` checkpoint (with `model/model.yaml`, `model/baseline.yaml`, `model/model.msgpack`)
10+
- i-PI installed (`pip install ipi`)
11+
- LOREM i-PI driver installed: `lorem-install-ipi-driver`
12+
13+
## Files
14+
15+
- `input.xml` — i-PI configuration for driven dynamics (EDA-NVE). Adapt the E-field parameters (`amp`, `freq`, `peak`, `sigma`), cell size, and `start.xyz` to your system.
16+
- `start.xyz` — Starting structure in atomic units. Replace with your system's geometry.
17+
- `run.sh` — Launch script. Set `MODEL_PATH` to your trained checkpoint folder.
18+
19+
## Running
20+
21+
```bash
22+
# Edit run.sh to set MODEL_PATH, then:
23+
bash run.sh
24+
```
25+
26+
i-PI outputs will be written to `i-pi.*` files. Key outputs:
27+
- `i-pi.properties.out` — energy, kinetic energy, E-field strength over time
28+
- `i-pi.positions_*` — nuclear trajectories
29+
- `i-pi.bec{x,y,z}_*` — BEC tensor components over time
30+
31+
## Adapting to your system
32+
33+
1. Replace `start.xyz` with your starting geometry (atomic units).
34+
2. Set the cell size in `input.xml` (large for isolated systems, physical for periodic).
35+
3. Set `pbc='True'` in `<ffsocket>` if your system is periodic.
36+
4. Adjust E-field parameters to match your target excitation.
37+
5. Point `MODEL_PATH` in `run.sh` to your trained LoremBEC checkpoint.

examples/md-ipi/input.xml

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,55 @@
1+
<simulation verbosity='high'>
2+
3+
<total_steps>500</total_steps>
4+
5+
<ffsocket mode='unix' name='driver' pbc='False'>
6+
<address>lorem</address>
7+
</ffsocket>
8+
9+
<output prefix='i-pi'>
10+
<trajectory filename='forces'>forces</trajectory> <!-- total forces (a.u.): nuclear + q_e Z @ E(t) -->
11+
<trajectory filename='Eforces'>Eforces</trajectory> <!-- driven forces (a.u.): q_e Z @ E(t) -->
12+
<trajectory filename='momenta'>momenta</trajectory> <!-- nuclear momenta (a.u.) -->
13+
<trajectory filename='velocities'>velocities</trajectory> <!-- nuclear velocities (a.u.) -->
14+
<trajectory filename='positions'>positions</trajectory> <!-- nuclear positions (a.u.) -->
15+
<trajectory filename='becx'>becx</trajectory>
16+
<trajectory filename='becy'>becy</trajectory>
17+
<trajectory filename='becz'>becz</trajectory>
18+
<trajectory filename='dip' extra_type='dipole'> extras </trajectory>
19+
<properties filename='properties.out' verbosity='low' > [ step, time, conserved, kinetic_md, potential, Efield, Eenvelope ] </properties>
20+
</output>
21+
22+
<prng>
23+
<seed>10545</seed>
24+
</prng>
25+
26+
<system>
27+
<forces>
28+
<force forcefield='driver' />
29+
</forces>
30+
31+
<initialize nbeads='1'>
32+
<file mode='xyz' units='atomic_unit' > start.xyz </file>
33+
<velocities mode='thermal' units='kelvin'> 0 </velocities> <!-- no initial velocities: excited by E-field only -->
34+
<cell mode="manual" units="atomic_unit">[ 50,0,0, 0,50,0, 0,0,50 ]</cell>
35+
</initialize>
36+
37+
<ensemble>
38+
<temperature units='kelvin'>0</temperature>
39+
</ensemble>
40+
41+
<motion mode='driven_dynamics'>
42+
<driven_dynamics mode='eda-nve'> <!-- Electric Dipole Approximation (EDA) on top of NVE -->
43+
<timestep units='femtosecond'> 1 </timestep>
44+
<efield>
45+
<amp> [ 0 , 5e-2 , 0 ] </amp> <!-- amplitude of the electric field pulse in a.u. -->
46+
<freq units="thz"> 115 </freq> <!-- frequency of the electric field pulse in THz -->
47+
<peak units="picosecond"> 0.2 </peak> <!-- peak of the electric field pulse -->
48+
<sigma units="picosecond"> 0.05 </sigma> <!-- width (std of gaussian envelope) -->
49+
</efield>
50+
<bec mode="driver"/> <!-- BECs computed by LOREM, not read from file -->
51+
</driven_dynamics>
52+
</motion>
53+
54+
</system>
55+
</simulation>

examples/md-ipi/run.sh

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
#!/bin/bash
2+
# Driven dynamics with LOREM BEC model via i-PI.
3+
#
4+
# Prerequisites:
5+
# - i-PI installed (`pip install ipi`)
6+
# - LOREM driver installed (`lorem-install-ipi-driver`)
7+
# - Trained LoremBEC checkpoint
8+
#
9+
# Adapt MODEL_PATH and start.xyz / input.xml to your system.
10+
11+
MODEL_PATH="/path/to/checkpoint"
12+
13+
# Install LOREM driver into i-PI (idempotent)
14+
lorem-install-ipi-driver
15+
16+
# Start i-PI server
17+
i-pi input.xml > i-pi.out &
18+
echo "i-PI started"
19+
20+
# Wait for socket to be ready
21+
sleep 5
22+
23+
# Start LOREM driver
24+
i-pi-driver -a lorem -u -m lorem \
25+
-o model_path=${MODEL_PATH},template=start.xyz \
26+
> driver.out &
27+
echo "LOREM driver started"
28+
29+
wait
30+
echo "Simulation complete"

examples/md-ipi/start.xyz

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
3
2+
positions in atomic_units (ground state water, PES from pswater)
3+
O 3.77946e+00 3.41870e+00 3.77945e+00
4+
H 3.77945e+00 4.52675e+00 5.21072e+00
5+
H 3.77945e+00 4.52675e+00 2.34819e+00

src/lorem/calculator.py

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -126,6 +126,15 @@ def calculate(
126126
elif key == "stress":
127127
raise KeyError
128128

129+
# BEC passthrough: when model outputs "apt" (e.g. LoremBEC),
130+
# expose as "BEC" in (3*natoms, 3) layout for i-PI compatibility
131+
if "apt" in results:
132+
apt = np.array(
133+
results["apt"][self.batch.sr.atom_mask].reshape(-1, 3, 3),
134+
dtype=np.float32,
135+
)
136+
actual_results["BEC"] = apt.reshape(-1, 3)
137+
129138
if self.add_offset:
130139
energy_offset = np.sum(
131140
[self.species_weights[Z] for Z in atoms.get_atomic_numbers()]
@@ -136,7 +145,7 @@ def calculate(
136145
return actual_results
137146

138147
def get_property(self, name, atoms=None, allow_calculation=True):
139-
if name not in self.implemented_properties:
148+
if name not in self.implemented_properties and name != "BEC":
140149
raise PropertyNotImplementedError(f"{name} property not implemented")
141150

142151
self.update(atoms)

src/lorem/dynamics/ase.py

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,7 @@
1-
class LOREMCalculator:
2-
@classmethod
3-
def from_checkpoint(cls, model_path):
4-
pass
1+
from lorem.calculator import Calculator
2+
3+
4+
class LOREMCalculator(Calculator):
5+
"""LOREM ASE calculator with BEC support for i-PI integration."""
6+
7+
pass

src/lorem/dynamics/ipi.py

Lines changed: 0 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -21,11 +21,3 @@ def __init__(self, template, model_path, device="cpu", *args, **kwargs):
2121
def check_parameters(self):
2222
super().check_parameters()
2323
self.ase_calculator = LOREMCalculator.from_checkpoint(self.model_path)
24-
25-
def compute_structure(self, cell, pos):
26-
pot_ipi, force_ipi, vir_ipi, extras = super().compute_structure(cell, pos)
27-
28-
if "BEC" not in extras and "apt" in extras:
29-
extras["BEC"] = extras["apt"]
30-
31-
return pot_ipi, force_ipi, vir_ipi, extras

tests/test_calculator.py

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
import numpy as np
2+
3+
from ase.build import bulk
4+
5+
from lorem.calculator import Calculator
6+
from lorem.models.bec import LoremBEC
7+
from lorem.models.mlip import Lorem
8+
9+
10+
def test_calculator_energy_forces():
11+
model = Lorem(cutoff=5.0, num_features=8, num_spherical_features=2, num_radial=4)
12+
calc = Calculator.from_model(model)
13+
14+
atoms = bulk("Ar") * [2, 2, 2]
15+
calc.calculate(atoms)
16+
17+
assert "energy" in calc.results
18+
assert "forces" in calc.results
19+
assert calc.results["forces"].shape == (len(atoms), 3)
20+
assert "BEC" not in calc.results # Lorem (not LoremBEC) has no BEC
21+
22+
23+
def test_calculator_bec():
24+
model = LoremBEC(cutoff=5.0, num_features=8, num_spherical_features=2, num_radial=4)
25+
calc = Calculator.from_model(model)
26+
27+
atoms = bulk("Ar") * [2, 2, 2]
28+
calc.calculate(atoms)
29+
30+
assert "energy" in calc.results
31+
assert "forces" in calc.results
32+
assert "BEC" in calc.results
33+
34+
natoms = len(atoms)
35+
assert calc.results["BEC"].shape == (3 * natoms, 3)
36+
assert calc.results["BEC"].dtype == np.float32
37+
38+
39+
def test_calculator_bec_get_property():
40+
model = LoremBEC(cutoff=5.0, num_features=8, num_spherical_features=2, num_radial=4)
41+
calc = Calculator.from_model(model)
42+
43+
atoms = bulk("Ar") * [2, 2, 2]
44+
calc.calculate(atoms)
45+
46+
bec = calc.get_property("BEC", atoms)
47+
assert bec.shape == (3 * len(atoms), 3)
48+
49+
50+
def test_lorem_calculator_import():
51+
from lorem import LOREMCalculator
52+
53+
model = LoremBEC(cutoff=5.0, num_features=8, num_spherical_features=2, num_radial=4)
54+
calc = LOREMCalculator.from_model(model)
55+
56+
atoms = bulk("Ar") * [2, 2, 2]
57+
calc.calculate(atoms)
58+
59+
assert "BEC" in calc.results
60+
assert calc.results["BEC"].shape == (3 * len(atoms), 3)

0 commit comments

Comments
 (0)