Skip to content

Sulfur atom typing incomplete for sp3 hybridization #63

@kbsezginel

Description

@kbsezginel

See the trace for UFF4MOF:

Traceback (most recent call last):
  File "/Users/kutay.sezginel/anaconda3/envs/flex/bin/lammps-interface", line 33, in <module>
    sys.exit(load_entry_point('lammps-interface', 'console_scripts', 'lammps-interface')())
  File "/Users/kutay.sezginel/kutay/git/lammps_interface/lammps_interface/cli.py", line 17, in main
    sim.assign_force_fields()
  File "/Users/kutay.sezginel/kutay/git/lammps_interface/lammps_interface/lammps_main.py", line 441, in assign_force_fields
    param = getattr(ForceFields, self.options.force_field)(**attr)
  File "/Users/kutay.sezginel/kutay/git/lammps_interface/lammps_interface/ForceFields.py", line 3144, in __init__
    self.compute_force_field_terms()
  File "/Users/kutay.sezginel/kutay/git/lammps_interface/lammps_interface/ForceFields.py", line 53, in compute_force_field_terms
    self.compute_atomic_pair_terms()
  File "/Users/kutay.sezginel/kutay/git/lammps_interface/lammps_interface/ForceFields.py", line 62, in compute_atomic_pair_terms
    self.pair_terms(n, data, self.cutoff, charges=charges)
  File "/Users/kutay.sezginel/kutay/git/lammps_interface/lammps_interface/ForceFields.py", line 3152, in pair_terms
    data['pair_potential'].eps = UFF4MOF_DATA[data['force_field_type']][3]
KeyError: 'S_3'

There are 3 different S_3 atom types for UFF4MOF:

    "S_3+2": (1.064, 92.1, 4.035, 0.274, 13.969, 2.703, 0.484, 1.25, 6.928, 4.486, 1.047),
    "S_3+4": (1.049, 103.2, 4.035, 0.274, 13.969, 2.703, 0.484, 1.25, 6.928, 4.486, 1.047),
    "S_3+6": (1.027, 109.47, 4.035, 0.274, 13.969, 2.703, 0.484, 1.25, 6.928, 4.486, 1.047),

which correspond to the following geometries

S_3+2 ch3sch3 Bent two-coordinated sulfur (sp3), two single bonds
S_3+4 socl2 Pyramidal three-coordinated hypervalent sulfur
S_3+6 so2cl2 Tetrahedral four-coordinated hypervalent sulfur

It seems like in UFF S_3+6 is selected by default, I am not sure why that's the case. I think both of these force fields should have a special case for the sulfur atom. Would something like this work? (Assuming lone pair electrons are not counted as neighbors)

if data['element'] == "S":
    if self.graph.degree(node) == 4:
        data['force_field_type'] = "S_3+6"
    elif self.graph.degree(node) == 3:
        data['force_field_type'] = "S_3+4"
    elif self.graph.degree(node) == 2:
        data['force_field_type'] = "S_3+2"

I would be willing to make that change and open a PR if you think the solution above makes sense.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions