(many)body potentials in lammps: splitting a bulk system into 1,2,3 bodies?

I am totally new to lammps, could you kindly help me with (many)body potentials?

What i want to do is as follows:

Let us suppose i have a simulation box of bulk water that i want to use for MD simulations.

I want to define all possible combinations of dimers and trimers within this box, and apply a 2-body and 3-body potential, i.e., during MD i want to calculate the potential energy of my system as a sum of contributions 1body+2body+3body, where a body is defined molecularwise (a body = 1H2O) and a body is not rigid (the structure of H2O is allowed to change during MD).

Is it possible to use BODY or MANYBODY package (or any other) for MD considering the aforementioned desired conditions, i.e.,

1.splitting a bulk system into 1,2,3 bodies.

2.calc pot ene as a sum of 1,2,3 body contributions.

3.potential that allows the structure of H2O to be changed (not rigid).

I would be grateful for your help.

Yours sincerely,

Yevhen Horbatenko

Research Professor,

Room 219, IBS CMSD,

R&D Center at Korea University,

Anam-ro 145, Seoul, Republic of Korea

Tel: + 82-10-2779-5440

E-mail: papilio.p…@gmail.com

None of the existing potentials in LAMMPS operates like this, they are all based on interactions of individual atoms. You would have to implement your model as a new potential from scratch in C++.

Since you’re dealing specifically with water, I wonder if you could decompose / rewrite your model into a combination of intramolecular forces (which are then normal bond / angle potentials) and a many-body force between the oxygens. After all oxygen is close enough to the centre of mass of water that you should be able to approximate a whole-molecule force as a force on the oxygens.

I see.
Thank you so much for your reply.

Thank you for the suggestion.
I have never thought this way.
Let me think if can do that for my model.

Could you kindly indicate the place (.cpp and a function within it) where the reading of xyz coordinates occurs from input within lammps? Also, the place where the coordinates get updated during MD?
It would be very helpful and prevent me get lost among source files.

There is not a single place for that. Coordinates can be input in different ways with different commands.

Again, there is not a single place for that. There are multiple ways how positions are updated and there is a bit of a hierarchy of how this is addressed.

Please also note that LAMMPS is highly modular and takes advantage of polymorphism in C++ to have common procedures that can be overridden and thus the flow of control is rarely confined to a single file, but would climb up and down the class hierarchy. You need to develop a “big picture” view of the overall flow and then focus on the details.
For example, to implement a pair style that computes forces based on atom positions from a neighbor list, you can focus on the computation of the forces, since almost everything else will be transparently provided. Of course, the more you stray from the common workflow, the more you will have to augment those procedures with custom code.

Have you run MD simulations in general before? And simulations with LAMMPS?
If not, you should become familiar with the basic workflows of LAMMPS to conduct molecular simulations before you should start programming. Otherwise you will be lost all them time.

Then you can get some orientation by reading through 4. Information for Developers — LAMMPS documentation and 3. Modifying & extending LAMMPS — LAMMPS documentation

Thank you for your suggestions and for pointing out that LAMMPS uses the polymorphism of C++.
It was not evident to me.

Answering your questions:
I have done MD before, although with different software (GAMESS-US, VASP). I have the understanding of how MD usually works.
Simulations with LAMMPS, I have never tried them before. This is the first time I need LAMMPS, although I have done a couple of examples among those LAMMPS provides.

If you allow me to be more specific regarding the place where xyz corrds are read,
I need the place that reads xyz corrds from a data file which is specified in an input file with read_data keyword.

For example, there is in.gap file that contains the following line “read_data data_gap”.
There is also data_gap file that has xyz coords.
If someone could indicate me the place within the code that reads data_gap file, it would be great.
At least i will have the place to start from, while becoming more familiar with the basic LAMMPS workflows.

The read_data function has its source code in the src/read_data.cpp and src/read_data.h files. Sadly, we are not very creative in how we name our source files.

Studying the read_data command is probably the worst way to start learning about the flow of control in LAMMPS. It is one of the more convoluted and badly designed pieces of LAMMPS (mostly for legacy reasons and how it was later extended to do things it was never designed for in the first place and because its underlying file format is badly designed as well. Well, “designed” is probably not the right word, like so many legacy pieces of software that date back such a long time).

Based on your description, the best way to start is to read the recent LAMMPS paper and then continue with the (entire!) parts of the documentation I already pointed out to you. I don’t think that reading source code will help much, but if you put me on the spot, I would suggest reading the PairLJCut::compute() function and the PairSW::compute() function in the src/pair_lj_cut.cpp file and the src/MANYBODY/pair_sw.cpp file, respectively.

Thank you so much!

“we are not very creative in how we name”
From the other hand it is very informative. One does not need to spend much time if one wants to find the same info I was asking :slight_smile:

Thanks a lot for your explanation and suggestions!
Let me follow them as well.