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.
HTH,
Axel.
# 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 """
from lammps import lammps, LMP_STYLE_GLOBAL, LMP_TYPE_SCALAR, LMP_TYPE_ARRAY, LMP_TYPE_VECTOR
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
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.
@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}"
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.