Pyiron automatically assigns bonds in the input files for lammps for O and H whenever they exist in the structure

[sudarsan-surendralal]

Issue from StackOverflow: https://stackoverflow.com/questions/71846448/pyiron-automatically-assigns-bonds-in-the-input-files-for-lammps-for-o-and-h-whe

Whenever I purposefully put H atoms inside a structure (with Fe and O in the structure as host atoms) as interstitials, I expect that it will be defined as an isolated atom with no bonds between it and the surrounding host atoms. However, depending on the location of the H interstitial, pyiron sometimes defines a bond between the additional H and the original O for lammps calculation.

This can be useful for the automatic detection of bonds. However, how can one control this feature whenever not needed?

I also believe there is a bug with the .define_bonds() function. The inputs for the functions are:

species="O" element_list=["H"] cutoff_list=[2.0] max_bond_list=1 bond_type_list=1

While I believe the max_bond_list and bond_type_list are supposed to be lists, it only works error-free if both are defined as integers instead because of the way they are defined. However, defining them as integers then messes up running the jobs because they should be iterables.

[sudarsan-surendralal]

@ahmedabdelkawy Can you post the exact code snippet that geives you this error please

[ahmedabdelkawy]

Sure!

############################### Code reproducibility setup ###################################

pr = Project(path='bond_defintion')

#structure and interstatial definition 
structure_fe2o3_pug = ase_to_pyiron(ase.io.read('Fe2O3_Primitive_cell.txt',format='vasp'))
H_position = np.array([ 0.2519    , -0.14543453,  0.91813333])      #originaly obtained from a mesh that is generated according to cell dimensions
H_interstatial = pr.create.structure.atoms("H ", positions=[H_position], cell=structure_fe2o3_pug.get_cell())
structure_int = structure_fe2o3_pug.copy()
structure_int.append(H_interstatial)
structure_int.set_repeat([3,3,3])       # duplicated the cell for lammps minimum image convention

#Potential definition
#I use custom openkim potential that does not exist in pyiron
custom_potential = pd.DataFrame({
  'Name': ['Fe_O_Reax'],
  'Filename': [[]],
  'Model': ['Custom'],
  'Species': [['Fe','O','H']],
  'Config': [[ 'atom_style full\n',  # I use 'full' here as atom_style 'charge' gives the same result
              '## create groups ###\n',
              'group Fe type 1\n',
              'group O type 2\n',
              'group H type 3\n',
              'set group Fe charge 2.000\n',
              'set group O charge -2.000\n',
              'set group H charge 1.000\n',
              'kim interactions Fe O H\n'              
             ]]
})

#Lammps job
job_lammps_fe2o3 = pr.create_job(job_type=pr.job_type.Lammps, job_name='lammps_fe2o3_',delete_existing_job=True)
job_lammps_fe2o3.input.control.read_input(file_name='config.txt')
job_lammps_fe2o3.structure = structure_int.copy() 
job_lammps_fe2o3.potential = custom_potential
job_lammps_fe2o3.calc_static()
job_lammps_fe2o3.server.queue = "cmti"
job_lammps_fe2o3.run(delete_existing_job=True, run_mode='queue')

############################### Issues ###################################

job_lammps_fe2o3.bond_dict;               #Bonds are defined automatically 

job_lammps_fe2o3.decompress()
job_lammps_fe2o3['structure.inp'];       #Automatically defined bonds can also be seen in the input file

#Here if you define max_bond_list, and bond_type_list as lists: [1] it gives an error about their length
job_lammps_fe2o3.define_bonds(species="O", element_list=["H"], cutoff_list=[2.0], max_bond_list=[1], bond_type_list=[1])

#If you define them as integars it works 
#job_lammps_fe2o3.define_bonds(species="O", element_list=["H"], cutoff_list=[0.0], max_bond_list=1, bond_type_list=1)

#Now if you try to run it again based on the second definition of job_lammps_fe2o3.define_bonds it will crash because pyiron expects bond_type_list to be iterable 

job_lammps_fe2o3 = pr.create_job(job_type=pr.job_type.Lammps, job_name='lammps_fe2o3_',delete_existing_job=True)
job_lammps_fe2o3.input.control.read_input(file_name='config.txt')
job_lammps_fe2o3.structure = structure_int.copy() 

job_lammps_fe2o3.define_bonds(species="O", element_list=["H"], cutoff_list=[0.0], max_bond_list=1, bond_type_list=1)

job_lammps_fe2o3.potential = custom_potential
job_lammps_fe2o3.calc_static()
job_lammps_fe2o3.server.queue = "cmti"
job_lammps_fe2o3.run(delete_existing_job=True, run_mode='queue')

############################### Files ###################################

Fe2O3_Primitive_cell.txt config.txt

[ahmedabdelkawy]

If you commented out the set_repeat so that you only have one cell in the structure and tried structure_int.get_all_distances()[-1] where it prints out all distances between all the atoms in the structure and the additional hydrogen atom (last atom) you get array([ 3.98568602, 10.86495244, 1.11334241, 7.96598722, 7.39611106, 5.21884493, 7.32919853, 5.04723849, 7.27603826, 5.14392108, 0.])

All distances are more than the default cutoff of 2A except for the third one of 1.11334241 A. However, the two defined bonds in lammps are of length 5.21884493A and 5.04723849A

[sudarsan-surendralal]

'config.txt' seems to be another POSCAR structure. Perhaps you meant another file?

[ahmedabdelkawy]

Sorry, that is a mistake, here is the correct file!

config.txt

[sudarsan-surendralal]

I can confirm that this is a bug which I think I've fixed (see https://github.com/pyiron/pyiron_atomistics/pull/592). Can you run the following snippet (with atom style charge) using that branch to confirm?

pr = Project(path='bond_defintion')

#structure and interstatial definition 
structure_fe2o3_pug = ase_to_pyiron(ase.io.read('Fe2O3_Primitive_cell.txt',format='vasp'))
H_position = np.array([ 0.2519    , -0.14543453,  0.91813333])      #originaly obtained from a mesh that is generated according to cell dimensions
H_interstatial = pr.create.structure.atoms("H ", positions=[H_position], cell=structure_fe2o3_pug.get_cell())
structure_int = structure_fe2o3_pug.copy()
structure_int.append(H_interstatial)
structure_int.set_repeat([3,3,3])       # duplicated the cell for lammps minimum image convention

#Potential definition
#I use custom openkim potential that does not exist in pyiron
custom_potential = pd.DataFrame({
  'Name': ['Fe_O_Reax'],
  'Filename': [[]],
  'Model': ['Custom'],
  'Species': [['Fe','O','H']],
  'Config': [[ 'atom_style charge\n',  # I use 'full' here as atom_style 'charge' gives the same result
              '## create groups ###\n',
              'group Fe type 1\n',
              'group O type 2\n',
              'group H type 3\n',
              'set group Fe charge 2.000\n',
              'set group O charge -2.000\n',
              'set group H charge 1.000\n',
              'kim interactions Fe O H\n'           
             ]]
})

# Setting the individual charges for the atoms
charges = np.zeros(len(structure_int))
charges[structure_int.select_index("O")] = -2
charges[structure_int.select_index("Fe")] = 2
charges[structure_int.select_index("H")] = 1

structure_int.set_initial_charges(charges)
structure_int.get_initial_charges()

# Running the job
job_lammps_fe2o3 = pr.create_job(job_type=pr.job_type.Lammps, job_name='lammps_fe2o3_',delete_existing_job=True)
job_lammps_fe2o3.input.control.read_input(file_name='config.txt')
job_lammps_fe2o3.structure = structure_int.copy() 
job_lammps_fe2o3.potential = custom_potential
job_lammps_fe2o3.calc_static()
job_lammps_fe2o3.run()

[ahmedabdelkawy]

When I run using atom style charge, I get this error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [28], in <cell line: 7>()
      5 job_lammps_fe2o3.potential = custom_potential
      6 job_lammps_fe2o3.calc_static()
----> 7 job_lammps_fe2o3.run()

File /u/system/SLES12/soft/pyiron/dev/anaconda3/lib/python3.8/site-packages/pyiron_base/generic/util.py:213, in Deprecator.__deprecate_argument.<locals>.decorated(*args, **kwargs)
    203     if kw in self.arguments:
    204         warnings.warn(
    205             message_format.format(
    206                 "{}.{}({}={})".format(
   (...)
    211             stacklevel=2,
    212         )
--> 213 return function(*args, **kwargs)

File /u/system/SLES12/soft/pyiron/dev/anaconda3/lib/python3.8/site-packages/pyiron_base/job/generic.py:704, in GenericJob.run(self, delete_existing_job, repair, debug, run_mode, run_again)
    702     self._run_if_repair()
    703 elif status == "initialized":
--> 704     self._run_if_new(debug=debug)
    705 elif status == "created":
    706     self._run_if_created()

File /u/system/SLES12/soft/pyiron/dev/anaconda3/lib/python3.8/site-packages/pyiron_base/job/generic.py:1417, in GenericJob._run_if_new(self, debug)
   1415     print("job exists already and therefore was not created!")
   1416 else:
-> 1417     self.save()
   1418     self.run()

File /u/system/SLES12/soft/pyiron/dev/anaconda3/lib/python3.8/site-packages/pyiron_base/job/generic.py:1268, in GenericJob.save(self)
   1266 if self._check_if_input_should_be_written():
   1267     self.project_hdf5.create_working_directory()
-> 1268     self.write_input()
   1269     _copy_restart_files(job=self)
   1270 self.status.created = True

File /u/system/SLES12/soft/pyiron/dev/anaconda3/lib/python3.8/site-packages/pyiron_atomistics/lammps/base.py:366, in LammpsBase.write_input(self)
    364 if self.structure is None:
    365     raise ValueError("Input structure not set. Use method set_structure()")
--> 366 lmp_structure = self._get_lammps_structure(
    367     structure=self.structure, cutoff_radius=self.cutoff_radius
    368 )
    369 lmp_structure.write_file(file_name="structure.inp", cwd=self.working_directory)
    370 version_int_lst = self._get_executable_version_number()

File /u/system/SLES12/soft/pyiron/dev/anaconda3/lib/python3.8/site-packages/pyiron_atomistics/lammps/base.py:1248, in LammpsBase._get_lammps_structure(self, structure, cutoff_radius)
   1245     return lammps_structure
   1247 if structure is not None:
-> 1248     lmp_structure.structure = structure_to_lammps(structure)
   1249 else:
   1250     lmp_structure.structure = structure_to_lammps(self.structure)

File /u/system/SLES12/soft/pyiron/dev/anaconda3/lib/python3.8/site-packages/pyiron_atomistics/lammps/structure.py:236, in LammpsStructure.structure(self, structure)
    234     input_str = self.structure_bond()
    235 elif self.atom_type == "charge":
--> 236     input_str = self.structure_charge()
    237 else:  # self.atom_type == 'atomic'
    238     input_str = self.structure_atomic()

File /u/system/SLES12/soft/pyiron/dev/anaconda3/lib/python3.8/site-packages/pyiron_atomistics/lammps/structure.py:591, in LammpsStructure.structure_charge(self)
    587 atoms = "Atoms\n\n"
    589 coords = self.rotate_positions(self._structure)
--> 591 el_charge_lst = self._structure.charge
    592 el_lst = self._structure.get_chemical_symbols()
    593 el_alphabet_dict = {}

File /u/system/SLES12/soft/pyiron/dev/anaconda3/lib/python3.8/site-packages/pyiron_atomistics/atomistics/structure/atoms.py:2155, in Atoms.__getattr__(self, item)
   2153 if item in self._tag_list.keys():
   2154     return self._tag_list._lists[item]
-> 2155 return object.__getattribute__(self, item)

AttributeError: 'Atoms' object has no attribute 'charge'

[sudarsan-surendralal]

Are you still using the default pyiron version installed on the cluster? If so, the changes I made would not be reflected in the pyiron source code since it has not yet been merged and updated. You could wait till I do this

or..

1. Go to a directory of your choice on cmmc and clone pyiron_atomistics (git clone https://github.com/pyiron/pyiron_atomistics.git)
2. Checkout the general_bonds branch (git checkout "general_bonds")
3. Then add this directory to the '$PYTHONPATH' environment variable in the .bashrc file
4. Restart your jupyter server, open your notebook and: import pyiron_atomistics; print(pyiron_atomistics.__file__) to check that this points to your cloned directory. Then try running the notebook.

[ahmedabdelkawy]

Lammps now does not automatically generate the bonds in the structure.inp file. However, when I try to enforce the bond definition on any of the species I have, it does not get reflected in the structure.inp file.

[sudarsan-surendralal]

I guess in atom style charge you won't be able to define bonds. Could you post how your ideal structure.inp should look like here? Or at least how you'd like the bonds and charges to be set. Unfortunately, it is very difficult to generalize this unless we have concrete examples.