How to calculate energies spherically around each molecule in the system

Hi everyone,

I’m working on problem that I was hoping someone might have a suggestion for.

For each molecule in my system (a lipid bilayer, but the system info isn’t super important to the question) I need to create three groups:

  1. the molecule itself,
  2. atoms within a spherical cutoff around that molecule, like 10 angstroms (but are not part of that central molecule), and
  3. all other molecules.

I then want to use compute group/group to get all of the energies between each of these groups. I need to repeat this process for each molecule in the system (basically treating each molecule separately as a central molecule), and for every configuration of the system.

The question I have for the list is the following: is there an easy way of doing this? One option would be to calculate the distances separately, divide into groups, and then write a new data file ordered by distance from the molecule, and then do the group command by atom id, and just run many energy calculations. That would create a lot of files, so I was hoping someone might have a suggestion for a more efficient suggestion for defining these groups that relies on lammps commands rather than rewriting the data file many times.

Thanks in advance!

Compute group/group is a very expensive operation, especially when you have long-range electrostatics.

I would think about doing this in post-processing, where you can easily substitute long-range coulomb with a cutoff coulomb (your computation is going to be rather approximate anyway and you follow an accurately computed trajectory!), then you could consider using compute pe/mol/tally to get the inter-molecule contributions to the potential energy. That is not exactly what you describe, but since this is going to be a semi-quantitative property (it will depend significantly on the choice of cutoff) I think it would be still representative and effective to implement with available facilities in LAMMPS. With this you would still have to do one rerun per molecule. However, that can be easily parallelized. As an alternative, you could create a per-atom variant of the compute and then use compute chunk/atom and compute reduce/chunk to get everything in one go.

Anything more complex would then probably be best

1 Like

Axel -

Thank you for your suggestion - I’ll definitely take a look at compute pe/mol/tally, because that looks like it would definitely help with what I am trying to do. I agree that parallelizing and rerunning are probably the best way of doing this, I just wanted to make sure before I went through and did that.

Thank you!