Question on implementing a new sampling method

[vineetbansal]

Hello,
Sorry if this is not the correct place for this question, but I'm working on implementing Vendi Sampling on LAMMPS (and hopefully GROMACS too), from a prototype Python implementation.

The basic idea is that in a multi-replica simulation, the atomic positions in the replicas are inspected, evaluated for "diversity" (w.r.t centered/rotated reference positions), and go through an algorithm that gives us a scalar value called the "Vendi Score", the collective potential energy of the ensemble.

I'm trying to understand how to approach this from a development perspective. My understanding so far on reading on LAMMPS and the Colvars module is, that a suitable plan for development might be:

• Implement a new analysis method (say vendi_sampling), a class that is a child of colvarbias, much along the lines of metadynamics, with file-based communication between replicas, and something that outputs energy in its calc_energy method. (calc_forces cannot, and need not return anything meaningful since we do not have derivatives of energy w.r.t. the atomic positions).

The resulting colvar configuration might look something like this:

colvarsTrajFrequency 5

colvar {
  name atom_positions

  cartesian {
    name xyz
    atoms {
      name ala
      atomNumbersRange 1-22
      centerToReference on
      rotateToReference on
      refPositionsFile alanine-dipeptide-implicit.xyz
    }
  }

}

vendi_sampling {
  name vendi
  colvars atom_positions

  multipleReplicas on
  replicaUpdateFrequency 1
  replicasRegistry vendi-replicas.txt
  outputEnergy on
}

• Use this colvar configuration as a fix in LAMMPS and use the fix_modify command to run the simulation. The documentation on fix_modify states:

The fix_modify energy option is supported by this fix to add the energy change from the biasing force added by Colvars to the global potential energy of the system as part of thermodynamic output

So this might look something like:

  fix myfix all colvars my.colvars
  fix_modify myfix energy yes
  minimize 0.001 0.001 1000 1000

I'm a software engineer with no background in MD before this, so I'm wondering if someone can review this general plan and see if makes sense? Or perhaps there's a simpler way?

[jhenin]

Hi @vineetbansal
Converting this Issue to a Discussion, as I think this is a more suitable, open-ended forum.

Your general approach sounds correct. From scanning the paper, I see a Vendi force in the atomic equation of motion of eq. 12. This looks like it should be a derivative wrt. atomic positions.

[vineetbansal]

Thanks @jhenin - in the Python implementation of the approach that I'm looking to port over, we're using the TorchForce module, something that can autodiff the energy calculations:

you create a PyTorch model that takes a (nparticles,3) tensor of particle positions (in nanometers) as input and produces energy (in kJ/mol) or the value of the collective variable as output.
The TorchForce provided by this plugin can then use the model to compute energy contributions or apply forces to particles during a simulation.

In the paper, eq (10) gives us the energy, but my understanding is that the calculation of Vendi Force used in eq (12) relies on an autodiff and cannot be computed analytically.

So I guess my question is that can I just supply calc_energy and ignore calc_forces since there is no autodiff mechanism available in Colvars, and rely on LAMMPS for the correct time-evolution of the system?

[HanatoK]

No. I think you have to provide the biasing forces to make the dynamics correct. You can add libtorch following the same idea in #570 and use it in Colvars via the torchscript C++ API.