Lammps constraints

[mattwthompson]

Description

Provide a brief description of the PR's purpose here.

Checklist

• Add tests
• Lint
• Update docstrings

[codecov]

Codecov Report

Merging #897 (9258d6b) into main (e7f0298) will increase coverage by 0.02%.
The diff coverage is 100.00%.

❗ Current head 9258d6b differs from pull request most recent head 65dfd80. Consider uploading reports for the commit 65dfd80 to get more accurate results

Additional details and impacted files

[mattwthompson]

I think I'm going about this wrong. I'm trying to use the fix shake command to constrain all hydrogen-to-heavy-atom bonds, but it's choking on even the simplest case of methane. Clearly SHAKE isn't the algorithm that's meant to be used. Here are the files and logs:

out.lmp:

Title

5 atoms
4 bonds
6 angles
0 dihedrals
0 impropers

2 atom types
1 bond types
1 angle types
-1.00390625 38.99609375 xlo xhi
-0.4235839844 49.57641602 ylo yhi
-0.8071289062 59.19287109 zlo zhi
0.0 0.0 0.0 xy xz yz

Masses

1	12.01078
2	1.007947


Pair Coeffs

1	0.10884061	3.3795318
2	0.015779483	2.6445434

Bond Coeffs

1 harmonic	359.82124644905	1.090139506109


Angle Coeffs

1 harmonic	37.541272178735	108.5839257083


Atoms

1	1	1	-0.10868	-6.5803528e-05	-6.1392784e-06	2.1517277e-05
2	1	2	0.027170001	-0.056671143	1.0869141	-0.0859375
3	1	2	0.027170001	0.61962891	-0.3972168	-0.80712891
4	1	2	0.027170001	-1.0039062	-0.42358398	-0.069580078
5	1	2	0.027170001	0.44165039	-0.26660156	0.96289062

Bonds

1	1	1	2
2	1	1	3
3	1	1	4
4	1	1	5

Angles

1	1	4	1	5
2	1	3	1	4
3	1	2	1	4
4	1	2	1	5
5	1	2	1	3
6	1	3	1	5

tmp.in:

units real
atom_style full

dimension 3
boundary p p p

bond_style hybrid harmonic
angle_style hybrid harmonic
dihedral_style hybrid fourier
improper_style cvff
special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.8333333333 
pair_style lj/cut/coul/long 9.0 9.0
pair_modify mix arithmetic tail yes

read_data out.lmp

thermo_style custom ebond eangle edihed eimp epair evdwl ecoul elong etail pe

fix 100 all shake 0.0001 20 10 b 1
kspace_style pppm 1e-4
run 0

log.lammps:

LAMMPS (21 Nov 2023)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
  using 1 OpenMP thread(s) per MPI task
units real
atom_style full

dimension 3
boundary p p p

bond_style hybrid harmonic
angle_style hybrid harmonic
dihedral_style hybrid fourier
improper_style cvff
special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.8333333333
pair_style lj/cut/coul/long 9.0 9.0
pair_modify mix arithmetic tail yes

read_data out.lmp
Reading data file ...
  triclinic box = (-1.0039062 -0.42358398 -0.80712891) to (38.996094 49.576416 59.192871) with tilt (0 0 0)
  1 by 1 by 1 MPI processor grid
  reading atoms ...
  5 atoms
  scanning bonds ...
  4 = max bonds/atom
  scanning angles ...
  6 = max angles/atom
  reading bonds ...
  4 bonds
  reading angles ...
  6 angles
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0.5     
  special bond factors coul:  0        0        0.8333333333
     4 = max # of 1-2 neighbors
     3 = max # of 1-3 neighbors
     3 = max # of 1-4 neighbors
     4 = max # of special neighbors
  special bonds CPU = 0.000 seconds
  read_data CPU = 0.014 seconds

thermo_style custom ebond eangle edihed eimp epair evdwl ecoul elong etail pe

fix 100 all shake 0.0001 20 10 b 1
Finding SHAKE clusters ...
ERROR: Shake cluster of more than 4 atoms (src/RIGID/fix_shake.cpp:978)
Last command: fix 100 all shake 0.0001 20 10 b 1

[mrshirts]

Yeah, it looks like it's trying to treat them all as the same cluster, which is not what one wants to do generally.

Can you look at: https://docs.lammps.org/molecule.html in the "shake atoms" and "shake bond types" section?

I think methane is actually a special case that isn't allowed, since there are 4 hydrogen bonds attached to one heavy atom. Try ethane, which should be supported.

We need a LAMMPS expert on this . . .

[mattwthompson]

That's precisely it! (And, now that my memory is jogged like a student looking at the answer key, I remember this isn't the first time I've come across this corner case: openmm/openmm#3856) The idea of methane is that, since all of its bonds are constrained, the total bond energy should be zero. Ammonia handles that well enough, though.

I haven't implemented rigid water yet, but this is looking okay for some simple molecules:

In [16]: def run_lammps(smiles: str) -> EnergyReport:
    ...:     molecule = MoleculeWithConformer.from_smiles(smiles)
    ...:     molecule.name = "MOL"
    ...:     topology = molecule.to_topology()
    ...:     topology.box_vectors = Quantity([4, 4, 4], "nanometer")
    ...:
    ...:     return smiles, get_lammps_energies(sage.create_interchange(topology)).energies

In [17]: [run_lammps(smiles) for smiles in ["N", "CCO", "CC(=O)OC1=CC=CC=C1C(=O)O"]]
Out[17]:
[('N',
  {'Bond': 0.0 <Unit('kilocalorie_per_mole')>,
   'Angle': 0.62354148 <Unit('kilocalorie_per_mole')>,
   'Torsion': 0.0 <Unit('kilocalorie_per_mole')>,
   'vdW': -7.6032464e-05 <Unit('kilocalorie_per_mole')>,
   'Electrostatics': -0.01761300000000432 <Unit('kilocalorie_per_mole')>}),
 ('CCO',
  {'Bond': 0.045460236 <Unit('kilocalorie_per_mole')>,
   'Angle': 1.1155101 <Unit('kilocalorie_per_mole')>,
   'Torsion': 1.4631677 <Unit('kilocalorie_per_mole')>,
   'vdW': 0.15680472977 <Unit('kilocalorie_per_mole')>,
   'Electrostatics': -4.220942000000001 <Unit('kilocalorie_per_mole')>}),
 ('CC(=O)OC1=CC=CC=C1C(=O)O',
  {'Bond': 1.6112717 <Unit('kilocalorie_per_mole')>,
   'Angle': 4.8003716 <Unit('kilocalorie_per_mole')>,
   'Torsion': 11.484984206976 <Unit('kilocalorie_per_mole')>,
   'vdW': 9.6975537526 <Unit('kilocalorie_per_mole')>,
   'Electrostatics': -65.17742600000001 <Unit('kilocalorie_per_mole')>})]