[lammps-users] recompute force after perturbation

Dear all,

I'm attempting to write a compute currently which, among other things,
can recalculate the set of forces (for all coordinates and atoms) upon
a perturbation in atom position.

Certainly the initial positions and forces are available in atom->x
and in the force->pair, force->bond, etc., and therefore it is trivial
for me to access the full set of forces for the current timestep for
any atom and coordinate.

Can anyone think of a way to accomplish such a thing? I have done this
in my own codes, but in that case I have a full description of the
current phase point which is used as an input to a function which
returns all the forces for the system. All I do in that case is rerun
the same function for a slightly perturbed phase space input.

The point of this is a simple finite difference function to compute
curvatures in certain regions of the phase space.

Best,
Anthony

dear antony,

Dear all,

I'm attempting to write a compute currently which, among other things,
can recalculate the set of forces (for all coordinates and atoms) upon
a perturbation in atom position.

Certainly the initial positions and forces are available in atom->x
and in the force->pair, force->bond, etc., and therefore it is trivial
for me to access the full set of forces for the current timestep for
any atom and coordinate.

Can anyone think of a way to accomplish such a thing? I have done this
in my own codes, but in that case I have a full description of the
current phase point which is used as an input to a function which
returns all the forces for the system. All I do in that case is rerun
the same function for a slightly perturbed phase space input.

i would look into the replica package and implement the
perturbation as a (modified) replica of the original. then
all you need to do is to write a fix that would propagate
any replica coordinates according to the master copy,
and then apply the perturbation, compute the forces with
the normal pair styles and then collect the results.

that approach would have the additional advantage to be
embarrassingly parallel and applicable to non-pairwise
additive and long-range potentials.

cheers,
    axel.

All that is available in atom->x and atom->f
is the current position and total force. Nothing
like this is stored in force->pair or force->bond.
I.e. there is no force due to bonds stored anywhere.

The only thing your fix can do is to change atom->x
and re-evaluate all (or some) of the forces. Same cost
as doing an entire timestep, except you may also need
to re-neighbor if the perturbation is large (and then
reneighbor again) when you undo the perturbation.

If you
are just changing one (or a few) atom coords, you
could get the pairwise part of this by invoking
pair->single(). This is what compute group/group does.
there is also a bond->single(), but you'd have to
find what bonds are affected by your perturbation.
There is no angle/dihedral/improper->single().

Steve

Yeah I realized that only after I wrote my email. atom->x and atom->f
is exactly what I need. I like Axel's suggestion of a replica type
method, but it's significantly more complex than simply reevaluating
the forces. If I change atom->x, is there a simple way within a
compute or fix to reevaluate atom->f. The cost is perfectly fine for
me. I will probably just store atom->x and atom->f in 2 new objects,
modify the new atom->x. Then I'd like to recompute the forces, and do
the analysis I need to do, and then move on in the calculation
discarding the modified atom positions and forces.

It's not clear to me that there is a force call that can accomplish
this simply.

Thanks,
Anthony

anthony,

> All that is available in atom->x and atom->f
> is the current position and total force. Nothing
> like this is stored in force->pair or force->bond.
> I.e. there is no force due to bonds stored anywhere.

Yeah I realized that only after I wrote my email. atom->x and atom->f
is exactly what I need. I like Axel's suggestion of a replica type
method, but it's significantly more complex than simply reevaluating
the forces. If I change atom->x, is there a simple way within a
compute or fix to reevaluate atom->f. The cost is perfectly fine for
me. I will probably just store atom->x and atom->f in 2 new objects,
modify the new atom->x. Then I'd like to recompute the forces, and do
the analysis I need to do, and then move on in the calculation
discarding the modified atom positions and forces.

the easiest way to implement this would probably be a modified
verlet run style.

check out the Verlet::run() method in verlet.cpp
that will call the compute() method in all active
pair, bond, angle, etc. styles to compute the forces.

cheers,
    axel.

Axel's is the right answer. There is no single call that
computes all forces, but Verlet shows how to do it in the
line after

// force computations

in the run() routine.

Steve

Axel, Steve,

Thanks for your advice with this. I was able to implement the
calculation without much trouble, and look at everything from little
pieces of the landscape all the way to full hessians. If you think
people would be interested I'd be happy to post it as soon as I make
sure I haven't made too many errors ...

In the case of dumping the information to file, though, is there a way
you can think of to use the dump a global matrix of arbitrary
dimension? It's obviously not local or per-atom global, the only flag
I gave my compute routine was the array_flag.

Best,
Anthony

In the case of dumping the information to file, though, is there a way
you can think of to use the dump a global matrix of arbitrary
dimension? It's obviously not local or per-atom global, the only flag
I gave my compute routine was the array_flag.

What do you mean by a global matrix of arbitrary dimension?
All computes produce specific kinds of outputs, and there
are other commands that can process all of those outputs
and write them to files. See section 4.15 of the manual
for an overview.

To tell you if this is of general use and something to
release in LAMMPS, I'd have to know what it does,
how it's called, and why people would want to
use it. I suggest you write a doc page for it,
using one of the compute txt files in the doc dir
as a starting point.

Steve

The end result of my compute, in the most extreme case, is the full 3N
x 3N Hessian stored in a 3N x 3N double array. The only flag I set in
the compute is "array_flag" and I assigned "array" to my final output
matrix.

I am reasonably certain that my use for this is so specialized that
there would be no utility for anyone else. Nevertheless I will write
something up when I have time.

Best,
Anthony

NxN is potentially huge, and you'd have to sum it
to get the same answer on every proc.

If you have N quantities per atom, you could
make it a per-atom array, exposed by a compute,
but that's still huge. if the "matrix" is sparse,
which I assume it is, then it's really a "local"
object in compute lingo, like compute pair/local
exposes for all pairwise interactions. If you
expose it that way, then it could be output
via dump local, for example.

Steve

Thanks very much for everyones help on this. I'm sure I'll have other
questions in the future but for now things are working quite well!