[gmx-users] using one thermostat for entire system ??

[chris.neale at utoronto.ca]

I'm not sure what the reviewer has in mind. Therefore I would split  
the response into a number of possibilities.

1. If you simulated a box of those solute chemicals in the absence of  
a large solute, and if those solutes do not phase separate, then the  
collisions will equipartition the kinetic energy. The solutes are all  
small and chemically similar insofar as our forcefields treat them.

2. If one then added a large solute, one may be concerned that the  
thermal fluctuations of the solute will partition separately from the  
solute. This is a real concern and this is why I use Langevin dynamics  
(the sd integrator). But you seem to be not very concerned about the  
thermal fluctuations of the solute, so this, on its own, seems to be  
not a concern.

3. The solvent beside a frozen solute: this seems to be a serious  
possibility of a problem. Did you freeze it entirely? If so, does  
gromacs discount this from the kinetic energy? If so, then all seems  
fine. But if you used harmonic restraints then does that somehow suck  
energy out of nearby solvent molecules? I don't know but I might be  
concerned about that if I was trying to draw kinetic conclusions about  
the system.

4. What are your conclusions? Do they even depend on the equipartition  
of energy? If so, how?

So my base answer is that yes, you do possibly have problems and these  
can all be solved with Langevin dynamics (but then kinetics are  
irrelevant and conclusions can not be drawn). But more applicably,  
might these problems affect your conclusions? and here I honestly  
doubt that you have any problem by coupling solvents together. For  
example, if you coupled a box of water together,
  would you expect problems with energy partitioning? I think not.

Finally, to the reviewer's probable question, you have a frozen or  
restrained solute so likely you are not possibly experiencing a frozen  
ice cube type problem. If the reviewer has a particular problem in  
mind, perhaps they would be kind enough to share that with you so that  
you could address it in specific.

Chris.

-- original message --

Thanks Mark for the advice.

I have just rerun a test simulation with each of my gas species
coupled separately to a thermostat and have got similar values for the
quantities of interest (diffusion coefficients). However I am not sure
that this will satisfy the reviewer without a bit more justification
of why I chose to use a single thermostat for my system.

I remembered reading some gromacs information on the application of
thermostats (on which I based my decision to apply one thermostat for
the entire system). I have just managed to locate it on the FAQ page:

http://www.gromacs.org/Documentation/Terminology/Thermostats

What Not To Do
Some hints on practices that generally not a good idea to use:
?	Do not use separate thermostats for every component of your system.
Some molecular dynamics thermostats only work well in the
thermodynamic limit.  A group must be of sufficient size to justify
its own thermostat.  If you use one thermostat for, say, a small
molecule, another for protein, and another for water, you are likely
introducing errors and artifacts that are hard to predict. In
particular, do not couple ions in aqueous solvent in a separate group
from that solvent.

My system consists of
55 molecules of CO2
38 molecules of N2
3 molecules of O2
One large molecule of MCM-41 (mostly frozen but with mobile surface
groups) consisting of 4284 atoms.

My take on this was that the gas molecules were in small supply
compared to the other part of the system so I should avoid using
separate thermostats. Was I mistaken in making this assumption?

Is there any feeling for what "sufficient size" as referred to above
is? Does anyone know of any papers I could reference in my explanation?

Any comments welcome

Thanks

Jenny






Quoting Mark Abraham <Mark.Abraham at anu.edu.au>:

> On 8/03/2011 2:50 AM, Jennifer Williams wrote:
>> Hi,
>>
>> I simulated the diffusion of small gases (CO2, N2) in a framework   
>> structure which was mostly frozen with some mobile surface groups.   
>> I applied a temperature thermostat to the entire system (i.e I   
>> didn't couple the gas molecules and framework separately). I have   
>> now been asked the following by a reviewer:
>>
>> "Please comment on any artefacts that might arise as a result of   
>> non-equipartition of energies.  For example, what is the calculated  
>>  temperature of each of the gas species and mobile species?"
>>
>> I have tried to explore the temperature of each species separately   
>> but g_energy will only give me the energy of the system as a whole.  
>>  Are there any other tools within gromacs to look at the  
>> temperature  of each species in turn from the md runs I already have?
>>
>> Does anyone have a suggestion as to how I can address the   
>> reviewer's comment and proove that using one thermostat for the   
>> whole system is OK (that is what I am hoping!).  Is there anything   
>> in particular that I should be looking closely at?
>
> g_energy can only report the global temperature as saved in the .edr  
>  file, computed during the simulation. The only way of decomposing   
> the temperature is to save velocities to the .trr file with nstvout.  
>  (Some colleagues had to do this recently.) Then, judicious use of   
> trjconv to take subsets and matching topology hacking (tpbconv   
> doesn't quite work) mean you can use mdrun -rerun on the subset .trr  
>  and matching .tpr to get a group-wise temperature. That can   
> demonstrate whether temperature differentials (ie. heat flows) exist.
>
> If you didn't save velocities, you're out of luck.
>
> Mark