[lammps-users] Bad potential or Bad simulation skills?

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I can’t reproduce a MD simulation result regarding to thermal expansion published in a paper and I am not sure it is due to bad inter-atomic potential I realized or bad simulation skilled employed during my simulation.

I want to use LAMMPS to study various properties of uranium dioxide. Apparently the potential is not covered in LAMMPS. But the analytical form of the empirical inter-atomic potential is presented in the paper mentioned earlier. I tabulated the potential and used the pair_write to verify whether I did it correctly.

I want to validate the tabulated potential by modeling the thermal expansion under NPT ensemble and comparing it to the published result. But the simulation result is completely messy and has no pattern at all.

Could anybody give me some suggestion about where I did wrong? Another thing that is worthy of mentioning is the paper claims the employment of Ewald summation which is absent from mine. I wonder could this be the key factor that I am missing.

BTW, I am using LAMMPS MAY-01-2010.

My input file is listed below:

UO2 which is iso-structural with fluorite

units metal

boundary p p p

atom_style atomic

variable T equal 10

lattice fcc 5.47065

region box block 0 4 0 4 0 4

create_box 2 box

lattice fcc 5.47065

create_atoms 1 box #Uranium

lattice custom 5.47065 a1 1 0 0 a2 0 1 0 a3 0 0 1 basis 0.25 0.25 0.25 basis 0.75 0.25 0.25 basis 0.25 0.75 0.25 basis 0.75 0.75 0.25 basis 0.25 0.25 0.75 basis 0.75 0.25 0.75 basis 0.25 0.75 0.75 basis 0.75 0.75 0.75

create_atoms 2 box #Oxygen

timestep 0.002

thermo 1000

pair_style table linear 1200 # Use tabulated potential

pair_coeff 1 1 Yakub.table Uranium 8

pair_coeff 2 2 Yakub.table Oxygen 8

pair_coeff 1 2 Yakub.table Dioxide 8

mass 1 238

mass 2 16

neighbor 0.5 bin

neigh_modify every 5 delay 0 check yes

velocity all create $T 825577 dist gaussian

fix all all nvt temp $T $T 0.01 # equilibrate at low temperature and then heat it up

thermo_style custom step temp vol press

dump all all custom 1 dump.atom id type xs ys zs

run 5000

unfix all

fix all all npt temp $T 4500 0.01 iso 0 0 0.5 drag 0.0 # thermal expansion

run 2000000

Any suggestion would be highly appreciated.
Thank you for your attention.
Bin

2

You can’t compare results with and without using long-range electrostatic calculations; the results are just not going to be comparable to one another; this is especially true if the potential was parameterized using long-range summations explicitly.

While the choice of long-range summation technique (Ewald versus PPPM) should not make a difference, you absolutely should see a difference between using long-range summations and ignoring them.

Note that this does not preclude the possibility that you also have errors in your simulation script; however, it also explains why even a “correct” script would still show disagreements with the earlier work.

–AEI

3

Another point - have you looked/plotted the potential
you are generating with pair_style table, by using the
pair_write command? You would want to be sure
the interpolation you are doing looks good.

Steve

4

Another follow-up to Ahmed's point. It doesn't look like your input
or potential define/use charges at all. So you might ask how
good can an oxide potential be that doesn't use charges. It may not
be any good for thermal expansion effects.

Steve

5

Dear Sir
Thank you very much for all your comments and valuable suggestions.
I did use “pair_write” to verify the tabulated potential and it is correct.
I have made an adjustment through taking the charge of the particles into account and introducing Ewald summation into the simulation.

My following question is after setting the charge of atoms via “set” command LAMMPS indicates my system is not charge neutral. But I suppose the total positive and negative charges should be equal since the multiplications of respective ionic charge and corresponding population are same. I went through the history of mailing list and found a similar case. It was mentioned reading-in the model through “read data” could avoid this warning. I wonder which software I should use to build the “read data” file.

Another question is regarding to how to combine tabulated potential and Ewald summation. I tried to include or exclude the columbic potential in the tabulated potential file. But neither case gives me the correct results.

Any hint or comment would be highly appreciated.

Thank you very much.

Bin

2010/5/31 Steve Plimpton <[email protected]>

6

The error not charge neutral simply means the sum of charges
is not 0.0. Should be easy to check your inputs to see why that
is the case.

Steve