[lammps-users] extending of compute centroid/stress/atom command

LAMMPS LAMMPS Mailing List Mirror
1

Dear all,

I am going to use compute centroid/stress/atom command. I found this in the manual: Note that the stress for each atom is due to its interaction with all other atoms in the simulation, not just with other atoms in the group. So, I wonder how can I just calculate the interaction with other atoms in the group.

My pair style is hybrid (lj/cut and lj/charmm/coul/charmm). There are three kinds of atoms in the box:

1-1 lj/charmm/coul/charmm

1-2 lj/cut

1-3 lj/cut

2-2 lj/cut

3-3 lj/cut

Then： group A type 2 3.

I only want to calculate the pair interaction of group A (2-2 and 3-3). I have read some source codes of compute_centroid_stress_atom. And I found:

if (force->pair->centroidstressflag == CENTROID_AVAIL) {

double **cvatom = force->pair->cvatom;

for (i = 0; i < npair; i++)

for (j = 0; j < 9; j++)

stress[i][j] += cvatom[i][j];

}

Then, I only found the zero accumulators of cvatom in pair.cpp:

if (cvflag_atom && alloc) {

n = atom->nlocal;

if (force->newton) n += atom->nghost;

for (i = 0; i < n; i++) {

cvatom[i][0] = 0.0;

…

but I don’t find how cvatom is calculated. Also, I don’t find cvatom in pair_lj_cut.cpp.

This is my first attempt to modify the LAMMPS code. I would like to know where and how to achieve this feature. As the final goal is to calculate the auto autocorrelation function of heat flux, I think it’s the best way to change the source code, rather than pre- or post-processing

Thank you so much for your help.

[

2296804970

[email protected]

](https://maas.mail.163.com/dashi-web-extend/html/proSignature.html?ftlId=1&name=2296804970&uid=2296804970%40qq.com&iconUrl=https%3A%2F%2Fmail-online.nosdn.127.net%2Fqiyelogo%2FdefaultAvatar.png&items=["2296804970%40qq.com"])

签名由 网易邮箱大师 定制

2

There are multiple things here that don’t make much sense:

  • if I remember correctly, centroid stress and virial stress are the same for pair-wise additive potentials like you are using.
  • you are looking in the wrong places for the computation of per-atom stress. the compute does not “compute” (that would be too time consuming) instead it signals that it should be accumulated during the regular force calculation. the code in pair.cpp is clearing the per-atom stress arrays before the forces (and stresses and energies) are accumulated. the actual stress (and energy) accumulation is done via calls to one (or more) of the ev_tally() or v_tally() functions within the actual pair style in the compute() method. this can all be inferred from carefully reading the source code and paying attention to the comments plus some of the general explanations of the process and flow of control in the manual and the LAMMPS paper(s) and particularly the papers explaining how the stress in LAMMPS is derived.
  • there is no need to use a hybrid pair style. you can use lj/charmm/coul/charmm throughout. just make certain the charges on atoms of type 2 and 3 are zero

also, do you have some reference in the published literature that the computation you want to do and the property you want to get is actually useful and meaningful?

3

Thank you so much for your help.

-I’m not sure this calculation makes some physical sense. This idea comes from <Heat conduction in chain polymer liquids: Molecular dynamics study>(https://doi.org/10.1063/1.3613648). In my system, there are two kinds of materials: A and B. The only connection between A and B is the pair interaction. I want to eliminate all the effects of B (pair potential energy in Pe/atom, pair force in Stress/atom), and only consider the difference of A between this system and the pure A system. But I couldn’t delete B because this will affect the structure of A in this system.

-There have many-body interactions in my system. So I used compute centroid/stress/atom.

-I have found compute() in lj_cut.cpp.

The reason I used the hybrid pair style is:

I found in pair_hybrid.cpp, vatom[i][j] and eatom[i] are accumulated through different sub-styles:

for (m = 0; m < nstyles; m++) {

…

for (i = 0; i < n; i++) eatom[i] += eatom_substyle[i]

…

for (i = 0; i < n; i++)

for (j = 0; j < 9; j++)

cvatom[i][j] += cvatom_substyle[i][j]

}

So I am thinking about: if I use different sub-styles to distinguish the contact. I could get the specific pe and stress. Like this:

pair_style hybrid lj/cut lj/charmm/coul/charmm lj/cut

pair_coeff 1 1 lj/charmm/coul/charmm

pair_coeff 1 2 lj/cut 2

pair_coeff 1 3 lj/cut 2

pair_coeff 2 2 lj/cut

pair_coeff 3 3 lj/cut

group A type 2 3

compute 1 A centroid/stress/atom NULL virial

compute 2 A pe/atom NULL virial

Then, I change the source of compute() in pair_hybrid.cpp:

for (m = 0; m < nstyles-1; m++) {

…

for (i = 0; i < n; i++) eatom[i] += eatom_substyle[i]

…

for (i = 0; i < n; i++)

for (j = 0; j < 9; j++)

cvatom[i][j] += cvatom_substyle[i][j]

I’m not sure it’s going to work. Or whether there are other effects of this change.

[

2296804970

[email protected]

](https://maas.mail.163.com/dashi-web-extend/html/proSignature.html?ftlId=1&name=2296804970&uid=2296804970%40qq.com&iconUrl=https%3A%2F%2Fmail-online.nosdn.127.net%2Fqiyelogo%2FdefaultAvatar.png&items=["2296804970%40qq.com"])

签名由 网易邮箱大师 定制

On 3/21/2022 01:27，Axel Kohlmeyer[email protected] wrote：

4

hacking pair style hybrid as you describe is about the worst possible approach that I can think of. it will render pair style hybrid unusable for anything else.
the LAMMPS’ way of doing this would be to extend compute pe/atom and/or compute stress/atom and/or compute centroid/stress/atom to have an option to access only the contribution from a hybrid substyle instead of the full (hybrid) pair style data. this could be done with adding some extra settings to the computes via compute_modify.

in the input you would refer to choose the substyle with, e.g.:

compute_modify ID pair lj/cut 2
or
compute_modify ID pair lj/charmm/coul/charmm

that would then augment how the corresponding compute processes the “pair” contribution to the property it computes through accessing only the selected substyle.

inside the compute you would then grab a pointer of the pair hybrid sub-style by calling force->pair_match()
see for example the compute pair command, which provides access to the pair sub-styles’ global energies.

5

Thank you again for your help!

With your help, I almost accomplished my goal.

I wrote three new compute styles (pe/sp/atom, stress/sp/atom, centroid/stress/sp/atom) similar to the old compute styles.

Then in pair_hybird.cpp, I create three new vectors (eatomsp, vatomsp, cvatomsp) like eatom, vatom, cvatom. And make sure that these three new vectors only accumulate under the specific substyles.

At last, in new compute styles, the energy[i] or stress[i][j] points to eatomsp, vatomsp, cvatomsp, instead of pointing to eatom, vatom, cvatom in the old compute styles.

I have examined the impact of these modifications. The temperature, potential energy, pair energy of the whole system is the same before and after modification. Moreover, I tested the energy of atoms using compute group/group. And I found the accumulation of the difference before and after modification is exactly half of compute group/group. So, I believe the modification is accurate.

The only flaw is that when I use compute pe/atom, I found pe/atom and pe/sp/atom show me the same results. But the results of centroid/stress/sp/atom and centroid/stress/atom are different. I don’t know what leads to the wrong calculation of compute pe/atom. But I can overcome this problem by using an unmodified LAMMPS to calculate pe/atom

[

2296804970

[email protected]

](头像签名)

签名由 网易邮箱大师 定制

On 3/21/2022 19:40，Axel Kohlmeyer[email protected] wrote：