Non standard Thermodynamic inplementation via compute_pair_local class

Dear Developers,
I am working on a problem where we would like to compute the free energy to create cavities in materials/complex fluids.
A simple approach would seem to consist of adding psuedo (Lennard Jones atoms with a cut-off of negligible initial diameter and slowly scale them up during the simulation, where they are included in the force and potential evaluations, but are not part of the time step integrating group. Note that the ti fix in lammp will not do the job, because it assumes a linear dependence in the perturbation parameter (lambda) (which is sigma in my case)

So far, making holes is easy, BUT my difficulty is extracting the derivative of the corresponding potential with respect to the particle diameter from lammps - lets call it dvds (step 4 below). Mathematically the problem is simple, but implementing it in LAMMPS seems tricky to
say the least.
My idea for extracting the derivative was to:
[1] modify ComputePairLocal::compute_pairs
to compute dvds in the similar way that it calculates the energy and pair force through the pair_lj_cut class (and declare the compute_pair_local as a friend of the pair_lj_cut class so its private members can be accessed)
[2] add an additional array *dvbuf to the buffers - double *dbuf,*ebuf,*fbuf; to accumulate the result
[3] only accumulate to dvbuf atom for pairs where one atom is a pseduo atom
[4] extract the relevant data in the pair_cut_lj class - in particular the private members lj1[i][j], lj2[i][j], lj2[i][j], lj4[i][j] and similar quantities
If I could access these members I could with a few lines compute the terms to be added to dvbuf

I have compiled and run steps [1] to [3] - all seems fine.
BUT step [4] is a head-banger. I believe I am confused regarding access of data from a private member of a class to a member of a different class.

Any suggestions ?
Thanks & Best Regards,
Donal
Dr Donal MacKernan
University College Dublin

Dear Developers,

dear donal,

I am working on a problem where we would like to compute the free energy to
create cavities in materials/complex fluids.
A simple approach would seem to consist of adding psuedo (Lennard Jones
atoms with a cut-off of negligible initial diameter and slowly scale them
up during the simulation, where they are included in the force and potential
evaluations, but are not part of the time step integrating group. Note
that the ti fix in lammp will not do the job, because it assumes a linear
dependence in the perturbation parameter (lambda) (which is sigma in my
case)

So far, making holes is easy, BUT my difficulty is extracting the
derivative of the corresponding potential with respect to the particle
diameter from lammps - lets call it dvds (step 4 below). Mathematically the
problem is simple, but implementing it in LAMMPS seems tricky to
say the least.
My idea for extracting the derivative was to:
[1] modify ComputePairLocal::compute_pairs
to compute dvds in the similar way that it calculates the energy and pair
force through the pair_lj_cut class (and declare the compute_pair_local as
a friend of the pair_lj_cut class so its private members can be accessed)
[2] add an additional array *dvbuf to the buffers - double
*dbuf,*ebuf,*fbuf; to accumulate the result
[3] only accumulate to dvbuf atom for pairs where one atom is a pseduo atom
[4] extract the relevant data in the pair_cut_lj class - in particular the
private members lj1[i][j], lj2[i][j], lj2[i][j], lj4[i][j] and similar
quantities
If I could access these members I could with a few lines compute the terms
to be added to dvbuf

you are making your life needlessly complicated.
why you a lennard jones potential, when a harmonic
one would do just as well.

second, what you describe is what is generally referred to
as a "collective variable" and there are several approaches
to work with them and use them to compute (free) energies.

there are currently two frameworks for this, that are available
for LAMMPS, there is the PLUMED package
( http://www.plumed-code.org/ ) and there is the COLVARS
module (originally implemented for NAMD) that i recently
integrated into LAMMPS-ICMS and soon hope to have merged
(together with other bugfixes, new features and enhancements)
into the mainline LAMMPS code.

if you don't want to mess with that, but rather roll your own,
i suggest you consider the following strategy.

1) implement a fix (check out fix indent as a template) that
  produces a spherical potential repulsive only harmonical
  potential (that'll make computing derivatives easy).
  doing lennard jones is technically possible, but what
  would be the benefit?
2) compute the repulsive force of the pseudoparticle on
  each atom and apply it to the atoms.
3) sum it up, post-process any which way you like.
4) use the compute_scalar(), compute_vector(), or
   compute_array() infrastructure to make your results
   available for writing to files (e.g. via print or thermo output)
5) rejoice.

axel.

You could also talk to Sai, who wrote the TI features
in LAMMPS and knows a lot about various ways
to compute free energies. He is CCd.

Steve