HCACF never decays and Thermal Conductivity in one direction shoots up

GPUMD
1

Hello everyone,

I am trying to simulate a superlattice structure with periodicity 1 and Periodicity 2 using a NEP potential trained on DFT data on GPUMD version 4.0

The potential has a low RMSE and also a good validation/testing plots for forces and energy with a r^2 value of 0.99. So I started doing Molecular dynamics on NEP-CPU implemented on LAMMPS ( As I am more familiar with LAMMPS, I find NEP-cpu + LAMMPS combination easy to proceed but I also use GPUMD for a few other property predictions like SEDs or DOS.

But the problem now is the plots that I have attached for HCACF and Thermal conductivity. Although I run the production run for more than 50ns with a time step of 1fs, the hacacf looks really weird to me with no signs of decaying. And the subsequent Thermal conductivity plot makes no sense.

What could be the potential problem and solution to this issue?

Please let me know if you need more information in this regards.

The input file for performing GK in LAMMPS is shown below.

( I couldn’t attach the input file so I am copying the script here )

dimension 3

#define units
units metal
atom_style atomic
boundary p p p

box tilt large
read_data P3_S311.lmp

pair_style nep
pair_coeff * * nep_P1P2_combined_aligned_Tsweep.txt A B

timestep 0.001 # 0.5 fs
thermo 10000
thermo_style custom step ke pe temp press density lx ly lz vol cpuremain

replicate 20 2 2

fix 1p12 all momentum 1 linear 1 1 1 angular

fix relaxx all box/relax aniso 0.0
minimize 0 1e-6 30000 10000
unfix relaxx

variable T equal 300

velocity all create $T 3627941 dist gaussian mom yes
velocity all scale $T

dump 3 all atom 100000 all.lammpstrj

fix 3 all npt temp $T $T 0.1 aniso 0.0 0.0 500
run 300000
unfix 3

undump 3

dump 4 all atom 10000 nvt.lammpstrj

fix 1 all nvt temp $T $T 1
run 200000
unfix 1

undump 4

reset_timestep 0
dump 2 all atom 5000000 dump_final.lammpstrj

thermo_style custom step ke pe temp press density lx ly lz vol cpuremain

#--------------------------------------------------------------------------
#kappa parameters
#--------------------------------------------------------------------------

variable p equal 6000 # correlation length
variable s equal 14 # sample interval
variable d equal $p*$s # dump interval
fix 1 all nve

thermo 10000

compute   		myKE all ke/atom
compute   		myPE all pe/atom
compute   		myStress all centroid/stress/atom NULL virial
compute   		flux all heat/flux myKE myPE myStress




variable  		Jxtot equal (c_flux[1]-c_flux[4])
variable  		Jytot equal (c_flux[2]-c_flux[5])
variable  		Jztot equal (c_flux[3]-c_flux[6])

fix tot all ave/correlate $s $p $d v_Jxtot v_Jytot v_Jztot start 300000 type auto file J0Jt_virial_all.dat ave running

run 20000000

Thank you

tc
tc689Ă—529 9.86 KB

hcacf
hcacf707Ă—529 10.5 KB

2

Hi @lammps_umd, sorry for the late reply. I haven’t done much work on heat conductivity, but I’ll throw my hat into the ring.

First things first, what material is it that you’re simulating? Have you visually checked that the MD is stable, and that all relevant properties (energy, temperature, simulation cell etc.) are kept stable throughout the simulation?

Second, when checking the accuracy of your NEP model, how have you validated the model beyond computing RMSEs? In my experience, a model that is somewhat unstable can still exhibit quite low RMSEs.

3

Hi @elindgren, Thank you for replying.

The material I am simulating is a superlattice. Yes, I have checked the simulation cell, temperature and energy and they are stable throughout the production run.

While checking the accuracy of NEP model, I test the prediction accuracy of the model to the unseen test dataset and I get r^2 value of 0.99. Also I de check the rdf of structure from aimd and MD and they overlap with each other. Are there any other way of validating the potential? As there are no experimental results available, I am not being able to validate it by predicting the property yet.

However, if I do the simulation at higher temperature, the hcacf decays. It’s only at and around room temperature that the hcacf is not decaying.

4

Other ways in which to validate the model can be to compute e.g. phonon dispersion or energy volume curves and comparing to DFT. You can find examples of how to compute various properties in the calorine documentation under Working with NEP models: Get started — calorine documentation

With that said, it sounds like your model is pretty converged if it agrees with AIMD. Maybe @erikfransson knows more about the specifics of the GK thermal conductivity calculations and how to converge them?

5

@lammps_umd Im not familar with LAMMPS GK so I can not tell from the input file if things are correct or not.

Are you sure the HACF is not suppose to be this long-lived in time? The signal looks very clean and clear and reasonable to me.
Do you know what the expected kappa values are ? Is it suppose to be isotropic? (You have `replicate 20 2 2` and the `kx` values is deviating alot from y and z so maybe something going on here?)

Maybe also worthwhile to try running this directly in GPUMD with GK and see if the strange results persists?

ping @freeriks maybe have thoughts?

6

Thank you for your response. I replicated it 20 2 2 because I want a converged thermal conductivity along x axis only. I couldn’t do that in GPUMD as it displays error like “the thickness is lower than the NEP cutoff” and when I replicate bigger to avoid this, I get unreasonable thermal conductivity ( isotropic and higher than expected but the system is anisotropic )

So I did the calculations in LAMMPS where I can replicate and perform MD.

I could only resolve the hcacf convergence issue at higher temperature but at 300K, it never worked.

I am not sure what the problem is yet.

I wish somebody knows and responds it.

Thank you.

7

If you are using periodic boundary conditions in all direction, it is better to use larger box sizes in the y and z directions. That way, you will make GPUMD run the job happily.

2 Likes
8

Hi lammps_umd, have you gotten any further with this problem?

If it’s still not solved, can you share the results you get when performing the GK calculation in GPUMD?

9

Hi Elindgren,

Thank you for following up.

I haven’t solved the issue yet.

The following are the input and output file and also the gk post-process matlab code.

Please let me know if we can solve the issue from these files.

NEP_gk.m (2.8 KB)

hac.out (269.5 KB)

run.in (482 Bytes)