I have written an implementation of a simple forward-difference
numerical hessian in a LAMMPS compute which returns a vector (I manage
the index addressing manually instead of using an array type). I have
attempted to do this by replacing elements of the atom class with
local versions of the data via pointer gymnastics and a lot of local
data that my compute manages. The forces then run and return into my
own arrays locally, and at the end I put everything back to the
original data.
The problem is, I get different results in the classical dynamics (and
rapidly atoms leave the simulation) when running the hessian compute,
so I must be doing something fundamentally wrong. I have not yet been
able to identify the source of the problem within my own code, so my
guess is that I am overlooking or unaware of something in LAMMPS
itself which is preventing this from working.
In principle, my code should return the LAMMPS atom class with exactly
no modifications, and it's the only LAMMPS memory I mess with.
I attach my full source for completeness. A brief description of the
mechanics of the methods are below.
The forces then run and return into my
own arrays locally,
If you’re re-invoking the forces from a compute, that could be
problematic. Things like pair/bond styles store internal
state that can be used later in the timestep (e.g. energy, virial),
and you have essentially no control over when during the timestep
your compute will be invoked, since various other commands
use fixes & variables that invoke computes at different stages
of the timestep. I can’t think of any other computes that invoke
OK that's disappointing. I figured there were things I had missed
(energies, virials, etc.), but the lack of explicit ordering may kill
my approach entirely, which was to just make a copy of everything
updated during forces and put it all back at the end.
Maybe the thing to do would be to just create a complete copy of the
LAMMPS class, do everything there, and destroy it at the end. Can you
think of anything implemented like this using copies of the LAMMPS
class currently that I could use as a reference? Or is this idea also
problematic for some other reason I haven't considered?
I have done some further reading and realized that there is no copy
constructor available for any of the components or for the LAMMPS
class. I was then able to deduce that the instability in the dynamics
was introduced by the call to comm->reverse_comm(). Without this, I
get normal dynamics for my particular problem. I also removed
forward_comm() with no effect.
For now, my problem is solved and my code basically works. I would
like to understand the situations where not having this communication
would be a major issue. My understanding was that, at the very least,
this would cause problems with newton set to on, however, it doesn't
in my case. In fact, I crash with newton set to off only (an issue to
be debugged later). Any insight would be very much appreciated.
I have done some further reading and realized that there is no copy
constructor available for any of the components or for the LAMMPS
interesting to mention this, i just sent a patch to steve that marks
the automatically generated default and copy constructors for the
LAMMPS class private, so that they are not called by accident, but
rather generate a compile error.
class. I was then able to deduce that the instability in the dynamics
was introduced by the call to comm->reverse_comm(). Without this, I
get normal dynamics for my particular problem. I also removed
forward_comm() with no effect.
For now, my problem is solved and my code basically works. I would
like to understand the situations where not having this communication
would be a major issue. My understanding was that, at the very least,
this would cause problems with newton set to on, however, it doesn't
in my case. In fact, I crash with newton set to off only (an issue to
be debugged later). Any insight would be very much appreciated.
forward_comm() communicates properties from local atoms to its
corresponding ghost atoms and reverse_comm() collects data from ghost
atoms to its original local atoms. have a look at pair_eam.cpp for an
example.
interesting to mention this, i just sent a patch to steve that marks
the automatically generated default and copy constructors for the
LAMMPS class private, so that they are not called by accident, but
rather generate a compile error.
When/if this comes up, I can always pull the default copy constructors
out of private for myself.
forward_comm() communicates properties from local atoms to its
corresponding ghost atoms and reverse_comm() collects data from ghost
atoms to its original local atoms. have a look at pair_eam.cpp for an
example.
Great, this explains my observed behavior. I intentionally throw away
any and all properties computed during the force calls anyway, so I
can ignore it and I shouldn't need to account for whatever it updates.
The forward comm is still important then to me. Thanks.
interesting to mention this, i just sent a patch to steve that marks
the automatically generated default and copy constructors for the
LAMMPS class private, so that they are not called by accident, but
rather generate a compile error.
When/if this comes up, I can always pull the default copy constructors
out of private for myself.
which would be a *very* bad idea, since that would create a shallow
copy and thus using that copied instance would not be what you want.
that is the reason to mark it private in the first place. ignoring
this kind of hint is asking for trouble.
forward_comm() communicates properties from local atoms to its
corresponding ghost atoms and reverse_comm() collects data from ghost
atoms to its original local atoms. have a look at pair_eam.cpp for an
example.
Great, this explains my observed behavior. I intentionally throw away
any and all properties computed during the force calls anyway, so I
can ignore it and I shouldn't need to account for whatever it updates.
The forward comm is still important then to me. Thanks.
i don't share your optimism. like steve said, what you are doing is
rather brutal and ignoring a lot of conventions in LAMMPS. if it works
it is more by chance than by construction.
which would be a *very* bad idea, since that would create a shallow
copy and thus using that copied instance would not be what you want.
that is the reason to mark it private in the first place. ignoring
this kind of hint is asking for trouble.
fair enough. i used the wrong language. i'd write some sort of
non-private deep copy which deals with what i want.
i don't share your optimism. like steve said, what you are doing is
rather brutal and ignoring a lot of conventions in LAMMPS. if it works
it is more by chance than by construction.
i agree it's not pretty. it is nonetheless functioning perfectly. if
you have another suggestion on how to approach implementing it, i'm
all ears.
for that you would first have to describe what it is that you are
doing exactly. or rather what it is that you need LAMMPS to do that is
not directly available. i don't have the time to reconstruct this from
some source code.
I am implementing a general forward difference numerical hessian. I
will drive LAMMPS as a library and extract the result of this
calculation as input to my own codes. I've implemented this many times
in my own research codes but I'm now looking to leverage the broad
capabilities of LAMMPS going forward, instead of continuing to build
MD from scratch. If you are more asking about the motivation, it is
for input to the calculation of a range of tangent space/Lyapunov
vector calculations.
I am implementing a general forward difference numerical hessian. I
will drive LAMMPS as a library and extract the result of this
calculation as input to my own codes. I've implemented this many times
in my own research codes but I'm now looking to leverage the broad
capabilities of LAMMPS going forward, instead of continuing to build
MD from scratch. If you are more asking about the motivation, it is
for input to the calculation of a range of tangent space/Lyapunov
vector calculations.
i meant more on a technical level.
what you describe sounds to me like it would be better implemented as
a command that operates similar to a run or a minimize style. i
definitely would not code it as a compute, as those are purely
"consumers" and not "modifiers". you can internally either use fix
store/state to retain a copy of your initial state or follow a scheme
similar to what the respa run style does (and then you would have
"inner" and "outer" steps and also can decide how to consistently
handle fixes that add/modify forces to the system).
i definitely would not code it as a compute, as those are purely
"consumers" and not "modifiers". you can internally either use fix
store/state to retain a copy of your initial state or follow a scheme
similar to what the respa run style does (and then you would have
"inner" and "outer" steps and also can decide how to consistently
handle fixes that add/modify forces to the system).
this is good info thanks, i will give this approach a whirl.
i guess i initially thought of it as a compute as i strictly want to
modify *nothing* about the classical trajectory or the internal state
of lammps, simply to output a matrix which is constructed via calling
the set of force functions, already set up for me in lammps, ndims *
natoms times. each successive call differs only by a small
perturbation in the system's positions. i need the difference between
the original lammps forces and the new set from the perturbed
positions to construct the hessian. for my use, the hessian is
strictly an "output", or static property of the system at a given
time, even if it temporarily has to do some internal modification to
actually do the calculation.
this is good info thanks, i will give this approach a whirl.
i guess i initially thought of it as a compute as i strictly want to
modify *nothing* about the classical trajectory or the internal state
of lammps, simply to output a matrix which is constructed via calling
the set of force functions, already set up for me in lammps, ndims *
natoms times. each successive call differs only by a small
perturbation in the system's positions. i need the difference between
the original lammps forces and the new set from the perturbed
positions to construct the hessian. for my use, the hessian is
strictly an "output", or static property of the system at a given
time, even if it temporarily has to do some internal modification to
actually do the calculation.
here is another idea worth looking into.
since typically, the force computation is the most time consuming
part, you could also parallelize the process into a multi-replica
method and then set up an integrator that would communicate the
positions after the initial integration step to all replica, then
apply the perturbations and compute the forces in parallel. collect
the results and output it on the original, which has computed the
forces to be used in the next integration.
there should be a few discussions in the mailing list archive with
nicola varini on using this kind of approach to implement an EVB
scheme into LAMMPS. there the differences between replica would be the
bond topology, not coordinates, but the scheme is the same: replica 0
is the "one true" MD run and the other replica are "slaves" that get
the position updates from replica 0 and their own propagation gets
ignored.
the nasty thing about this is setting up the intra and inter replica
communicators, but it is not impossible. it would definitely take the
most advantage of LAMMPS' strengths.
here is another idea worth looking into.
since typically, the force computation is the most time consuming
part, you could also parallelize the process into a multi-replica
method and then set up an integrator that would communicate the
positions after the initial integration step to all replica, then
apply the perturbations and compute the forces in parallel. collect
the results and output it on the original, which has computed the
forces to be used in the next integration.
there should be a few discussions in the mailing list archive with
nicola varini on using this kind of approach to implement an EVB
scheme into LAMMPS. there the differences between replica would be the
bond topology, not coordinates, but the scheme is the same: replica 0
is the "one true" MD run and the other replica are "slaves" that get
the position updates from replica 0 and their own propagation gets
ignored.
the nasty thing about this is setting up the intra and inter replica
communicators, but it is not impossible. it would definitely take the
most advantage of LAMMPS' strengths.
this sounds like the best suited approach to my problem by far, and
even if complicated, in the meantime i have a version that (maybe
miraculously) works. and as i think about it it may be simpler than
the EVB scheme -- each "replica" in my case is trivially parallel and
can be organized to populate its own row or column of the hessian
matrix, so the structure of the post comm isn't as tough i don't think
(though i'll have to look into the EVB implementation more). the only
question mark is whether the overhead is worth it -- even though my
method requires dn force calls, it is trivial in comparison to the
work done by my codes (which scale as ~dn^3). it may turn out that
parallelizing the forces does me no good with the limited system sizes
i'm able to achieve with my post analysis. still though, if it's
cleaner and allows for future possible scaling, it's worth
investigating.
If you are calculalting the Hessian thru the library
interface as an analysis tool, then you have control
over how the compute is used. As opposed
to just writing a compute and expecting
it to work in any manner that users might
invoke it. As I said before the “problem”
with a compute that does what yours
is doing (I think, don’t fully understand it),
is that there are many ways to use and invoke
a compute in LAMMPS, and I’m guessing in
some of those modes, what you are doing
wouldn’t work.
But if you create a tool that wraps LAMMPS thru
its lib interface and thus invokes the compute
exactly the way you need it to be invoked,
then it’s probably fine.