Out of atom range error in hydrated model

Hi all,

I am currently simulating the polymer-water system. While the dry model could work well, adding water molecules leads to the crash. The error usually came at the first few lines or just the first line of the run so sometimes I am unable to visualise the trajectory and locate the problematic atoms. I double-checked the force field coefficients and the initial geometry. They are good and there is no overlapping in the system. I also do the minimisation in the beginning.

The thing I could not understand is that not all the packing would fail. I created 5 different packings for my system. The only difference between these packings is the initial geometry since they are randomly distributed in the unit cell. I followed the same protocol to develop the polymer system and the polymer-water system. Sometimes only one out of five of the system could work but sometimes all of them failed with the same error. I think the successful one is very lucky and I can’t replicate the success. I am wondering why would this happen and if the initial geometry would be the key issue.

Another thing I’d like to mention is I tried to test some other polymer systems (different kinds of polymer, different functional groups). They almost all work well with the water. So I considered the force field coefficients may be one of the issues for the current polymer system. But if the cofficients are not good, it is difficult to explain there is still one sucsessful packing with those cofficients.

I would like to know if there is some tips or tricks I could consider rather than remodelling the whole system from scratch. Since there is one lucky case, the format and syntax of my input script is at least good enough to start the running.

Any comment or suggestion are highly appreciated.

rescale.in (1.4 KB)

It is difficult to give specific advice without being able to reproduce the simulations.

My guess would be that you may be deforming the system too quickly for the chosen settings and that your neighbor list settings are too aggressive. So either a shorter timestep and a longer run or using a neighbor list delay of 0 instead of 2 could help.

Hi Axel,

Thank you for your comment. I try to change the neigh_modify setting by ‘delay 0’ and reduce the timestep from 1fs to 0.5 or 0.2fs. It doesn’t eliminate but just delay the error. While other polymer systems work well with the water, I think there are some potenial issuses in my current system and the water is likely to highlight those issues.

It seems there is no way to debug the system by simply modifying the input script given my case. Thank you for your time.


Pay attention to the “Note” in the documentation for fix deform:

For non-equilibrium MD (NEMD) simulations using “remap v” it is usually desirable that the fluid (or flowing material, e.g. granular particles) stream with a velocity profile consistent with the deforming box. … One way to induce atoms to stream consistent with the box deformation is to give them an initial velocity profile, via the velocity ramp command, that matches the box deformation rate. This also typically helps the system come to equilibrium more quickly, even if a thermostat is used.

I suspect you have one (or several) H-O bonds across the shear boundary in your initial configurations that failed, and the configuration that ran successfully did not have H-O bonds initially across the shear boundary. So the failed configurations immediately had large velocity mismatches across the boundary, whereas the successful configuration came to a steady state before any bonds crossed the boundary.

Try adding in a velocity ramp – alternatively, you can get to a suitable steady state by combining successive runs, the first with a very small shear rate and the last with your desired shear rate.

1 Like

The other thing to verify is whether this happens specifically with TIP4P. Since you use the /tip4p styles, your input data file should be compatible with any 3-point water force field (I would try SPC/E and generally avoid hydrogen LJ centres), and redefining the water coefficients (LJ, bond and angle) and charges (using set type X charge Y) should give a quick check.

The other thing to verify is whether this happens specifically with TIP4P

Or asked differently: for which specific water model was your polymer parameterized?

Force fields parameters follow specific rules and conventions when they are derived and thus those parameters and rules would have to be changed if you change the water model or use parameters for an implicit solvent for explicit water.