8

I'm conducting a molecular dynamics simulation for silica. Some time ago I turned to the fluctuating dipole model, and after much effort I'm still having problems implementing it.

In short, all oxygen atoms in the system are polarizable, and their dipole moments depend on their position with respect to all the other atoms in the system. More particularly, I use TS potential (http://digitallibrary.sissa.it/bitstream/handle/1963/2874/tangney.pdf?sequence=2), where dipoles are found iteratively at each time step.

This means that when evaluating forces acting on atoms, I have to take into account this potential energy dependency on coordinates. Before, I was using simple pairwise potential models, and so I would set my program to compute forces using analytic formulas obtained by differentiating potential energy expression.

Now I'm at a loss: how to implement this new potential? In all the articles that I've found they only give you formulas, but not the algorithm. As I see it, when I compute forces, acting on a certain atom, I have to take into account the change of the dipole of this atom, the change of dipoles of all the neighboring atoms, then tge change of dipoles of still more atoms, and so on, as they depend on each other. After all, it is because of this interdependency that the dipoles are found iteratively at each time step. Clearly, I can't compute forces iteratively for each atom, because the computational complexity of the algorithm would be way too high. Should I use some simple functions to account for the change of dipoles? This doesn't look like a good idea either, cause dipoles are calculated iteratively, with high precision, and then, where it actually matters (computing forces), we would use crude functions?

So how do I implement this model? Also, is it possible to compute forces analytically, as I did before, or is it necessary to compute them using finite difference formula for derivative?

I haven't found the answer for my question in literature, but if you know some article, or site, or book, where this material is highlighted, please, direct me to that source.

Thank you for your time!

==================================================================================

UPDATE:

Thank you for your answer. Unfortunately, this was not my question. I didn't ask how to compute dipoles, but how to compute forces given that those dipoles vary with movement considerably.

I tried to compute forces in straightforward manner (not taking into account dipoles interdependance via their distances, just compute the dipoles on each steps, then compute the forces as if those dipoles are static), but the results I got were not physically correct.

To analyze the situation, I set up a simulation of a system consisting of just two atoms: Si and O. They have opposite charges, and so they oscillate. And the energy time dependence graphic looks like this:

enter image description here

The curve on the top represents kinetic energy, the one in the middle represents potential energy without taking dipole interaction into account, and the one at the bottom represents potential energy of the system, where dipole interaction was taken into account.

You can clearly see from this graphic, that the system is doing what it shouldn't do: climbing up the potential slope. So I decided this is due to the fact that I didn't take dipole moment coordinate dependence into account. For instance, at a given time point, we compute the forces, and they are directed so as to move both atoms toward each other. But when we do move them towards each other (even slightly), the dipole moment changes, and we find out that we actually ended up with higher potential energy than before! During the next time step the situation is the same.

So the question is, how to take this effect into account, cause what few ways I can think of are either way too computationally intense, or way too crude.

Denis
  • 163
  • 1
  • 8
  • 1
    * mind blown. To be honest, I'm not sure if you'll get an answer here (although I could be wrong). You may have better luck at physics.stackexchange.com, since you may find someone there with a knack at computational physics. – tay10r May 22 '13 at 07:53
  • Thank you for the link. I guess, I'll try asking there as well, though I'm not impressed by the sort of questions I found on the main page there) – Denis May 22 '13 at 07:59
  • I agree, they can be a bit general. In any case, I wouldn't mind setting a bounty on the question if you don't get an answer. – tay10r May 22 '13 at 08:02
  • 2
    This is probably a better fit for your question: [Computational Science SE](http://scicomp.stackexchange.com/). – Dan May 24 '13 at 07:53
  • Have you checked out `NAMD` http://www.ks.uiuc.edu/Research/namd/features.html? The source is available so you may be able to leverage some of their methods. – Shafik Yaghmour May 24 '13 at 16:36
  • @ShafikYaghmour Thank you for the link. Didn't notice your comment at first, sorry. Didn't find any mention of the TS potential on that site, though I might search it again later. A pity it's about biophysics, for the TS potential was designed specifically for silica which I study right now. Maybe I'll write to authors of the potential later... – Denis May 27 '13 at 13:52
  • @TaylorFlores Wow! Thank you for the bounty! I must admit, I'm confused: I've had an impression that users with more or less high score here want to increase it, not to give it away. Besides, I didn't see computational physics in the list of your hobbies in your profile. Thank you anyway. – Denis May 27 '13 at 13:58
  • @Dan Indeed. Asked the same question there as well. Still, this site is more lively, more views, more comments, more users. – Denis May 27 '13 at 14:00
  • 1
    yOU SAY "*In all the articles that I've found they only give you formulas, but not the algorithm*", can you point us to those formula? Then we might have a better chance of figuring out how to algorithmitize them. – RBarryYoung May 29 '13 at 14:39
  • @RBarryYoung Those formulas ar written in section 2 paragraph "B. The force field" of the article I've put the link to at the beginning of the text of my question. They have numbers 2 to 5. – Denis May 30 '13 at 12:53
  • Have you tried changing the size of the time steps? Sometimes, when you simulate a continuous process in discrete steps, the results deviate too much from the reality. In these cases, changing the size of the steps affects the results. Smaller steps lead to more realistic results... until the limited precision of the data, allied with cumulative errors, spoil the whole thing. – comocomocomocomo May 30 '13 at 14:51
  • @comocomocomocomo Yes, I've tried it. It did not help. If you just read my speculations on the essence of the problem, you'll see that even infinitely small time step won't rectify the matter (of course, that is, if my speculations are true). – Denis Jun 06 '13 at 10:49
  • Ok, but... is the result sensible to changes in the step size? (Even though no step size produces the correct result). If yes, then you need better formulas. Maybe the problem doesn't have a solution at all. You can't use anything near infinitely small steps because of limited precision and cumulative error. – comocomocomocomo Jun 06 '13 at 13:51

1 Answers1

1

Not sure I fully understand your question, but it sounds like you might be needing to implement a Markov chain type solution?

See this nice post for more info: http://freakonometrics.hypotheses.org/6803

EDIT. Reason I suggest this is that it sounds like you have a system where state of each atom depends on it's neighbors, and in turn the neighbors state depends on their neighbors and so on. Conceptually this could be modeled as a huge matrix and you iteratively update each value based on it's neighbors (???). This is intractable, but the article linked to shows how to solve a problem of very large transition matrices using Markov chains instead of computing the actual matrix.

Steve
  • 8,469
  • 1
  • 26
  • 37
  • I might know little about Markov chains, but I don't see how they could be fitted here. There's nothing random here. Sorry if that's just my ignorance, can't afford to study unrelated topics just to see if there is in fact some relation right now. – Denis May 23 '13 at 18:39
  • Definitely doesn't sound like MCs could be useful here. (I could, of course, be wrong) – Dan May 24 '13 at 07:54
  • I'll add little more explanation, but this is possibly going off in the wrong direction. – Steve May 24 '13 at 14:03
  • Oh, I see what you're saying. Usually with physical simulations one needs to know the forces exactly. Further, having the force be dependent on the outcome of a stochastic optimization procedure would lead to the method being unconservative of things like energy/momentum. Though this is sometimes desirable for simulating thermal or dissipative systems. – Dan May 24 '13 at 18:19
  • Yes, that is the essence of my problem, but I'm not yet ready to turn to stochastic methods. – Denis May 25 '13 at 07:58