Time-dependent forces/potentials

[Nicolas Künzel]

Hello,

as I already had a short email conversation with Michele Ceriotti I will put together the already discussed and then have some more questions.

What I want to do for my Masterthesis is to apply an external time-dependent force to a system which pulls the system over an energy barrier in order to calculate the Helmholtz free energy (via Jarzynski's equality).

When using i-PI I therefore get a problem due to the time-dependent force since the force calculation in the clients is just initialized once and then can't be changed anymore while the client is connected. At least not in a useful manner since i-PI doesn't send information about the time-step and the actual MD time to the client because that would violate the encapsulation.

So one idea would be to modify the driver.x in order to get the time-dependent force there but that still leaves the problem with the time sending. Also: if one wants to use this external force in different systems it is not useful to hardcode the external potential and forces into the driver since one considers different reaction paths for different reactions, so the force acts on different collective variables (CVs).
It also would be nice to use already existing packages which can deal with the calculation of forces due to CVs so one doesn't have to implement that by one self.

The idea of Michele Ceriotti now was that I should implement the Jarzynski pull (so the external force) as a new motion class inside i-PI.
But my question here is now: that still wouldn't allow me to use external packages and I would have to hardcode the forces inside there than, right? But maybe I just got the idea wrong.

Another idea I had now was: could one define a new ffsocket (so like ffdebye) which can be disconnected every timestep after its force calculation, then changes the input for the external force calculation and reconnects the client? I have the feeling that this would increase computation time a lot since it always needs to open and write to a file?

Or could I change the socket communication so that it sends the time-step? I know it would violate the encapsulation but maybe it is possible to do that with an if-clause so one has to enable that in the input.xml file?


Thanks for opening this forum and I hope my questions are not written too confusing.

Kind regards

Nicolas

[Michele Ceriotti]

There are many questions in here, and not a very simple answer. Something that I believe could be a solution to all your problems is that we're finalizing an interface to Plumed, that I believe has support for steered MD. For the moment, what can help you tremendously is that probably you want to apply the Jarzinsky force to the centroid of the ring polymer, so you actually will have a single call to the force evaluation per step, so (unless you stop and restart your calculation) you can basically update the timer independent of what i-PI is doing. 

To me, this is the best short-time solution until we get the plumed interface merged, then you'll see how you can pass to Plumed the current step (which is implemented already, although you'll see why that is special for the case of Plumed).

Meanwhile, I think it will make you good to spend some time studying the code and the examples, and perhaps also browse the literature to understand more in detail what it means to compute a free energy in a PIMD context.

Cheers

Michele 

[Nicolas Künzel]

Thanks for the fast answer.

That sounds actually like a nice thing. In what time scale do you think you will finish the interface to Plumed?

Actually I had the same idea of just applying the force to the centroid which should work well with LAMMPS at least. I am already trying to implement that behavior into the driver code for a simple harmonic oscillator and will try that. Should know soon if that is a useful way since the free energy for a harmonic oscillator in the classical and the quantum mechanical case can be calculated analytically.

Cheers
Nicolas

[Nicolas Künzel]

I just wanted to tell you that that procedure actually works. Since a few days now. One reason why it took my quite so long was, that I didn't know that i-PI fixes the centroids of the particles by default so they won't move in the external potential. It took me a long time until I figure that out and set fixcom to false.

Maybe it would be better to set the default value of fixcom to false instead of true? Or just write that clearly in the documentation so one knows that from the beginning?

Cheers
Nicolas

[Michele Ceriotti]

uhm. 99% of the time one wants to fix the center of mass in simulations, however you're right that the documentation needs to be improved.

if you find aspects -- such as this -- that could be clearer, feel free to propose changes to the docs, they're part of the repo.

Cheers

[Mariana Rossi]

A better documentation for the fixcom flag is something that a few other people have asked

and it is in my agenda. We will go through the whole manual soon ;)

[Barak Hirshberg]

Hello Everyone,

I hope it is OK I am jumping in a year later to this discussion.

I found this tutorial online (given by Mariana?)
https://th.fhi-berlin.mpg.de/sitesub/meetings/dft-workshop-2015/uploads/Meeting/i-pi-tutorial-slides.pdf

It says that "the estimators implemented in i-PI are designed for False option. Using True here leads to wrong results for small systems". Then I saw this discussion and I was wondering - have the default values changed since the original post? I could not find a lot of details about fixcom in the documentation. (It still says the default is True, so I am assuming it is correct)

Could you perhaps elaborate on what the expected errors? (Or point me in the write direction on where to read about it?)

Thanks,

Barak

[Mariana Rossi]

Hello Barak,

I have never seen these slides before :-D No clue who was talking about this in 2015 (from the website it seems it was Mateusz Marianski https://th.fhi-berlin.mpg.de/sitesub/meetings/dft-workshop-2015/index.php?n=Meeting.Program.)

In fact it seems to me that that phrase does not sound right, I am not really sure what it means. In i-Pi you can use true or false there and all should be good. There was a bug perhaps in the output of the kinetic energy (only the ouptut, internally all was fine) a very old version, but we fixed it long ago.

All the best,

Mariana

[Mariana Rossi]

P.S: Please use https://th.fhi-berlin.mpg.de/sitesub/meetings/dft-workshop-2017/index.php?n=Meeting.Program if you want to have a look in a tutorial I did last year with ab initio codes and also the tutorials and these: http://ipi-code.org/resources/tutorials/

[venkat kapil]

The dynamics has been fine. The estimator of temperature and the classical kinetic energy were not working with constraints on the center of mass and on the atoms. They have been fixed. 

[Barak Hirshberg]

Thank you Venkat and Mariana for your answers!

And thanks Mariana for the most recent tutorials you suggested! I will definitely use them.

All the best,

Barak