I am running a replica exchange simulation as shown in the attached file below. The temperature range of the replicas are at a high temperature (~700-800 K), and I am getting the “out of range atoms” error, typical of organic systems (with hydrogens) at high temperatures. I am able to run these individual simulations, by themselves, using a very small time step (0.1 fs), but during replica exchange the processors lose track of the atoms.
Is it possible to carry out a minimization after every temperature swap? I’m thinking that the higher densities at the lower temperatures, when suddenly raised to higher temperatures, are resulting in unrealistic (or at least problematic) pair energies, especially for the hydrogens. If it were possible to minimize between each temper exchange, maybe these problematic dynamics could be avoided.
If you have to use a 0.1fs timestep while using fix shake on all bonds involving hydrogen atoms, then there is something seriously wrong with your force field.
I also would have serious doubts about the utility of fix npt and temper/npt at such high temperatures.
And I would not use fix balance on such a system setup.
The forcefield doesn’t require both fix shake and the 0.1 fs timestep, I just left shake in for good measure. Probably not best practice.
My advisor recommended replica exchange at these temperature ranges to try to avoid getting caught in a local minima. Would you mind summarizing your concerns as to why temper/npt wouldn’t be appropriate? Just so I can be informed when I talk to my advisor about it.
I left fix balance in because, in the first 10,000 time steps or so, there are huge load imbalances spread throughout the simulation box. Would it be better to leave it in for the first part and then unfix it, or just take it out entirely?
You wrote you are getting out-of-range atoms, but that those simulations require a 0.1 fs timestep without the replica exchange. That is extremely small for a system with fix shake on bonds. You should be able to run 1 fs timestep with fix shake (even larger at 300K). So the fact that you need to use a small timestep is an indication that atoms move to fast and that means that atoms get to close and that can only happen with bad force field parameters. Are your force field parameters validated for 800K??
If you are having trouble, you need to reduce the number of variables, so that it will be easier to pinpoint the cause of the trouble. Fix npt (and temper/npt) versus fix nvt and temper are complications, same is true for fix balance. The general rule here is: first you have to be able to run it at all, then run it correctly, and only then you can worry about efficiency and minor details.
As for the rationale, this works the other way around: do you have a good reason to use fix npt and temper/npt? My approach is always, don’t use a feature unless you know it is required or you use it for a good reason.
That just means that you are starting from a non-equilibrium initial condition. I would not want to start a parallel tempering simulation from that, but first equilibrate it.
Also, your input seems to be based on moltemplate. You’re going to make your lift massively easier, if you write out a data file with all the force field parameters and charges using write_data with the pair ij option. That will simplify the startup (and make it much faster) and that also simplifies and speeds up debugging.
Fix shake - I see, this is redundant given the extremely small timestep.
Force field parameters validated at 800 K - Good question, I honestly hadn’t even thought of that. I’m using the OPLS distribution from tinker/provided with moltemplate.
Do I have a good reason for using npt or temper/npt - This dovetails with your other comment on fix balance; one of the reasons I use npt is because the box size is changing a lot… but it certainly makes sense to get the highly non-equilibrium starting conformation equilibrated first and then pass that to an npt run or temper run.
Temper vs temper/npt - Yes, makes sense, temper/npt is an added complication to temper. I could strip it down to its most basic elements and see if it works then to isolate the problem like you suggested.