[gmx-developers] Re: Various integrator questions

[Michel Cuendet]

Hi guys,
 
Michael said:

>By the way, there are other, more complicated issues with constraints
>when it comes to pressure control and scaling the system, (which would
>occur with any integration algorithm) as I understand (Michel, you
>could probably explain more!), we should be using the molecular
>virial, as opposed to the atomic virial, and scaling the distance
>between molecular centers, not simply atomic centers.  I'll read up on
>this a little more, but if anybody has looked into this issue more
>carefully, let me know 
>

I think molecular scaling has the advantage that it doesn't interfere 
with intramolecular constraints. The Kalibaeva et al. paper is one of 
those with more equations than text, but I think it leads to a quite 
usable result. The use of molecular virial doesn't seem to add 
complexity with respect to the atomic virial. I haven't really dissected 
all the details, but I see that they are using SHAKE. The point is to 
call it at the right stage of the Trotter decomposition, it seems.

I'm convinced that the molecular virial is the way to go for, say, a box 
of water. But in the case of a huge protein in the middle of a not so 
huge box, I can imagine how it could pose some problems at the boundary 
of the protein. If you think about a homotecy centered at the middle of 
the box, the first water shell would move a lot with respect to the 
protein surface, which is not scaled. I don't know about how the 
compressibility of the proteins should look like. Maybe one could think 
about a "residue virial" since we use constaints usually only inside a 
residue. But this is science fiction, sorry.

David said:

>>>Second, the Nose-Hoover thermostat gives the correct ensemble for your
>>>> > system plus "something else". In fact it therefore does not give you the
>>>> > correct ensemble for what you are interested in.
>>>      
>>>
>>> 
>>> I think Michel addressed this before -- for the physical variables,
>>> you get the correct canonical
>>> distribution.
>>> 
>>    
>>
>But not necessarily the correct fluctuations, as described in the
>GROMACS manual (probably written by Berk, but I'm not sure). This means
>that averages may be correct, and possibly also the amplitude of the
>fluctuations, but the time structure of the fluctuations is artificial.
>Higher moments may also be affected.
>

Well if you have the right distribution, you automatically have the 
right moments, and the right fluctuation amplitude. It is true that the 
NH thermostat shows its virtues only if you give it enough time. But 
this is already far better than the berendsen thermostat, for which no 
formal proof exists about what distribution it generates at all. It is 
true that the time structure of the fluctuations can be affected by 
persistent oscillations of the NH variable (as shown by Holian et al. 
Phys. Rev. E 52 : 2338,1995). But the addition of NH chains solves this 
problem.

Now whether the time structure of the fluctuation is artificial, this is 
rather a philosophical question. On one hand it has to be artificial in 
a way, because the isolated hamiltonian physical system alone would 
evolve at constant energy. On the other hand, how are the energy 
fluctuations induced in nature, when the system is coupled to a "real" 
heat bath (say an infinite solution) ?  Nobody really knows, but all we 
really care about is averages. I think the goal is to have a thermostat 
which is proved to give the right distributions, AND which perturbs the 
dynamics in as little as possible (continuous, homogeneous, and small 
magnitude of velocity scaling). For example in some codes, you can apply 
one NH thermostat per atom: this gives the right canonical distribution 
very quickly, but then the coupling forces outweight the forces 
resulting from the physical potential. This is not what we want, because 
then we are half way between MD and some kind of Monte-Carlo sampling, 
in the sense that the way you get from one point to the other has little 
to do with the physical dynamics.  Some authors come up with some 
"improved" thermostat, and call it better than NH, because it brings the 
system faster to thermal equilibrium. I tend to think that those 
thermostats are not neressarily better for production, if they perturb 
the sytem more than necessary.

Bye,

Michel

-- 
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Michel Cuendet, Ph.D. student
Laboratory of Computational Biochemistry and Chemistry
Swiss Federal Institute of Technology in Lausanne (EPFL)
CH-1015 Lausanne						
Switzerland                            Phone : +41 21 693 0324
lcbcpc21.epfl.ch/~mitch	                       +41  1 633 4194
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++