isotropic vs anisotropic pressure coupling for organic crystal pressure control

Hello everyone,

I would like to get some advice on the subtleties regarding anisotropic vs. isotropic pressure control when simulating crystals in lammps. I am attempting to measure the enthalpiy of an 8x4*8 alpha glycine crystal (1024 molecules). Using the geometry obtained from the CSD, and after applying the appropriate rotations to cartesian frame, I found that the experimental cell lengths are lx = 40.8160 ang, ly = 47.8836 ang, 40.3593 ang, with a tilt factor of fz = -16.1431, with all other tilts zero. These were the initial lengths I used in my input data file for lammps.

I initially ran an equilibration for 2ns in the npt ensemble, with the following fix applied:

fix 2 all npt temp 330.0 330.0 2000.0 aniso 1.0 1.0 1000.0

I noticed that half way through the run, the total energy dropped suddenly. I checked the geometry of the lattice, and have plotted the box lengths (lx,ly,lz) below in aniso.pdf. As you can see, the box lengths experience sudden changes, after having reached a plateau.

I suspected that this was connected to the pressure control, and ran the same simulation, but now with isotropic pressure coupling. The box lengths are plotted in iso.pdf. As you can, there is so abrupt structural change when maintaining the pressure this way.

At this point, I cannot tell which type of pressure control is more appropriate. The abrupt change with anisotropic pressure coupling to me is unsettling. However, it could just be that isotropic coupling is preventing a structural relaxation. Furthermore, the isotropic pressure coupling reproduced the equilibrium box lengths better than anisotropic coupling, which shows a flip in the ranking of the lx and lz box lengths. The force field I am using is GAFF + CNDO charges (D.W. Cheong, Y.D. Boon, Comparative Study of Force Fields for Molecular Dynamics Simulations of Alpha-Glycine Crystal Growth from Solution). In the cited work, its not clear what coupling the authors used with performing bulk crystal simulations.

Any insight regarding the selection of the coupling method, or comments about the physics I am observing, would be greatly appreciated.

aniso.pdf (21.8 KB)

iso.pdf (20.9 KB)

Dear Connor,

Can you provide the values for the tilt factor as well? What I think may be happening is that your simulation box is being flipped, which could potentially cause problems with the computation of the strain tensor due to it’s dependence on a reference box (h0).

Thanks,
-David

Hi David,

Thank you for your response. I have attached the isotropic and isotropic xz tilt factor data in the attached tilt.pdf image (the other tilts are zero). As with the lx,ly,lz box lengths, the xz tilt factor shows an abrupt change when using anisotropic pressure coupling.

tilt.pdf (17.3 KB)

Thanks, it seems that my idea was off base since the change in tilt doesn’t indicate a flip. I suspect that there is a structural relaxation in your system involving tilting that is suppressed in isotropic relaxation. Perhaps you can visualize your system to get a better idea about what is going on.

-David

The original crystal structure appears to be stable under NPanisoT for quite a long time, albeit with lattice constants that are a little off. This would suggest that you have implemented the potential correctly in LAMMPS, although more evidence would be necessary to feel confident about this (compare with results published in Cheong and Boon). The problem you are seeing seems to be a structural transition to a more stable phase. You can examine this by calculating PE+pV for the different states. Since the changes in V appear small and p=1atm is small, the pV term is probably negligible. so it is likely that the PE for the new structure is substantially lower.

You should try changing aniso to tri. Sometimes suppressing box tilts can make a crystal less stable, although in generally the opposite is true. You should also examine what happens in a smaller cell. Generally, smaller cells are more stable than larger cells.

If you can’t make the instability go away, it might be okay to suppress it by using NPisoT. However, you should worry about what this instability will do in the context of your intended application.

Aidan

Thanks everyone,

I have tried visualizing in VMD for the past couple days, however the resulting series of image are not particularly revealing.

Aidan, your intuition is correct. I checked the total energy first, and noticed that it dropped precipitously. This is what made me suspect a structural relaxation. I have attached a plot of the total energy vs. time for two type methods. As I am trying to measure enthalpies, the results differ more than I would like between the two methods. The question in my mind is whether the structural relaxation is indeed something physical, or is rather an artefact of the algorithm employed. I.E., is aniso always regarded as more “physical.” If so, the instability isn’t really an instability, but really a necessary part of the equilibration. I was leaning towards the iso algorithm, as it did a better job of reproducing the equilibrium box lengths.

I was thinking about trying the tri command a couple days ago, however I was slightly confused about what pressure to specify for the tilts. The manual states:

"Using “tri Pstart Pstop Pdamp” is the same as specifying these 7 keywords:

x Pstart Pstop Pdamp
y Pstart Pstop Pdamp
z Pstart Pstop Pdamp
xy 0.0 0.0 Pdamp
yz 0.0 0.0 Pdamp
xz 0.0 0.0 Pdamp
couple none"

I could not decide if the pressure for the tilts (xy,yz,xz) should be 1.0 atm (the total pressure I am using for my simulation), or if it should be kept at zero, as in the manual example. 

tote.pdf (17.8 KB)

The fact that the aniso PE is higher than the iso PE is a clear sign that your aniso structure is damaged in some way, probably during the initial relaxation. Instead of burning vast numbers of flops on a large system running for many nanoseconds, you should take a few steps back. Start out by doing some short simulations with the smallest acceptable supercell. Once you know the true lattice constants for that system at P=1, you can start your big system with the correct cell dimensions. In that way, you will avoid any strange behavior caused by a fast initial relaxation.

Aidan

For anyone who comes across this thread for help in the future,

Aidan’s insight was indeed correct. Asking the system to equilibrate in the NPT ensemble using anisotropic coupling right off the bat was simply too much to ask. The proper equilibration procedure is (1) NVT ensemble (2) NPT ensemble isotropic coupling (3) NPT ensemble anisotropic coupling. To verify this, I have attached a plot of total energy vs. time. The first 2ns were an equilibration in the NPT ensemble with isotropic coupling. At 2 ns, the fix was switched to anisotropic coupling. The system total energy relaxes, and shows smooth dynamics for the remainded of the run :slight_smile:

Thanks for your help everyone!

tote_iso_aniso.pdf (23.5 KB)

1 Like

Sorry, just realized the legend of the proceeding graph was off. The blue data is isotropic coupling, and the red data is anisotropic coupling. Sorry for the confusion!