Hi,
My goal is to simulate irradiation under applied stress. I want to be able to simulate multiple PKAs to observe a change in the microstructure. system: carbon atoms with a machine learning potential. periodic boundary conditions.
At first, I tried NPT and I am using an adaptative timestep. The issue is that the temperature takes a very long time to come back to the initial temperature and I am going to run into computational time limit if I want to simulate more than 10 PKAs. Moreover, I noticed that the pressure increases a lot too and doesn’t decrease!
Now, I am thinking of using NVE + langevin thermostat (to have the thermal spike decreasing more rapidly) + berendson barostat (to come back to expected applied pressure). My concern is: does this physically make sense? And is that the way to go about it?
Or should I do NVE + temp/rescale (temperature control) followed by NPT? and repeat?
In summary, I need to find a meaningful but also efficient way of controlling those parameters, because for now I would probably be able to control them well in longer times, which would reduce the number of PKAs I can introduce (so dpa simulated).
And also mixing some ensembles (NVE) with other commands (temp/rescale) make me wonder about the physics meaning and correctness
Any advice would be great
You can always write a restart and continue from there.
You have a significant disruption to a tiny system (compared to a bulk experiment), so it can take a long time to relax. It is a solid after all, right?
What is your objective here: have better statistics by looking at more events, compressing time, or simulate a very strong irradiation?
Never use fix temp/rescale
for any meaningful production simulation.
Yes it is a solid
yes the objective is to simulate strong irradiation event. And my concern is that if I want to do that and claim I am simulating irradiation at a specific temperature and under a specific applied pressure, then I feel I should come back close to those specified temperature and pressure after each PKA event. And I was wondering which way to go about that
The physical approach is to wait. If you have strong irradiation in an experiment, that sample will heat up accordingly. You probably have not estimated how strong an irradiation your sample would have to undergo in an experiment. Length scales and time scales in MD are painfully short. The whole idea of using a thermostat is unphysical unless the coupling is very weak or you want the thermostat to also represent an implicit solvent.
As I mentioned, you can extend simulations with restarts.
If you need better statistics, you can trivially parallelize this. From your original, already equilibrated sample, you perform a piecewise trajectory and save a restart every 10ps or so.
If you get 5 restarts you can then multiply them by making, say, 10 copies of each, assigning a new random velocity with different seeds and simulating those for a few picoseconds. After that you have 50 decorrelated samples that you can run your studies on and run those independently, and each can have one cascade and then relax, restart, and continue.
You can compress time by using a Langevin thermostat instead of a Nose-Hoover thermostat and in both case experiment with the relaxation time, but then you cannot really make any more assessments of the dynamics of your system, since your time axis is unknown.
The overall issue is that the accessible time and length scales are not sufficient to study this problem well in the first place. But suppressing the symptoms is not going to help you much.
You just turn a computer simulation into a computer animation.
Ok thanks.
If I was to use Langevin thermostat, I assume I would use NVE. But in this case, I need a separate barostat, like Berendson? Otherwise pressure isn’t controlled right?
Again, I would look at this from the opposite direction. Should the volume change? And should it change quickly? So it may be more physically correct to just use a larger system, so the impact of the local disorder on the total system is less. Fix press/berendsen is designed for liquids, so I don’t think it is a good choice. In principal, both, the Berendsen thermostat and the Berendsen barostat have known problems and they are exaggerated on solids. In an experiment, the system’s volume will also change slowly, if at all, so you likely have heightened local stresses. I am not enough of a domain expert on judging this for your specific use case.
A big problem is that people often do questionable things, and get them published, too, and sometimes even if those simulation results do not correlate very well with experiments.
At least from the outside, there is not much advice than can be given beyond that the physics should always determine how you set up simulations and not how fast those simulations can run and how much resources you have available without making an additional effort.
I would not use a barostat at all for this simulation for two reasons.
Firstly, from my own personal experience, some ML force fields (especially the faster ones) are poorly tailored for barostatting, as they may not correctly “taper off” near the cutoff, or more generally have the correct long range behaviour. (I am sure most ML force fields do not suffer this, but I have seen some that do.) A barostat attempts to pressurise the system by changing its density, and in doing so it must also change the average number of particles within the cutoff around a typical particle, making it particularly sensitive to any errors near or at the cutoff.
Analytical force fields address this either through the tail corrections of pure LJ terms or through long-range electrostatics (which will dominate LJ dispersion for meaningfully polar or charged systems), but an ML FF parameterised without a long-range “embedding” doesn’t really have a principled way to address this. Far simpler to perform a constant volume simulation.
Besides, similar to what @akohlmey said, why shouldn’t a local vacancy temporarily have a different local pressure from its bulk surroundings? A particularly silly thermodynamicist might even posit that the material hasn’t truly re-equilibrated to pressure until … atomistic reorganisation has simply refilled the vacancy and restored the old lattice. Possibly on geological timescales. Neither feasible nor interesting.
Even without invoking a thermodynamicist’s ghost, solid materials can sustain inhomogeneous stresses for quite a while and remain recognisably stable and intact. The Prince Rupert’s Drop is a particularly dramatic example.
So – unless you know that your material of interest undergoes pressure-driven reorganisation at relevant pressures and timescales (in a diamond anvil cell, maybe) I would leave out the pressure altogether. Simply start with the known bulk density and delete atoms as you like until your computer time runs out
Check D.Tramontina et al., Comp.Mat.Sci. 227 (2023) 112304