Thermal Expansion using AIMD

[Fermi1976]

Dear VASP Forum,

This message is a continuation of the issue I am encountering with thermal expansion using MLFF. Please refer to the previous post: https://www.vasp.at/forum/viewtopic.php?t=20163

I have been investigating how the c-axis of a 6x6x2 graphite supercell changes with temperature using AIMD (NPT ensemble, TEBEG=20K, TEEND=420K) and various Van der Waals corrections. Specifically, I used the PBE+ functional (IVDW=20, IVDW=263), SCAN+ (rVV10, IVDW=12), and R2SCAN+rVV10. As shown in the first figure, there is a significant variation in the lattice constant, with SCAN+rVV10 exhibiting the least change compared to the other functionals.
1. If AIMD is unable to accurately predict the lattice constant's behavior as a function of temperature, then training the force fields would likely be ineffective. This could be the issue I am encountering in my training, unless there is something wrong with my setup.

Additionally, I performed AIMD calculations with fixed temperatures (TEBEG=TEEND) at 50 K and 300 K. In both cases, the lattice constant decreases as the MD step progresses (See figs 3 & 4).
2. I am attaching the INCAR and the ICONST files for reference. Any suggestions or recommendations to improve the lattice constant behavior as a function of temperature would be greatly appreciated.

INCAR_ICONST.zip

Thank you,

[jonathan_lahnsteiner2]

Dear Fermi1976,

Your calculations seem to be alright in my opinion. It just looks like that you are using a starting structure which is not equilibrated at the conditions you are running your MD simulation. You have to do a equilibration NPT run and then start your production run with the box in equilibrium.

All the best Jonathan

[Fermi1976]

Dear Jonathan,
Thank you for your reply. I really appreciate your help.
Maybe I misunderstood your message, but the first figure (NPT calculations between 20 and 420K) tells me that the structure keep contracting instead of expanding.
What should I do to reach equilibrium with NPT calculations?

Thank you for your support.

[jonathan_lahnsteiner2]

Dear Fermi1976,

Yes this is true. I took this into account. It seems to me that your starting volume is very far from the equilibrium volume of your system. Therefore for all temperatures you are simulating the boxes will contract first. You have to determine the equlibrium volume for all temperatures separately. To tell to which volume/ lattice parameter c the cells will evolve is not possible to tell from the data you sent because the MD runs are quite short. I would run the simulations for longer and then check to which values your lattice parameters converge for the different temperatures. Maybe you can also find a better starting box by doing a static calculation with IBRION=1 or 2 and setting ISIF=3.

All the best Jonathan

[Fermi1976]

It seems to me that your starting volume is very far from the equilibrium volume of your system. Therefore for all temperatures you are simulating the boxes will contract first. You have to determine the equlibrium volume for all temperatures separately. To tell to which volume/ lattice parameter c the cells will evolve is not possible to tell from the data you sent because the MD runs are quite short. I would run the simulations for longer and then check to which values your lattice parameters converge for the different temperatures. Maybe you can also find a better starting box by doing a static calculation with IBRION=1 or 2 and setting ISIF=3.

Hi Jonathan,
Thank you for your prompt response.

The starting volume I was using was based on 6x6x2 supercell (288 atoms) optimization. For the optimization process, I used IBRION= 2 and ISIF =3 (see attached INCAR). This resulted in a 2c value of 13.3 A compared to the experimental value of 13.4 A. Do you think the optimized lattice constant is far from equilibrium value?

I was also extending the NPT running (TEBEG=TEND=300K) by adding 1500 MD steps more. The figure below shows the first 1500MD I shared last time and the new 1500 MD steps. The lattice constant 2c is still decreasing, I will add 2000 MD steps more and let you know how the 2c evolves as a function of the MD steps.

One more question, for the MD simulations I tried ICONST with the all LA are either activated or not, do you have better recommendations? please see the attached ICONST file.

Thank you so much for your support and patience.

[jonathan_lahnsteiner2]

Dear Fermi1976,

Even when you did a optimization to obtain a starting structure for your MD run there is no guarantee that you will reach the true ground state. The structure optimization algorithms can easily be trapped in local minima if the starting structure was chosen unlucky, I think from the additional MD data you provided you can see that the 2c parameter starts to equilibrate to a value around 12.85A but additional MD data is needed to confirm this. I would just run the MD for a little longer and check that the 2c parameter really starts oscillating around this new equilibrium value. And then pick one of this structures as a future starting structure for further simulations. In principle you could take the one with the lowest energy and then relax this structure further and use this as a future starting structure.

Your ICONST file looks fine to me if you want to constrain the lattice angles that is. If you don't constrain them like this the box might loose it's characteristic shape which you want to impose. But maybe it is worth a try to run an MD starting from the equilibrated starting structure and try if the simulation by itself will maintains it's shape.

All the best Jonathan