Re: [ESPResSo-users] Individual particle energies

[Rudolf Weeber]

Hi Rhys,
> Thanks, This method does work but I'm getting noisy results so will need to 
> take more samples than anticipated which could be quite slow.
> 
> How complicated would it be (for me?) to try implement the feature for 
> non-bonded interactions only? I'm not super familiar with the the back end of 
> espresso but would it be possible to access the verlet list for the particle 
> and then sum using the ```calc_non_bonded_pair_energy()``` from 
> ```core/energy_inline.hpp``` function? I think that I would run into issues 
> with whether or not the verlet lists are suitable when doing widom tests etc, 
> are the verlet lists recompiled every time you add/move particles manually?
If it is only about non-bonded and bonded interactions, but not about 
electrostatics, the solution would roughly look like this:


Take src/core/energy.cpp:energy_calc() as a starting point
Modify the kernels in the short_range_loop() call in that function
short_range_loop([](Particle &p) { add_single_particle_energy(&p); },
                   [](Particle &p1, Particle &p2, Distance &d) {
                     add_non_bonded_pair_energy(&p1, &p2, d.vec21.data(),
                                                sqrt(d.dist2), d.dist2);
                   });
In case you are unfamiliar with the syntax, please checkout C++ lambda 
functions.
The short_range_loop() function takes care of the cell system and verlet list 
details and just executes the singe particle kernel on all particles and the 
pair kernel on al pairs within the Verlet list or neighboring cells.

You can just include an if-statement which skips the pair kernel unless one of 
the two particles is the one you  are interested in. Particle ids are stored in 
p.p.identity
when p is the particle.

Bonds are more complicated, because they are only stored on one of the 
particles being bound.
This means that, if you need to take into account bonds, you would need to 
write a custom
energy_inline.hpp:add_bonded_energy() which skips bonds if any of the partners 
has the id you are interested in.

On a single core, you would then expose your 
double single_particle_non_bonded_energy(int id) 
function to python by means of a
cpef extern e.g., in src/python/espressomd/analysis.pyx.

On more than one MPI node, you would additionally need code to call the 
function on the other nodes and an mpi_reduce for the energy result.
If you decide to pursue the project, I can provide more details.

Hope that helps
Regards, Rudolf

--
Dr. Rudolf Weeber
Institute for Computational Physics
Universität Stuttgart
Allmandring 3
70569 Stuttgart
Germany
Phone: +49(0)711/685-67717
Email: address@hidden
http://www.icp.uni-stuttgart.de/~icp/Rudolf_Weeber