3
$\begingroup$

I want to read a lammps data file in ASE, for instance the TIP3P or the SPC model. The lammps examples provide the input script and the data file, which work perfectly.

ASE is not able to read the .mol file: it cannot parse neither Shake *** nor Shake *** section of the lammps-data files. Thus I removed them, but even then UnboundLocalError: local variable 'mol_id_in' referenced before assignment. Finally I wrote my own file, which can be read - see it at the bottom.

Once the file is read, I set and apply a LammpsLib calculator with exactly the same inputs as in the Lammps tutorial, and I receive a Exception: ERROR: Boundary command after simulation box is defined (src/input.cpp:1414).

I think that the parsing messes up the order of information.

Question: How can I see what is the input file that is fed to lammps?

# Read file

import ase.io
from ase.calculators.lammpslib import LAMMPSlib

water = ase.io.read('water.mol', format='lammps-data')

lammps_orig_header = [
    'units real',
    'atom_style full',
    'region box block -5 5 -5 5 -5 5',
    'create_box 2 box  bond/types 1 angle/types 1 extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2'
    ]

lammps_orig_commands = [
    'mass 1 15.9994',
    'mass 2 1.008',
    
    'pair_style lj/cut/coul/cut 10.0',
    'pair_coeff 1 1 0.1553 3.166',
    'pair_coeff 1 2 0.0    1.0',
    'pair_coeff 2 2 0.0    1.0',
    
    'bond_style zero',
    'bond_coeff 1 1.0',
    
    'angle_style zero',
    'angle_coeff 1 109.47',
    
    'molecule water spce.mol',
    'create_atoms 0 random 33 34564 NULL mol water 25367 overlap 1.33',
    
    'timestep 1.0',
    'fix rigid     all shake 0.0001 10 10000 b 1 a 1',
    'minimize 0.0 0.0 1000 10000',
    'velocity all create 300.0 5463576',
    'fix integrate all nvt temp 300.0 300.0 100.0',
    
    'thermo_style custom step temp press etotal density pe ke',
    'thermo 1000',
    'run 20000 upto',
    'write_data spce.data nocoeff'
    ]

lammps = LAMMPSlib( lammps_header=lammps_orig_header,lmpcmds=lammps_orig_commands,log_file='test.log')

# Set calculator and calculate energy
water.set_calculator(lammps)
water.get_potential_energy()

Error:

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Cell In[1], line 49
     47 # Set calculator and calculate energy
     48 water.set_calculator(lammps)
---> 49 water.get_potential_energy()

File ~/WORK/external_packages/ase/ase/atoms.py:737, in Atoms.get_potential_energy(self, force_consistent, apply_constraint)
    734     energy = self._calc.get_potential_energy(
    735         self, force_consistent=force_consistent)
    736 else:
--> 737     energy = self._calc.get_potential_energy(self)
    738 if apply_constraint:
    739     for constraint in self.constraints:

File ~/WORK/external_packages/ase/ase/calculators/abc.py:24, in GetPropertiesMixin.get_potential_energy(self, atoms, force_consistent)
     22 else:
     23     name = 'energy'
---> 24 return self.get_property(name, atoms)

File ~/WORK/external_packages/ase/ase/calculators/calculator.py:537, in BaseCalculator.get_property(self, name, atoms, allow_calculation)
    534     if self.use_cache:
    535         self.atoms = atoms.copy()
--> 537     self.calculate(atoms, [name], system_changes)
    539 if name not in self.results:
    540     # For some reason the calculator was not able to do what we want,
    541     # and that is OK.
    542     raise PropertyNotImplementedError(
    543         '{} not present in this ' 'calculation'.format(name)
    544     )

File ~/WORK/external_packages/ase/ase/calculators/lammpslib.py:360, in LAMMPSlib.calculate(self, atoms, properties, system_changes)
    359 def calculate(self, atoms, properties, system_changes):
--> 360     self.propagate(atoms, properties, system_changes, 0)

File ~/WORK/external_packages/ase/ase/calculators/lammpslib.py:381, in LAMMPSlib.propagate(self, atoms, properties, system_changes, n_steps, dt, dt_not_real_time, velocity_field)
    379     self.start_lammps()
    380 if not self.initialized:
--> 381     self.initialise_lammps(atoms)
    382 else:  # still need to reset cell
    383     # NOTE: The whole point of ``post_changebox_cmds`` is that they're
    384     # executed after any call to LAMMPS' change_box command.  Here, we
   (...)
    391     # apply any shrink-wrapping *after* it's updated the cell
    392     # dimensions
    393     if 'pbc' in system_changes:

File ~/WORK/external_packages/ase/ase/calculators/lammpslib.py:668, in LAMMPSlib.initialise_lammps(self, atoms)
    666             break
    667     else:
--> 668         self.lmp.command('boundary ' + self.lammpsbc(atoms))
    670 # Initialize cell
    671 self.set_cell(atoms, change=not self.parameters.create_box)

File ~/miniconda3/envs/pgel/lib/python3.10/site-packages/lammps/core.py:617, in lammps.command(self, cmd)
    614 if cmd: cmd = cmd.encode()
    615 else: return
--> 617 with ExceptionCheck(self):
    618   self.lib.lammps_command(self.lmp,cmd)

File ~/miniconda3/envs/pgel/lib/python3.10/site-packages/lammps/core.py:49, in ExceptionCheck.__exit__(self, exc_type, exc_value, traceback)
     47 def __exit__(self, exc_type, exc_value, traceback):
     48   if self.lmp.has_exceptions and self.lmp.lib.lammps_has_error(self.lmp.lmp):
---> 49     raise self.lmp._lammps_exception

Exception: ERROR: Boundary command after simulation box is defined (src/input.cpp:1414)
# Water molecule. TIP3P geometry

3 atoms
2 bonds
1 angles

Coords

1    0.00000  -0.06556   0.00000
2    0.75695   0.52032   0.00000
3   -0.75695   0.52032   0.00000

Types

1        1   # O
2        2   # H
3        2   # H

Charges

1       -0.834
2        0.417
3        0.417

Bonds

1   1      1      2
2   1      1      3

Angles

1   1      2      1      3
$\endgroup$
1
  • $\begingroup$ +1 for a well described question. It feels great to see a question with relevant details such as error codes, sample files. $\endgroup$ Commented Nov 12, 2023 at 7:43

0

You must log in to answer this question.