Lammps vector variable calculation

I seem to be having trouble understanding some of the variable computations involving vectors or arrays.

As an example, I have a system of identical molecules. I use

compute molecules all chunk/atom molecule
compute 1 all vcm/chunk molecules

and I now want to find the average translational kinetic energy of the molecules every 100 timesteps. My first attempt was to use

variable vx2 equal 0.5*mass*ave(c_1[1]^2)

to find the average contribution to the kinetic energy in the x-direction. My understanding was that this would square the x-component of the center of mass velocity of each molecule and then take the average of the squares. However this doesn’t seem to be the case. The output of variable vx2 remains near zero and is sometimes positive and sometimes negative. Can someone clarify what is actually being calculated by variable vx2 here?

For an example behaving as I would expect, if I do something similar with the radius of gyration using the following commands

compute rg all gyration/chunk molecules
variable ave_rg equal ave(c_rg)

it seems to output the average radius of gyration at a particular timestep, as I expect.

Any clarification on how the vx2 example works would be appreciated. Or if there is a better way to tackle the same type of calculation.

I think you are trying to do something that is beyond the capabilities of the LAMMPS variable mechanism. The data from compute vcm/chunk is a local array and there is no way to process that unless you access individual elements like in c_1[1][1], c_1[1][2] and c_1[1][3] for x, y, and z of the first chunk. In addition, I suspect you are skipping over the unit conversions.

To do more sophisticated data processing, you need to use the embedded python interpreter which can then access data from fixes and computes for all cases.

here is an example for the data.spce system with 1500 SPC/E water molecules from examples/PACKAGES/tally with some comments added to explain what what does.


# standard system setup
units           real
atom_style      full

read_data       data.spce

pair_style      lj/cut/coul/long 12.0 12.0
kspace_style    pppm 1.0e-4

pair_coeff      1 1 0.15535 3.166
pair_coeff      * 2 0.0000 0.0000

bond_style      harmonic
angle_style     harmonic
dihedral_style  none
improper_style  none

bond_coeff      1 1000.00 1.000
angle_coeff     1 100.0 109.47

special_bonds   lj/coul 0.0 0.0 1.0

neighbor        2.0 bin

fix             1 all shake 0.0001 20 0 b 1 a 1
fix             2 all nvt temp 300.0 300.0 100.0

# declare python code as inline text.
# this gets the pointer to the current LAMMPS class passed as an argument so it can access the data.

python source here """

def get_avg_ke_mol(lmpptr):
   # initialize python lammps class from C++ pointer
   lmp = lammps(ptr=lmpptr)

   # conversion factor for (mass)*(velocity)*(velocity) to energy
   mvv2e = lmp.extract_global('mvv2e')
   # compute chunk/atom returns the number of chunks as global scalar
   numchunk = lmp.extract_compute('molchunk', LMP_STYLE_GLOBAL, LMP_TYPE_SCALAR)
   # compute vcm/chunk returns the center of mass velocities as global array[numchunk][3]
   molvcm = lmp.extract_compute('molvcm', LMP_STYLE_GLOBAL, LMP_TYPE_ARRAY)
   # we need the mass per chunk to compute the energy
   molmass = lmp.extract_compute('molmass', LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR)

   # for computing the average
   avg_ke = 0.0
   # loop over chunks
   for i in range(0,int(numchunk)):
     # kinetic energy of one chunk: 0.5*m*v^2 in 3-d.
     mol_ke =  0.5*molmass[i]*mvv2e*(molvcm[i][1]*molvcm[i][1]+molvcm[i][2]*molvcm[i][2]+molvcm[i][3]*molvcm[i][3])
     # print("mol[%d] ke = %g" % (i, mol_ke))
     avg_ke += mol_ke
   avg_ke /= numchunk
   #print('avg ke: ', avg_ke)
   return avg_ke

# define chunks and compute vcm/chunk
compute molchunk all chunk/atom molecule
compute molvcm   all vcm/chunk molchunk

# we need the mass per chunk. compute from per-atom mass with compute reduce/chunk
compute mass     all property/atom mass
compute molmass  all reduce/chunk molchunk sum c_mass

# define python style variable to get access to the return value of the python function
variable avg_ke python get_avg_ke_mol
# declare interface to python function: we have one input (the pointer to the LAMMPS class ->SELF) and one return value thus "format" must describe to types.
python get_avg_ke_mol input 1 SELF return v_avg_ke format pf

# output computed average as returned by the python variable (which invokes the python function)
fix p all print 10 "avg ke/mol: ${avg_ke}"

timestep        2.0
thermo          10

run 50

Working with per-atom variables and computes can be more reliable. You can spread per-chunk quantities into each atom in the chunk with compute spread/chunk/atom, and combine per-atom quantities back into per-chunk quantities with compute reduce/chunk.

These transformations open up many possibilities, and one little additional trick is to designate a group of “molecular representatives”: use reduce/chunk to find each molecule’s unique atom with min (or max!) ID, spread back to all atoms, and then the group of molecular representatives just contains all atoms whose mol-min-ID is the same as their own ID. Then the “count” of that group is automatically the count of all molecules in the system, and spreading molecular quantities back into that group and summing or averaging over only that group automatically gives you sums or averages over molecules instead of atoms.

Here’s a generic script snippet illustrating the key concepts. I’ve implemented the total dipole (for which there’s an easy no-molecules-involved alternative so that you can verify it for yourself), and implemented something like what you might be trying to do in your code snippet. You’ll have to further customize, test, and debug it for yourself, in particular if you have mono-atom “molecules” in the mix.

# molecular chunk setup
compute molchunk all chunk/atom molecule
variable id atom id
compute minid    all reduce/chunk molchunk min v_id
compute dipole   all dipole/chunk molchunk
compute vcmmol   all vcm/chunk molchunk
compute spreads  all chunk/spread/atom molchunk c_minid c_dipole[1] c_vcmmol[1]
variable isrep atom (id==c_spreads[1])
run 0 # otherwise variable isrep isn't ready for next line
group molrep variable isrep
print "number of molecular representatives = $(count(molrep))"

# get total x-dipole in two different ways:
# from summing dipole/chunk over molecules
# or directly summing over atoms
# compute totdip1 only includes group "molrep"
# if group "all" is used, overcounting occurs
compute totdip1  molrep reduce sum c_spreads[2]
compute xu all property/atom xu # need unwrapped x in PBC
variable dipatom atom q*c_xu
compute totdip2  all    reduce sum v_dipatom
# "thermo ... c_totdip1 c_totdip2" will confirm they give same result

# reduce-sum atomic vcm-x-kinetic energy into molecular totals
variable kex atom 0.5*mass*c_spreads[3]^2
compute mol_kex all reduce/chunk molchunk sum v_kex
# spread back to mol-representatives
compute spread_mol_kex all chunk/spread/atom molchunk c_mol_kex
# and average:
compute ave_mol_kex molrep reduce ave c_spread_mol_kex

By the way:

When I run this on my LAMMPS (2Aug2023) it errors out with “Atom vector in equal-style formula”, which it should, since mass is a per-atom quantity which can’t be interpreted in an equal-style formula. So I’m not sure what exactly happened on your machine.

Thank you for the reply, Axel. Normally I just do all calculations in post processing using python. I was attempting this calculation within LAMMPS to get a better understanding of how some of LAMMPS’ built-in functionality works, so the responses have been helpful. Am I understanding correctly that compute vcm/chunk outputs a global array where the columns are x, y, and z and the rows are the chunks? Reading the Sec 8.3.1 - Scalar/vector/array data, I was not positive about the array indexing, but I took it to mean that if I took my compute c_1 and indexed it with c_1[1], this essentially cast the first column as a vector. Did I misunderstand what referencing c_1[1] is doing? If it is treating c_1[1] as a vector of length N_chunk, I would think I could do calculations on it the same as the output of compute gyration/chunk, which is also a vector of length N_chunk. I was able to average c_rg with an equal-style variable.

In these examples, I’m essentially just leaving everything in LJ units and setting most of the physical parameters to one, since it was just a toy system to see how some of the computes and variable styles work. So yes, I wasn’t really paying attention to unit conversions as all I was interested was whether the computes and variables were outputting the results I expected based on my understanding of the documentation.

Thanks for the guidance.


@srtee Thank you for the responses. That method casting per-chunk quantities to per-atom quantities and using “representative” atoms for the molecule is a really clever way to access some extra functionality. In fact, in my toy example, I didn’t even multiply by the mass or 0.5 as I wanted something simple to compare to an external calculation, so in reality my variable was just

variable vx2 equal ave(c_1[1]^2)

I probably should have just left them off my example to avoid confusion. Good catch on the mass though.


I was searching for a way to count the number of molecules in a sample, and ended up here. I hope you don’t mind reviving this post. Brilliant suggestion by @srtee!
In case someone else needs it, here is another way to count the molecules:

variable id atom mol
compute mxid all reduce max v_id
variable nmol equal c_mxid
run 0 post no
variable nmol equal $(v_nmol) # store the value before deleting the compute.
uncompute mxid
print "number of molecules = ${nmol}"

One thing to watch out for, the molecule ids can start at a number other than one, so the max molecule ID isn’t guaranteed to be equal to the count.

1 Like

For a gap-tolerant version:

compute molchunk all chunk/atom molecule nchunk once ids once compress yes
variable molnum equal $(c_molchunk) # for immediate evaluation to static integer

Note you may need to run 0, use thermo_style custom to ensure the chunk compute is invoked, etc.

1 Like

There is always something to learn from you guys! Thanks for pointing out to the fallacy in my approach.