Adding a new forcefield

[Simone Sturniolo]

Hi all,

I'm trying to make my own fork of i-Pi to test an algorithm that I'm working on. This is supposed to be an alternative to PIMD, working with a similar concept (multiple classical beads coupled through a fictitious force that reproduces the quantum properties). So my idea was to add an implementation of it as an alternative and then run simulations with it using DFTB+ as an engine for interatomic forces. Unfortunately, I'm having trouble finding where are the PIMD spring forces effectively computed, and so I can't figure out how to edit or reproduce that part of the code. Could someone guide me to the right module? Thanks!

Cheers,

Simone

[venkat kapil]

Hi Simone, the spring forces are computed in ipi/engine/beads.py but in my understanding they are not used to integrate the equations of motion w.r.t the spring term of the Hamiltonian since it cleaner and more efficient to do it analytically in terms of normal modes of the free ring polymer (ipi/engine/normalmodes.py).  Perhaps Michele can correct me if I am wrong. 

[Simone Sturniolo]

I see - is there any other example that is straightforward enough to transform by replacing the forces instead? Maybe the NEB version of the forces?

[venkat kapil]

You can look at the implementation of Suzuki-Chin (SC) PIMD (<dynamics mode= 'sc'>)  where an extra term that depends on the forces used to integrate. This extra force is a depend object and is computed in forces as fsc. So what you can do is that define your extra force in forces as f_whaever and then use it to update your momenta and positions by creating pwhateverstep and qstep. You will have to disable the normal mode propagation though. Have a look at the regular NVT Integrator and the SC Integrator. Should become more clear.  

[Michele Ceriotti]

If you can cast the problem in the normal modes basis, it is much better as it'll probably be possible to integrate using a much larger time step. If not this is a bit tricky because the ForceField object is designed to get one bead at a time (or many, but each independent of each other). Otherwise, in my opinion the simplest path to get a working and relatively self-contained implementation would be to implement a new Dynamics object, that computes whatever fancy bead-bead forces you need as a subroutine of the integrator module.

The only downside is that it'll be hard to "propagate" the information on these forces to the output module, but I don't think this is a problem for a proof of principle.

Cheers

Michele

[Simone Sturniolo]

I don't think normal basis mode is really an option, as the forces are largely anharmonic. So I guess this could be my way to go. Thanks! I'll take a look and write again if I need more details.

Simone