NVT with gas reactions

[Zack_Gainsforth]

Hi All,

I’m simulating some organometallic reactions in gas at various temperatures (in this case let’s use 300K). In order to produce the gas I started with a solid that had the right elements in it (about 400 atoms, Fe, C, and H). I heated it up to 10000K to produce a gas, and then cooled it back down to the desired temperature. (Maybe there is a better way to do this! I’m all ears!) I’m using a Reaxff potential.

I then hold at the desired temperature for a few nanoseconds:

units real
timestep 0.25
…
print “Equilibrate at 300K”
fix 2 all nvt temp 300 300 $(100.0*dt)
run 10000000
unfix 2

What I notice is that molecules tend to freeze after a few thousand time steps. Occasionally, many molecules will pick up speed and subsequently freeze again – also on the same time scale. I made an animation of the simulation to show what I mean:

Plotting the temperature from the log file shows that the average temperature is steady mostly between 200 and 400K which should be reasonable given that I only have 400 atoms in the simulation. However, there are sharp spikes sending the temp up to 800K, or down to 150K and these seem to correspond to freezing or jumping of molecular species:

I suspect the following is happening (though obviously I don’t know, or I wouldn’t be writing this post):

Molecules are undergoing reactions and the energy added/removed by the reactions are sufficiently large to overload the thermostat’s ability to keep an equilibrium. If this is the case, then simply running longer won’t help since I would expect reactions to continue to occur and continue to create short-term non-realistic temperature fluctuations. I believe my timesteps are sufficiently short. I suspect what I need to do us run a much bigger simulation (maybe 10000 atoms?) so that the differences in energy from each reaction are not enough to influence the thermostat dramatically. Is this correct? Is there another tactic to stabilize the simulation without breaking my CPU, or am I doing this correctly and I should now scale it up and move it to a server?

Cheers,

Zack

[sjplimp]

You could try using fix langevin instead of fix nnt. A Langevin thermostat is applied on a per-atom
basis and not a collective basis like fix nvt. So it might cool down hot atoms more efficiently.

Steve

[Zack_Gainsforth]

Hi Steve,

I’m running it with Langevin and sure enough it is much more stable now. It looks like the solvent term causes trajectories that seem odd (molecules changing direction without a collision), but perhaps if I run long enough I may sample enough configurations to find all the molecular reactions. Once the reactions settle down, I may also try again with nvt and see if I get a similar result.

This was helpful, thank you!

Zack