Hi all,
I am trying to perform NPT ensemble simulations of crystalline TATB, which is triclinic using a 6a x 6b x 6c super cell (where a,b,c are the edge vectors of the triclinic unit cell). The fix NPT command that I use is of the form below:
fix 1 all npt temp Tstart Tstop 1000.0 tri Pstart Pstop 1000.0 drag 3.0 flip no nreset 10000
The pressure is constant throughout the simulation and the temperature is ramped up to perform super heating at a rate of 2x10^11 K/s. For simulations at 1000 atm (Pstart and Pstop) , the super cell appears to lose periodicity at 950K. The variation of box parameters and visuals of the simulation cell are shown in the attached file (slide1).
From the behavior of the box parameters and visuals of the simulation cell, it appears that the box shape changes, the material loses periodicity and develops an interface. Running multiple realizations (with different seeds for drawing initial velocities) does not change the behavior.This also happens in simulations at 1 atm.
When using a larger 10a x 10b x 12c super cell, the behavior manifests at 800K (1000 atm). The details for the 10a x 10b x 12c super cell is in the file slide2.
Has anyone seen this behavior before with triclinic cells? I would appreciate any help to understand why this is happening.
Thanks,
slide1.pdf (636 KB)
slide2.pdf (735 KB)
Have you dumped snapshots and examined the box size/shape
as a function of time in the LAMMPS dump file? There is nothing that prevents a material
from partitioning into 2 phases or creating an internal surface
within a periodic box. Esp if your box dynamics are bad at high
temp or press.
Steve
Some comments:
1-) Why don’t you explain what you expect the system to do? I personally have no idea if TATB goes a phase transform or else at those high temps/rates etc.
2-) Have you checked that the units fit the potential parametrization, the timestep is correctly chosen, etc?
3-) What do you mean by the box loosing periodicity? Do you mean loosing crystalline structure?
4-) What happens if you go with the “aniso” option instead of the “tri” one? Again, are you expecting your box tilt factors to change?
Carlos
Thank you Carlos and Steve, for the comments.
in answer to Carlos’ comments:
-
I expect the system to melt, since the experimental melting point (for TATB, pressed powder) is 723 K and the temperatures that my system is heated up to is more than the expected super-heating limit of 1.2 times the melting point. There could be a solid state phase transition before melting in TATB, but that is yet to be found out.
-
The units are set to “real” and it is consistent with the potential. The time step is 0.25 fs. This value is chosen so that there are roughly 40 time steps per period of the highest frequency eigen mode in crystalline TATB, which is the antisymmetric amine stretch (time period of roughly 10 fs).
-
The box remains crystalline. I had attached two files with my previous email which shows the snap shots of the simulation cell and the variation of box parameters. One of them is attached here again with this email. In the snapshots shown in the file named slide1, once can see that there is a “gap” which develops between the super-cell and it’s periodic images. This is what I meant by “losing periodicity”.
-
My understanding is that the “tri” keyword should be able to accommodate phase transitions (if any) with changes in the unit-cell angles. Also, as one can see in the attached file, the tilt factors do change with increase in temperature and this indicates that “tri” is the appropriate setting.
in answer to Steve’s comments:
The snap shots and changes in the box parameters are shown in the attached file. One of my concerns is that, for the systems that I had tried, the partitioning always happens at the boundary of the periodic box. It appears as if the box size/shape changes but the crystal does not “fill” the box. And as I had explained in the previous email, it happens at different temperatures for different sized cells - this gives an indication that it is probably not physics.
I am trying to see if there is a way to look at it in detail and understand the cause. Any pointers will be very much helpful.
Further details of the potential and the simulation, if it helps:
Non-polarizable, flexible force field given in J. Chem. Phys. 139, 074503 (2013). All bonds, angles and improper dihedrals are modeled as harmonic oscillators. Dihedrals are single cosine terms. Non-bonded interactions are modeled using the Buckingham potential and Coulomb interactions with fixed partial charges. Intra molecular non-bonded interactions are excluded. Kspace style: pppm with relative error in forces set at 1e-7. LAMMPS version: Dec. 7, 2013.
Thanks, again.
slide1.pdf (636 KB)
I have looked at the pic {Didn’t do it before 
}
That feels to me as if the atomic motion is getting decoupled from the box motion. Almost as if the box was lagging behind when chasing the atoms. Have you look at the cell fluctuations while varying the drag keyword? What happens with drag=0.0? Have played with the damping parameter of the barostat as well?
Carlos
Something else to add. This is your research and my time is limited but I looked (FAST) at the original paper where the FF parametrization is done and don’t see anything about high temps there. 400K was the max they tried. You are blasting 1000K to a molecular crystal that seems to use a parametrization that may not work in the range you are trying to cover. Anyways, something to keep in mind. Disregard my comment if I missed something due to the rush.
Carlos
Thanks, Carlos. I’ll check the effect of the drag parameter and get back.
You are right that the potential has not been tested at high temperatures. It might also contribute to this issue, but my thought was that simulation would probably be inaccurate at high temperatures and not just ‘fail’ as in this case.
Yes, the mess should not be only the result of a bad potential but could help to enhance the problem. Pure LJ solids tend to be quite brittle and yours is not much different from such cases.
Carlos
Follow up:
I tried changing the drag parameter and the barostat coupling term. The loss of periodicity (explained in the previous emails) occurs irrespective of the value specified for the parameters.
I have attached plots of the evolution of Yz tilt parameter vs time, to show the effect of changing the drag parameter (effectofDrag.pdf) and the coupling constant (effecofCouple.pdf). Only this tilt parameter is shown since the ‘divergence’ which correlates with the loss of periodicity can be clearly seen. The onset of divergence is marked using arrows with colors corresponding to the colors of the curves. The temperature ramp is also shown in the plots. It can be seen that the temperature at which the problem occurs is sensitive to the change in the coupling and drag parameters, but the problem persists nevertheless.
Thanks, again for taking a look and any thoughts/comments are much appreciated.
effectofCouple.pdf (140 KB)
effectofDrag.pdf (129 KB)
Can you plot as well your pressure fluctuations?
Carlos
A bit more to add. I am afraid I don’t know well the inner workings of the lammps fix npt. However, keep in mind that you are trying to simulate a special system in the sense that molecular crystals have very anisotropic normal modes (the soft modes could be orders of magnitude larger than the high frequency ones. And so NPT MD would be sensitive to this fact. This is actually a challenge also for techniques like metadynamics where one needs to make the width of the bias gaussian potential a function of the eigenmode wavelength in order to avoid having to deposit a whole of gaussians. Anyways, the two plots you showed agree with my statement in the sense that retarding the box motion delays the breakup point.
Here is something else to try: I have the feeling the the unit cell of your crystal is such that you have the molecule (or molecules) located inside the cell and no bonds initially cross the periodic boundaries (you probably used the replicate command to create supercells). Given you have observed the breaking happening right at the boundaries I wonder what would happen if you start by using a cell with origin so that the molecule(s) has(ve) bonds crossing the periodic boundaries (not enclosed in the cell). If that is not your case then disregard my comments and wait for one the guys who really know the stuff (core developers) to provide some insight. I may well be crying wolf for no reason.
Carlos
Thanks (again) for the comments. I have attached two files which shows the the system pressure and the Sigma_yz as reported by LAMMPS, for the different values of drag (effectofDrag-Stress.pdf) and couple (effectofCouple-Stress.pdf) parameters that I have tried. The plots are moving averages, where the average is taken over twice the value of the barostat coupling parameter. The arrows correspond to the location of the apparent ‘divergence’ in the Yz tilt parameter. The plots are noisy but some features can be observed in Sigma_yz at these locations.
Regarding the periodic boundaries, my initial configuration does have bonds straddling the periodic-boundaries. But it will be still be interesting to see what happens when the cell is translated to the middle of the box. I’ll get back once I test that.
effectofCouple-Stress.pdf (473 KB)
effectofDrag-Stress.pdf (514 KB)
I’m still not clear from your various plots on
why you think the box size is changing (in an incorrect way)
or some kind of non-periodicity is resulting.
You are running NPT with box flips turned off. Your plots
of various tilt factors show no sudden changes in
tilt factor (e.g. due to a box flip) as they shouldn’t b/c you
turned flipping off. What your are calling “non-periodic”
is based on the molecules pulling away from the box boundaries,
not anything to do with the box shape itself?
Are you driving the system (via NPT and the temperature) to
a tilt that is far beyond where LAMMPS would normally
want to flip? If so, I would try running while allowing flips
and see if you get different behavior. Conceptually there
should be no difference between the two (aside from round-off) and flipping will
maintain a less-skewed box which will be better computationally.
It just might be harder to analyze/visualize your results. We only
allow the non-flip for that reason.
But it’s possible there is some setup code in the neighboring
that assumes the box will not become too skewed, and if
it does, then not enough ghost atoms are being communicated.
I’d need to check into those details. If you get different
answers with flip vs no-flip that would be a good hint there
is some assumption which we didn’t realize large skews would
break.
Steve
Steve,
Yes, by “non-periodic” I mean that the molecules pull away from box-boundaries. As I had shown, the box parameters also show a change close to the point where this ‘non-periodicity’ is observed.
The only reason why I used ‘flip no’ was for ease of analysis, as you had mentioned. I’ll perform a run with the flip allowed and see if it changes anything. Will get back.
It appears like allowing LAMMPS to flip the box changes the results. Steve guessed it right - I am driving the system to a point where LAMMPS actually wants to flip the box and the results appear to be different beyond that point. I have attached two files which shows the details of two simulations - with and without flip. The drag is set to 1.0 and couple is set to 1 ps in both the simulations. The plot shows the time evolution of the tilt parameter Yz, which is driven beyond the 0.6*Ly limit (for NPT) for flipping with some snapshots, as explained below.
The box-flipping can be easily observed in the file FlipYes.pdf. The point is also marked with arrows. I have shown snapshots of my super-cell with periodic-box and periodic-images drawn using VMD. The snapshots on the left side are two projections of the super-cell 0.575 ps before the flip. In one of the projections, there is a gap between the periodic images where molecules appear to be pulling away from the box-boundaries. Looking at the values of the Yz tilt parameter, I can see that this occurs when 0.5Ly < Yz < 0.6Ly. The images on the right shows the same projections 0.675 ps after the flip and it can be seen that there are no ‘gaps’ once the flip is performed.
The corresponding plot and snap-shots when ‘flip no’ is used, is given in the file FlipNo.pdf. The arrows show the time when flip would have happened. The snapshots are taken at same time steps as those shown in FlipYes.pdf. Comparing the snapshots on the right side in both files (when Yz>0.6*Ly), it can be seen that the ‘gap’ persists in the case with ‘flip no’. Snap-shots taken later in the simulation also show the gap when flip is not allowed, which suggests that there is some inconsistency between the two cases.
FlipYes.pdf (320 KB)
FlipNo.pdf (898 KB)
So what is the maximum skew your simulation gets to with flips turned off?
Or can you tell what the skew is at the point your results start to go bad?
What we will probably do is put an error check in the code to say you
have skewed too far and ghost atoms will not be fully correct.
I.e. the solution for you is not to turn flipping off.
Steve
Steve, as far as I have looked, I can see this problem occurring when the tilt >0.5*L.
In my case this translates to Yz > 0.5Ly. Hence, I observe it only for a short while when 0.5Ly < Yz < 0.6*Ly for NPT simulations with flips allowed.
Thanks again, for the help. I’ll go with flip turned on.
I looked at how we do the tilt/flip code. The comm of ghost neighbors
is independent; it should work with any tilt. The attached test
script works for any bizarre tilt factors you want to specify in the
region prism (mutiples of 2), and always gives the same number
of neighbors and correct thermo output and dynamics.
So basically I still have no idea what is going wrong
in your simulation. The difference with flips on/off suggests
it’s a LAMMPS problem, but it could still be some issue
with your model or your interpretation of the output.
Steve
in.test (471 Bytes)
Two ideas to test to see if your model
is causing the problem.
a) you are also applying external pressure
in the x direction, which will affect the box
length in x. If that enforces a box length
that is non-commensurate with the number
of molecules in the x direction, then the
system might want to separate. Maybe
that happens at a certain value of tilt.
What happens if you change the external
pressure in x, maybe make it 0.0?
b) what if you shift the box origin (again in x) so
that it doesn’t start between 2 molecules, but
in the middle of a column of molecules?
Does the separation still happen, but maybe
in a different location, like in the middle of the
box? Again, that would indicate there is nothing
special about the box edges and PBC.
Steve
I translated the cell by 3*a, where a is the unit-cell edge vector aligned with x with same box definitions as before. In this case, I observe the separation in the middle of the box, as shown in the file TranslatedFinal.pdf. The initial configuration is shown in TranslatedInitial.pdf. If it was something with the neighboring/PBC I would expect it to happen again at the cell boundaries. I am quite puzzled…
About the system pressure: All the plots that I had shown was for simulations done at 1 atm. I start off my heating at 670 K, with the super-cell created using lattice parameters obtained from (an independent) NPT simulation at 670K, 1 atm. I can check 0 atm, but I am not very clear as to why the system would separate when (super-) heated in NPT instead of accommodating the thermal expansion.
Another thing that I haven’t checked yet is the effect of neighbor skin distance. I am using a distance of 1.0 Angstroms, which is smaller than the default for ‘real’ style units. This setting worked well for my low temperature simulations but may be not for the current simulations.
I will get back once I perform the checks with the zero pressure and neighbor skin.
TranslatedFinal.pdf (113 KB)
TranslatedInitial.pdf (24.5 KB)