Dear LAMMPS,
I am writing a pair of computes which rely on the spherical harmonic functions to calculate parameters of local and global order, and would like to prepare them for submission to LAMMPS as either a standard compute or user package.
Presently, I use BOOST to evaluate these functions, but this introduces a large external dependency to my computes. I could also use recursion relations to generate them in place, or code them directly (I will only need l <= 10), but this would not be optimal performance wise. How do you recommend I proceed?
Thanks,
Anthony Frachioni
Dear LAMMPS,
I am writing a pair of computes which rely on the spherical harmonic
functions to calculate parameters of local and global order, and would like
to prepare them for submission to LAMMPS as either a standard compute or
user package.
Presently, I use BOOST to evaluate these functions, but this introduces
a large external dependency to my computes. I could also use recursion
please check, if this is a "header-only" feature of BOOST. most of
BOOST is supposed to be like that and thus rather than having LAMMPS
depend on an external package, you could just take a tested copy of
the relevant header files and write a small interface library to it
and place it into the lib directory. this is for example already done
with the BLAS/LAPACK support for the AtC and AWPMD packages. i've also
done the same for the voro++ library in LAMMPS-ICMS, but that has not
yet migrated into the upstream version.
this way you can have your cake and eat it.
relations to generate them in place, or code them directly (I will only need
l <= 10), but this would not be optimal performance wise. How do you
recommend I proceed?
how expensive is your analysis in the first place? and how much of
that is due to evaluating spherical harmonics (which tends to be
expensive, i give you that)? there is little use worrying about
optimizing something that is not taking up a significant amount of CPU
time.
axel.
First off, it sounds like you are coding up the Steinhardt
bond-orientational order parameter Q_l and W_l. These are very useful
and it would be great to have them in LAMMPS as a compute. Please
send it to us when it is ready.
The SNAP and GAP potentials that I am working with use Wigner U
functions, which are analogs of the spherical harmonics on the
3-sphere. For given arguments (theta0, theta, phi) we need the values
of all of the functions up to some maximum order. Hence it is much
more efficient to code up the recursion relations rather than use lots
of calls to a library (if one existed). Even if this is not true in
your case, there are probably specific aspects of your calculation
that you can use to eliminate duplicate computation. I am guessing it
would not be hard for you to write your own function for specific
values of l'm that is faster than the BOOST version.
But if BOOST is really faster, then it would be nice to provide a
platform-configurable interface in the LAMMPS lib directory as Axel
mentioned, but also provide the option to compile with built in code,
because *many users and platforms* will not have access to BOOST. Of
course, if it is possible and legal to put the relevant BOOST source
files in lib, that would be ideal.
Aidan
[...]
that you can use to eliminate duplicate computation. I am guessing it
would not be hard for you to write your own function for specific
values of l'm that is faster than the BOOST version.
the way boost works is that it expands the corresponding functions
from templates in the compiler at compile time. most fast
implementations that i know either have a small generator tools that
is compile during compile time and will then generate all spherical
harmonics (and derivatives, if needed) up to a maximum l, or people
have generated that code with mathematica, maple or similar for all
supported l and m and then use if statements to pick the specific
function required.
But if BOOST is really faster, then it would be nice to provide a
i would expect the boost version to be very fast, since it basically
uses the compiler as source code generation tool.
platform-configurable interface in the LAMMPS lib directory as Axel
mentioned, but also provide the option to compile with built in code,
because *many users and platforms* will not have access to BOOST. Of
course, if it is possible and legal to put the relevant BOOST source
files in lib, that would be ideal.
the way to do this would be to have a generic API for spherical
harmonics (or "special math" functions in general) and then copy only
the absolute minimum number of headers. boost can easily lead us into
"dependency hell", so it is essential to "quarantine" it. we
definitely do *not* want to bundle all of boost.
the boost license is BSD-like and GPL compatible, so we can include it.
http://en.wikipedia.org/wiki/Boost_Software_License
axel.
Axel,
I see. Sound like this is your chance to create ../lib/boost. I say
"Bring it!" :-).
Aidan
Axel,
I see. Sound like this is your chance to create ../lib/boost. I say
"Bring it!" :-).
am currently (for several weeks, actually) working on a lib/qmmm
(based on work from people here in trieste), so lib/boost will have to
wait its turn. 
axel.
I still think that calling the BOOST spherical harmonic function will
be slower than application-specific code. For example, if you know you
need Y^l_m(theta, phi) for all values of m, then BOOST will have to
perform 2l+1 calculations of exp(i*m*phi), whereas custom code
requires just a single call to exp(), one division and 2(l-1)
multiplies:
z1 = exp(i phi)
z1p = 1.0/z1
for m = 2,l {
zm *= z1
zmp *= z1p
}
This is similar to what is done in Ewald::eik_dot_r().
Aidan
...and there is another factor to consider: code generated with meta
programming for recursive algorithms can become quite bloated. while
each individual code path may be fast, it creates significant pressure
on the L1 cache for instructions. if the code doesn't fit into it, it
will become much slower. compact hand-written and tuned code, hits
this threshold much later.
so in the end it boils down to what the actual implementation needs.
...and as i mentioned in my initial response, it also depends on how
much time of the total time is spend on evaluating the spherical
harmonics overall and how much time it takes relative to the overall
time per timestep. we may be chasing a while goose here...
axel.
Dear Lammps users/developers,
I am slightly confused as to why the graph of the potential energy (rough graph sketch attached) of a Kremer-Grest polymer chain during relaxation after deformation fluctuates so much...
Initially, I scale the box 2 times in x direction, keeping y and z constant. After this, I relax the chain and keeping the box size extended in x, let the chain relax to obtain the new equilibrium. I use fixed boundary conditions for both deformation and relaxation, and nve+langevin for both. During relaxation, I also impose a soft wall potential. Is this due to wall interactions perhaps... Also, I cannot really explain the dip around ts 10000... Any suggestions are welcome, thanks a lot in advanced.
Anna
Dear Lammps users/developers,
I am slightly confused as to why the graph of the potential energy (rough
graph sketch attached) of a Kremer-Grest polymer chain during relaxation
after deformation fluctuates so much...
Initially, I scale the box 2 times in x direction, keeping y and z
constant. After this, I relax the chain and keeping the box size extended
in x, let the chain relax to obtain the new equilibrium. I use fixed
boundary conditions for both deformation and relaxation, and nve+langevin
for both. During relaxation, I also impose a soft wall potential. Is this
due to wall interactions perhaps... Also, I cannot really explain the dip
around ts 10000... Any suggestions are welcome, thanks a lot in advanced.
Your relaxation procedure sounds quite convoluted. Are you trying to relax
one polymer chain? If yes, why do you need to strain the box? And there is
a wall in the system as well? I don't do polymers but your whole
explanation is very confusing. Perhaps you should start by better
describing your system and goals. Next, a minimal input deck would help...
Carlos
Thank you Carlos, I agree, I did not give a clear explanation, will do so soon and provide an input file!
With best wishes,
Anna