Small modification of hybrid/molecular

LAMMPS LAMMPS Development
1

Dear Lammps developers,

I am working on some polymers which yield bad behaviors from a pairwise potential (e.g. OPLS) with a special_bond exclusion. In a short word, I need to exclude/soften some 1-6 and 1-7 interactions to get correct energy barrier of some dihedral’s rotation compared with DFT outcome. However, since I still wish to keep other non-bond interactions, a very straight-forward idea is using Molecule ID of each atoms as a cutoff indicator, which brings me to modify pair hybrid/molecular to make it could define “INTRA” and “INTER” by a threshold.

To fulfill the purpose, I made following modification on lammps version 29Aug2024:

in pair_hybrid_molecular.cpp and .h:

1,Copy/paste PairHybridMolecular::settings from PairHybrid::settings and made it could read an additional integer, e.g. pair_style hybrid/molecular 0 lj/cut 15.0 lj/cut 15.0. The added integer were stored as a variable named “MolecularIDCutoff”.

2,in PairHybridMolecular::init_style()
Add a function to pass the stored variable “MolecularIDCutoff” to class neigh_request.
The set->molskip(NeighRequest::INTRA) and (NeighRequest::INTER) were changed to (NeighRequest::INTRA_P) and (NeighRequest::INTER_P) respectively, which indicate the new molecule ID cutoff neighbor list building strategy.
Changes are made accordingly in neigh_request.cpp and .h

3,In PairHybridMolecular::single()
Change “if (atom->molecule[i] != atom->molecule[j]) m = 1” to “if (std::abs(atom->molecule[i] -atom->molecule[j])>MolecularIDCutoff) m = 1”

in neigh_list.cpp and .h:
4, in NeighList::post_constructor(), changes are made to bypass MolecularIDCutoff to NeighList.

in npair_skip.cpp and .h:
5, in NPairSkipTemp::build(), changes are made to handle (NeighRequest::INTRA_P) and (NeighRequest::INTER_P) after handdling (NeighRequest::INTER) as followed:
if ((molskip == NeighRequest::INTRA_P) && (abs(molecule[i] - molecule[j]) > MolecularIDCutoff )) continue;
if ((molskip == NeighRequest::INTER_P) && (abs(molecule[i] - molecule[j]) <= MolecularIDCutoff )) continue;

The modified code was compiled and I run some tests on the modified lmp executable. The test system has 2 atoms with 0.5 A distance between.

The intra pair_coeff is “lj/cut 1 1 0.1”
The inter pair_coeff is “lj/cut 2 100 1”

With pair_style hybrid/molecular 0 lj/cut 15.0 lj/cut 15.0, MID = 2,3 → strong repulsion
With pair_style hybrid/molecular 1 lj/cut 15.0 lj/cut 15.0, MID = 2,3 → weak attraction
With pair_style hybrid/molecular 2 lj/cut 15.0 lj/cut 15.0, MID = 2,3 → weak attraction
With pair_style hybrid/molecular 0 lj/cut 15.0 lj/cut 15.0, MID = 2,2 → weak attraction
With pair_style hybrid/molecular 0 lj/cut 15.0 lj/cut 15.0, MID = 2,3 → strong repulsion

It seems like the behavior is correct. However, I have following questions:
1, Did I made all the required modification?
2, What sort of tests I need to conduct to make sure I did not break anything?
3, In PairHybridMolecular::single(), there is an if statement to control the variable m by comparing the molecularID. However, I find that the statement did not affect the simulation outcome. In the beginning, I unintentionally put a wrong condition with this particular if statement but the result seems correct (at least for potential energy). I was quite confused about it.

Thanks, Best regards.
Jiawei Zhao, BRICI, Beijing, China

2

Your changes are quite invasive, so it is not likely that someone will even want to look at it from such a vague description that is not a ready-to-compile git branch. Even then it is quite unlikely, since the modification is in your interest and not anybody else’s and thus you have the most incentive to validate it. This doesn’t look like a change that could be added to the LAMMPS distribution.

LAMMPS comes with a large library of integrated tests. You should verify you can run those successfully with an unmodified version and then check out your modified version. Also, there are plenty of examples with reference outputs. You could use the python scripts in tools/regression-tests for running those in an automated fashion. Please note that this is not yet fully compatible: [Issue] Errors reported by regression tests · Issue #4347 · lammps/lammps · GitHub

During a regular simulation, the Pair::single() method is not called. It is used by certain functionality in LAMMPS that needs to compute the force/r and energy for a single pair of atoms, e.g. the pair_write command or compute group/group.

3

Dear Axel.

Thanks for your reply, suggestion and answer. I will check myself out with the regression test tool as you mentioned.

Jiawei Zhao

4

This is only the second best option. The best is the aforementioned integrated testing. See: 3.11. Development build options — LAMMPS documentation

5

If you are mainly concerned about imposing correct dihedral-dihedral interactions, you may wish to look at the CMAP potential frequently used in protein simulations for that precise purpose, implemented in LAMMPS as fix cmap.