peteboyd / lammps_interface

automatic generation of LAMMPS input files for molecular dynamics simulations of MOFs
MIT License
125 stars 63 forks source link

Sulfur atom typing incomplete for sp3 hybridization #63

Open kbsezginel opened 1 year ago

kbsezginel commented 1 year ago

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.