Hello all,
I am trying to write a compute (or series of computes), but I’m not sure if what I want to do is possible inside lammps. I have two types of atoms in my system, say type 1 and type 2. For each atom in type 2, I want to know what is the aggregate force contributed by all of the atoms in type 1. So basically, if there are n atoms in type 2, I am looking for n forces. (Let’s assume all forces are pairwise.)
I know I can use “compute pair/local” to get the pairwise forces of all atoms on all other atoms, but is it possible to use “compute reduce” or something similar to simplify all that data for me? Or should I just do the post-processing separately? Thank you.
Jonathan
jonathan,
Hello all,
I am trying to write a compute (or series of computes), but I'm not sure
if what I want to do is possible inside lammps. I have two types of atoms
in my system, say type 1 and type 2. For each atom in type 2, I want to
know what is the aggregate force contributed by all of the atoms in type 1.
So basically, if there are n atoms in type 2, I am looking for n forces.
(Let's assume all forces are pairwise.)
I know I can use "compute pair/local" to get the pairwise forces of all
atoms on all other atoms, but is it possible to use "compute reduce" or
something similar to simplify all that data for me? Or should I just do the
post-processing separately? Thank you.
many MD codes have an option to 'rerun' a precomputed trajectory for this
purpose. typically you would use this to compute some averaged properties
and then you don't need every time step, as those are highly correlated and
thus contribute little to the overall statistics.
thus instead of a compute, you would then implement a fix that takes
as arguments
a list of file names in your preferred trajectory format, then this
would be used
instead of fix nve/nvt/npt/... to propagate your system and you would
set all force
field parameters that you don't care about to zero, or don't define them at all.
the fix would then read and distribute the coordinates frame by frame, you
would compute the forces (and energies and ...) and then dump them into
a new trajectory file.
BTW: this approach could also be used to approximate analysis methods
that require pairwise additive potentials for interactions that are
not pair-wise
additive, provided that is a reasonable approximate potential. this way the
trajectory would be "correct" and only the analysis would be approximate.
cheers,
axel.
If you look at compute group/group or compute pair/local,
you will see that they use a neighbor list and they double loop
over pairs of atoms and call pair->single() to get a force
or energy. Compute group/group accumulates this
into one force; compute pair/local leaves it as N*M forces
where M = neighbors/atom. You are asking for something
in between, to get N forces. So you could use either
of those as templates and sum up the results as you wish.
Or you could write a reduce-like compute to take the output of compute
pair/local and sum N*M values into N.
Or you could post-process the output of pair/local written
to a dump file via dump/local.
Steve
Hi Axel and Steve,
Thanks for your suggestions. I will consider what’s the best way to proceed.
If I am using pair_style hybrid with lj/cut and airebo, will I be able to pick out the pair->single() function from the lj/cut specifically? I tried doing the pair/local blindly, and it returned with an error. Thanks.
Jonathan
If you just call the pair hybrid's single() function it will
pick the right one for the sub-style, e.g. lj/cut. I think
all you need to do is turn compute group/group into
compute group/group/atom. Which could be done
by keeping all the group/group code and simply
tallying up the returned energy/force values into
a per-atom array instead of into a single scalar.
Steve