2
$\begingroup$

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

lammps_header=[
    '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'
]

lammps = LAMMPSlib( lammps_header=lammps_header,lmpcmds=lmpcmds )
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
> ...

```
$\endgroup$
3
  • $\begingroup$ I don't see bond and angle terms -- without those, there is nothing stopping the oxygens and hydrogens from crashing straight into each other. $\endgroup$ Commented Nov 10, 2023 at 10:13
  • $\begingroup$ 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. $\endgroup$ Commented Nov 10, 2023 at 10:32
  • $\begingroup$ 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. $\endgroup$ Commented Nov 10, 2023 at 10:41

0

You must log in to answer this question.