Accessing per-atom quantities using umbrella code

LAMMPS LAMMPS Mailing List Mirror
1

Hello, I am running the umbrella code examples/COUPLE/simple code using a C code driver. I am trying to understand how to extract computed per-atom quantities, specifically the potential energy of any given atom. Below is the in.lj file from the example with the compute line added for potential energy:

# 3d Lennard-Jones melt

units lj
atom_style atomic
atom_modify map array

lattice fcc 0.8442
region box block 0 4 0 4 0 4
create_box 1 box
create_atoms 1 box
mass 1 1.0

velocity all create 1.44 87287 loop geom

pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5

compute peratom all pe/atom

neighbor 0.3 bin
neigh_modify delay 0 every 20 check no

fix 1 all nve

variable fx atom fx

dump 1 all atom 1 dump.lj

run 100

I added the following to the C code, but am getting errors. Which command should I use to get the c_peratom values?

double *u = (double *) lammps_extract_atom(lmp,“c_peratom”);
printf(“PE on 1 atom via extract_variable: %g\n”,u[0]);

Thanks,

Anne

2

Hello, I am running the umbrella code examples/COUPLE/simple code using a C code driver. I am trying to understand how to extract computed per-atom quantities, specifically the potential energy of any given atom. Below is the in.lj file from the example with the compute line added for potential energy:

​[…]​

I added the following to the C code, but am getting errors.

​nobody will be able to provide meaningful help unless you tell us what those errors are! nobody here has a crystal ball or can read minds.​

Which command should I use to get the c_peratom values?

​please have a look at src/library.cpp which has detailed explanations on the multiple library interface APIs provided by LAMMPS.​

​axel.​

3

Thanks Axel. I looked more closely into library.cpp and found “lammps_extract_compute”. According to the file:

void *lammps_extract_compute(void *ptr, char *id, int style, int type)

Where

id = compute ID
---------------------> I am assuming this is “pe” for potential energy

style = 0 for global data, 1 for per-atom data, 2 for local data

----------------------> assuming “1” for per-atom data
type = 0 for scalar, 1 for vector, 2 for arr

----------------------------> assuming “1” for vector

Using this, I added the lammps_extract_compute line in the following block of code in the original simple.c example:

// extract force on single atom two different ways

if (lammps == 1) {
double **f = (double **) lammps_extract_atom(lmp,“f”);
printf(“Force on 1 atom via extract_atom: %g\n”,f[0][0]);

double *fx = (double *) lammps_extract_variable(lmp,“fx”,“all”);
printf(“Force on 1 atom via extract_variable: %g\n”,fx[0]);

double *ptr = (double ) lammps_extract_compute(lmp,(char ) “pe”,1,1);
**printf(“POTENTIAL ENERGY on 1 atom via extract_variable: %f~~~~~~~~~\n”,ptr[1]);

}

After trying to run this, I get the following error:

**Force on 1 atom via extract_atom: 21.0465** **Force on 1 atom via extract_variable: 21.0465** __[anne-HP-Notebook:20852] *** Process received signal ***__ **[anne-HP-Notebook:20852] Signal: Segmentation fault (11)** **[anne-HP-Notebook:20852] Signal code: Address not mapped (1)** **[anne-HP-Notebook:20852] Failing at address: 0x8** **[anne-HP-Notebook:20852] [ 0] /lib/x86_64-linux-gnu/libpthread.so.0(+0x11670)[0x7faf2fd78670]** **[anne-HP-Notebook:20852] [ 1] simpleC(+0x673be)[0x56315d93d3be]** **[anne-HP-Notebook:20852] [ 2] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf1)[0x7faf2ee293f1]** **[anne-HP-Notebook:20852] [ 3] simpleC(+0x66dea)[0x56315d93cdea]** __[anne-HP-Notebook:20852] *** End of error message ***__
4

Thanks Axel. I looked more closely into library.cpp and found “lammps_extract_compute”. According to the file:

void *lammps_extract_compute(void *ptr, char *id, int style, int type)

Where

id = compute ID
---------------------> I am assuming this is “pe” for potential energy

​wrong. it is the ID of the compute you created. in your posted input example case it would be “peratom”. see the docs for the compute command.

style = 0 for global data, 1 for per-atom data, 2 for local data

----------------------> assuming “1” for per-atom data

​correct.​

type = 0 for scalar, 1 for vector, 2 for arr

----------------------------> assuming “1” for vector

​correct.
​

Using this, I added the lammps_extract_compute line in the following block of code in the original simple.c example:

// extract force on single atom two different ways

if (lammps == 1) {
double **f = (double **) lammps_extract_atom(lmp,“f”);
printf(“Force on 1 atom via extract_atom: %g\n”,f[0][0]);

double *fx = (double *) lammps_extract_variable(lmp,“fx”,“all”);
printf(“Force on 1 atom via extract_variable: %g\n”,fx[0]);

double *ptr = (double ) lammps_extract_compute(lmp,(char ) “pe”,1,1);
**printf(“POTENTIAL ENERGY on 1 atom via extract_variable: %f~~~~~~~~~\n”,ptr[1]);

}

unfortunately, ​it is not as simple as this. for efficiency reasons, computations of the potential energy are not performed explicitly, but done as part of the force computation. but for that to be triggered, there has to be a “consumer” for the compute defined (e.g. a dump style or a fix), so that LAMMPS can determine that it needs to accumulate per-atom potential energy. so to access the data of a compute between runs, you need to do something like this:

if (lammps == 1) {
int i,j;
int *nlocal;
double *pe;
int *id = lammps_extract_atom(lmp,“id”);
lammps_command(lmp,“fix xxx all ave/atom 1 1 1 c_peratom”);
lammps_command(lmp,“run 0”);

nlocal = lammps_extract_global(lmp,“nlocal”);
pe = lammps_extract_compute(lmp,“peratom”,1,1);
for (i = 0; i < nprocs_lammps; ++i) {
MPI_Barrier(comm_lammps);
if (i == me) {
for (j = 0; j < *nlocal; ++j)
printf(“PE of atom %d: %g\n”,id[j],pe[j]);
}
}
lammps_command(lmp,“unfix xxx”);
}

​please also note, that per-atom data is distributed across processors (hence the loops over nlocal and mpi ranks). only for a single MPI rank, do you have “nlocal == natoms”.

axel.​

5

Great, thank you Axel. I used your code and it worked on one processor.

Regarding your last note, if I was to use more than one processor, what modifications would I need to make to account for per atom data being distributed?

6

Great, thank you Axel. I used your code and it worked on one processor.

Regarding your last note, if I was to use more than one processor, what modifications would I need to make to account for per atom data being distributed?

​my code fragment already handles that case. the issue would be with your calling code, about which i know nothing. if you look carefully, you’ll see that “nlocal”, “id”, and “pe” are all local per MPI rank, and thus will contain different data, when running LAMMPS across multiple processors. please also note, that my code fragment is to illustrate the situation, it is not the best way to do this kind of thing in a setting where efficiency matters.​

​axel.​

7

Hi Axel. For my calling code, I just need to extract the potential energy of one atom of a known ID in a lattice. For example, the 36th atom. I understand the code snippet and how I could write some “if” statements to accomplish this, but that may not be the most computationally-efficient method.

What would be your suggestion be to most efficiently extract the PE of a single atom of known ID? I am expecting the amount of processors to vary as I develop, test, and use my code; is there any way to force the potential energy PE of a given atom ID to be calculated on a particular rank?

Thanks,

Anne

8

Hi Axel. For my calling code, I just need to extract the potential energy of one atom of a known ID in a lattice. For example, the 36th atom. I understand the code snippet and how I could write some “if” statements to accomplish this, but that may not be the most computationally-efficient method.

What would be your suggestion be to most efficiently extract the PE of a single atom of known ID?

​you need to avoid the extra “run”, so put a “consumer” of the per-atom compute into your regular run, and have it access the compute only on the steps, that you need the data on.​

I am expecting the amount of processors to vary as I develop, test, and use my code; is there any way to force the potential energy PE of a given atom ID to be calculated on a particular rank?

​no.​ you can only ask to tally the per-atom potential energy for either all atoms or none. …and that data is available on whatever MPI rank owns that atom. if you want that information on a different MPI rank, you need to use MPI communication.

axel.

9

Hi Axel,

What exactly is a “consumer”; could you please provide a simple example line?

Is there any way to execute the compute with zero timesteps? I just need the potential energy, without even needing to create any timesteps. I am relying on LAMMPS to evaluate an EAM potential of a single atom in a fixed configuration that my driver code uses for its own computations.

Thank you,

Anne

10

Hi Axel,

What exactly is a “consumer”; could you please provide a simple example line?

​see the use of fix ave/atom.

to repeat, the potential energy compute is different from many other computes, as it does not compute anything by itself, but only gives access to data that is already tallied during a force computation. you cannot change that.
the consequence is, that something, e.g. a fix or a dump have to reference this compute to tell LAMMPS before the force computation, that the per-atom potential energy is required. otherwise you get an error of the kind “energy not tallied on this time step”. (the same applied to per-atom stress, BTW).

the “run 0” command, that follows the definition of the fix, does exactly what you are asking for. it will only run the “setup phase” of an MD run, where the forces are reconstructed without actually advancing the positions.

Is there any way to execute the compute with zero timesteps? I just need the potential energy, without even needing to create any timesteps. I am relying on LAMMPS to evaluate an EAM potential of a single atom in a fixed configuration that my driver code uses for its own computations.

since the potential energy has contributions from multiple atoms, it is not possible to separate those. you will have to run a complete force computation (or write your own EAM code).

axel.​

11

Ok, I got it. Thanks for all of your help with this.

12

Here is my advice.

Step 1: write an input script that calculates
the per-atom energy of the single atom
you want (ID = 36) and prints it out with
thermo output. You could use compute pe/atom
with group = all and just output it with
thermo via a variable. This should work
with a “run 0”.

Step 2: now reproduce all those lines
in your calling code to do the same thing
thru the library interface. Define the compute
and the variable, etc. Invoke the run 0,
then invoke the variable and retrieve its
value.

Having both step 1 and 2 lets you compare
the results to insure you are doing the right thing.

Steve

13

Thanks Steve. Step 2 describes what I have been wanting to do. However, when I tried to implement it in the code, I got a bunch of errors, as seen in my previous emails. Currently I am relying on the code Axel provided. Would you be able to send example lines of code on how to implement Step 2? It would be greatly appreciated.

Anne

14

Thanks Steve. Step 2 describes what I have been wanting to do. However, when I tried to implement it in the code, I got a bunch of errors, as seen in my previous emails. Currently I am relying on the code Axel provided. Would you be able to send example lines of code on how to implement Step 2? It would be greatly appreciated.

​the sample code in ​examples/COUPLE/simple/simple.c already shows how to access data from a per-atom variable.

axel.

15

Hi Axel, I looked through simple.c several times, but I am only seeing how to access coordinates, velocities, and forces. However, I am specifically needing potential energy. From previous emails, it is apparent that potential energy has to be accessed in a different way, gathered from each of the processors and looping through each atom id.

Steve’s email says “Step 2: now reproduce all those lines
in your calling code to do the same thing
thru the library interface. Define the compute
and the variable, etc. Invoke the run 0,
then invoke the variable and retrieve its
value.” … this sounds like there is some direct way of accessing the potential energy of a particular atom id without looping through all of them. I am interested if there is a code that would make this possible.

Thanks
Anne

16

this sounds like there is some direct way of accessing the potential energy of a particular atom id without looping through all of them. I am interested if >there is a code that would make this possible.

Define an equal-style variable something like this:

compute eatom all pe/atom
variable one equal c_eatom[36]
thermo_style custom step temp pe v_one
run 0
print "Eatom = {one}"

Do it from a script first, then from the library interface.

Steve

17

Thank you Steve. I added those lines to in.lj in the COUPLE example, and am running just in LAMMPS without the driver code. I am having the following 2 issues, and a question:

  1. Whether I comment out thermo_style or not, I do not see any difference in the printed thermodynamic output in the terminal. This is what I get in either situation:
    Step Temp E_pair E_mol TotEng Press
    0 1.44 -6.7733681 0 -4.6218056 -5.0244179

  2. When I include the “print” line (print "Eatom = {one}"), I get the following error:
    ERROR: Compute used in variable between runs is not current (…/variable.cpp:1447)
    Last command: Eatom = ${one

  3. Assuming everything in the LAMMPS script will work properly, would this be the correct line in the C driver code to get the variable “one”?

double *pe = (double *) lammps_extract_variable(lmp,“one”,“all”);
printf(“PE on atom 36 via extract_variable: %g\n”,pe);

For reference, this is the in.lj file that I am using. I tried varying the order of commands around, but I keep getting the same error.

# 3d Lennard-Jones melt

units lj

atom_style atomic

atom_modify map array

lattice fcc 0.8442

region box block 0 4 0 4 0 4

create_box 1 box

create_atoms 1 box

mass 1 1.0

velocity all create 1.44 87287 loop geom

pair_style lj/cut 2.5

pair_coeff 1 1 1.0 1.0 2.5

compute eatom all pe/atom #added to example

neighbor 0.3 bin

neigh_modify delay 0 every 20 check no

fix 1 all nve

run 0 #added to example

variable fx atom fx

variable one equal c_eatom[36] #added to example

thermo_style custom step temp pe v_one #added to example

print "Eatom = {one}" #added to example

dump 1 all atom 1 dump.lj

18

anne,

you seem to be doing step 2 without having done step 1.
to make sense of all this, you need to understand the the affected commands. i.e. read the corresponding pages in the documentation and understand what the examples you are looking at or have been given do.

we can point you into the right direction, but we don’t have the time to give you all the details.

Thank you Steve. I added those lines to in.lj in the COUPLE example, and am running just in LAMMPS without the driver code. I am having the following 2 issues, and a question:

  1. Whether I comment out thermo_style or not, I do not see any difference in the printed thermodynamic output in the terminal. This is what I get in either situation:
    Step Temp E_pair E_mol TotEng Press
    0 1.44 -6.7733681 0 -4.6218056 -5.0244179

  2. When I include the “print” line (print "Eatom = {one}"), I get the following error:
    ERROR: Compute used in variable between runs is not current (…/variable.cpp:1447)
    Last command: Eatom = ${one

​yes, because you didn’t correctly configure a “consumer” for the compute (a custom thermo style in this case).

  1. Assuming everything in the LAMMPS script will work properly, would this be the correct line in the C driver code to get the variable “one”?

double *pe = (double *) lammps_extract_variable(lmp,“one”,“all”);
printf(“PE on atom 36 via extract_variable: %g\n”,pe);

​as a matter of principle, we don’t answer “what if” style questions (aka “re-insurance” questions). try it out, read the documentation (and code/example comments), study and de​bug it yourself.

​please note, that we have been rather generous in your help, since you seem to be new to this. however, you have to pick up the slack yourself and learn from what you have been told and make your own conclusions or else you will be get less and less help and at the end nothing.​

​axel.​

19

My example lines included a thermo_style custom that output the
value of the “one” variable. That is what Axel means
by a “consumer” of the compute. You seem to be using the
default thermo output. Hence the error.

Steve