# Water simulation in ASE/LAMMPS

I want to simulate some water in MD with ASE and a LammpsLib calculator.

I set up a minimal Lammps calculator and I try to calculate the energy per atom of an isolated molecule, and a 3x3x3 supercell. I get different energies. Question 1: How is it possible?

Moreover, when I try to optimize either system, the cell explodes reaching crazy high forces and energy oscillating of several orders of magnitude. Question 2: Is my Lammps calculator complete?

Perhaps some fix?

Thank you

from ase.calculators.lammpslib import LAMMPSlib

'units metal',
'boundary p p p',
'atom_style full'
]

lmpcmds = [
'mass 1 15.9994',
'mass 2 1.008',
'set type 1 charge -0.834',
'set type 2 charge 0.417',
'pair_style lj/cut/coul/long 8.0',
'pair_coeff 1 1 0.1521 3.1507',
'pair_coeff 2 2 0.0    1.0',
'kspace_style pppm 1.0e-4',
'minimize 0.0 0.0 1000 10000'
]


from ase.build import molecule

cell_size = 3.
water = molecule('H2O', vacuum=3.0)
water.set_cell([cell_size, cell_size, cell_size])
water.set_pbc([True, True, True])

# Set calculator and calculate energy
water.set_calculator(lammps)
water.get_potential_energy()/len(water.numbers) # = -2.901683325464138

import ase.build

# Define the supercell
sxcell_size = 3
sxcell_matrix = [[sxcell_size, 0, 0],
[0, sxcell_size, 0],
[0, 0, sxcell_size]]
sxcell = ase.build.make_supercell(water, sxcell_matrix)

# Set calculator and calculate energy
sxcell.set_calculator(lammps)
sxcell.get_potential_energy()/len(sxcell.numbers) # = -0.0899526348900734

import ase.optimize

dyn = ase.optimize.BFGS(water)
dyn.run(fmax=1)

> Step     Time          Energy          fmax
> BFGS:    0 10:40:43       -8.705050        5.815941
> BFGS:    1 10:40:43       -9.953564        7.020505
> ...
> BFGS:  176 10:40:45   -17347.959238 60056340.171811
> BFGS:  177 10:40:45    -4428.006468  3905232.089535
> ...

$$$$

• I don't see bond and angle terms -- without those, there is nothing stopping the oxygens and hydrogens from crashing straight into each other. Commented Nov 10, 2023 at 10:13
• Thank you. how could I know how much information do I need to provide lammps in advance? The bond and the angle could in principle be known by the molecule provided. Commented Nov 10, 2023 at 10:32
• I tried to insert bonds and angles, which lead to ERROR: Numeric index is out of bounds. I therefore added region box and a create_box commands, which in turn lead to ERROR: Boundary command after simulation box is defined`. It looks like ASE and Lammps are conflicting on who is defining the geometry. Commented Nov 10, 2023 at 10:41