Unexpected GPUMD Results

[fkavanagh03]

Hello everyone,

I have been running simulations in GPUMD to calculate the surface energy of Ice Ih (0001) using a NEP potential. I am computing surface energy as:
γ = (E_slab − N × E_bulk_per_molecule) / (2A)

where E_bulk_per_molecule is obtained from a separate bulk calculation of the same lateral dimensions, E_slab is the minimized slab energy, N is the number of molecules, and A is the cross-sectional area.

My workflow for both slab and bulk are:
• Coarse FIRE minimization
• NVT annealing (300K → 0.1K)
• NVE
• Fine FIRE minimization

The lateral dimensions are fixed at 180×180 Å and I vary the number of unit cells in the z-direction. In theory, surface energy should not depend on slab thickness once the slab is thick enough to decouple the two surfaces. However, I am seeing non-monotonic behaviour with no sign of convergence:

of unit-cells thick Surface Energy (meV/A2) Tot eng (eV) # of atoms

6 (bulk) 6.791 23621.495 124416
6 (slab) 6.791 23240.4130 124416
12 (bulk) 8.198 47736.763 (x2.021) 248832 (x2)
12 (slab) 8.198 47276.7096 248832
18 (bulk) 3.038 72259.783 (x3.059) 373248 (x3)
18 (slab) 3.038 72089.3214 373248

I am not sure why this is happening. Could this be a minimization convergence issue, or is there something else going on? Thanks!

[brucefan1983]

why you need NVT anealing? Could 300 K already make ice to melt?