Computing RDF with tersoff potential

Dear lammps-users,

I am trying to compute the Radial Distribution Function (RDF) for liquid Si by using Tersoff potential (Tersoff, Phys Rev B, 38, 9902 (1988)) but I am not getting the correct peak width and shape. Here is my code:

###################################Code################################
units real
dimension 3
boundary p p p
atom_style charge

lattice diamond 5.2763
region box block 0 1 0 1 0 1 units lattice
create_box 1 box

create_atoms 1 box
replicate 5 5 5
mass * 27.9861

pair_style tersoff
pair_coeff * * …/…/…/…/…/…/…/…/potentials/Si2.tersoff Si
#pair_style lj/cut 8.0
#pair_coeff * * 1 1 8.0

#compute myrdf all rdf 100 1 1
#fix 4 all ave/time 1 1 1 c_myrdf file tersoff2_3000Kl.rdf mode vector
#rerun dump2_tersoff2_nvt3000.dump first 30000 every 10 stop 50000 dump x y z

compute eng all pe/atom
compute eatoms all reduce sum c_eng

reset_timestep 0

thermo 10
thermo_style custom step pe density temp etotal

min_style sd
minimize 1e-10 1e-10 5000 10000

fix 2 all box/relax iso 0.0 vmax 0.001
min_style sd
minimize 1e-10 1e-10 5000 10000
unfix 2

dump 1 all atom 10 dump2_tersoff2_nvt3000.dump

timestep 0.01
fix 3 all nvt temp 3000.0 3000.0 1.0
run_style verlet
run 50000
###################################End###########################

RDF was calculated by running the code second time with changing the pair style, uncommenting all the lines for rerun command and deleting/commenting everything below the rerun line. Final RDF was computed as the average of the data between 30000 and 45000 timesteps .

I have attached the RDF thus plotted, along with temperature profile and the RDF published for Tersoff potential.

Would you please let me know if I am making any mistake in making liquid Si or computing its RDF? Let me know if you have any questions in understanding the code or my problem. Any suggestion or help is highly appreciated.

Thanks for your time,

Regards,
Vishank

plot_tersoff2_3000Kl.rdf.eps (19.1 KB)

Temp_profile_3000K.pdf (31.1 KB)

RDF_tersoff_3000_ref.jpg

You are using a 0.01ps timestep, which is probably too large for accurate results. There are probably other issues too. Before spending a lot of time dealing with the RDF, look at simpler quantities to validate your simulation, such as average energy, average pressure average density.

Dear lammps-users,

I am trying to compute the Radial Distribution Function (RDF) for liquid Si
by using Tersoff potential (Tersoff, Phys Rev B, 38, 9902 (1988)) but I am
not getting the correct peak width and shape. Here is my code:

why don't you first use any of the (tested) potentials bundled with LAMMPS?
the potential that you are referring to is available as "Si(C)" in
SiCGe.tersoff.
please note that those require the use of metal units.

axel.

Oh right, the script uses real units, so the timestep is 0.01 fs, which is far too small. You can probably get away with 1 fs. The simulation is only 500 fs, which is not enough time to properly equilibrate the melt from the crystal.

Dear Axel,

Thanks for your quick response. I tried it with the Stillinger-Weber potential and I was getting the correct RDF but not with the Tersoff potential. Moreover, the potential I am referring to was optimized for pure Silicon in the original paper of Tersoff.

Let me know if you find something else.

Thanks for your time to answer my query.

Regards,
Vishank

Dear Aidan,

Thanks for your quick response. I have left the system for 20 ps as well and got the same rdf, I reduced the run time for the example and the temperature profile was stable after 500fs so thought its okay. I can try with 1 fs and holding it for longer time. Moreover, I have computed other properties i.e. elastic properties, lattice constant, density and formation energy, which are coming in line with the published results, except the RDF.

Let me know if you find something else.

Regards,
Vishank

The thermostat tdamp is 1 fs. This is far too small and is probably producing artificial resonance. A typical value is 100 fs. These poor choices for simulation parameters are very basic errors that are easily avoided by reading the literature on MD simulations of liquid Si.

Aidan

also, it doesn't make much sense to me, to relax the box at
minimization (i.e. at 0K) and then continue the simulation with fix
volume at 3000K. better to set the density you want for your liquid Si
right away.

axel.

Dear Aidan,

I tried 5 tdamp values: 1, 0.1, 10, 100, 500 and tstep of 0.1 and 0.01 fs to make sure my parameters are correct. For time step 0.1 and 1 fs, I held the system at 3000K for 500k timesteps, i.e. for 50 and 500 ps, respectively, but the system remained in the crystalline phase when dump file was observed using Ovito even though the temperature profile was stabilized at 3000K.

However, for tstep=0.01 fs crystal symmetry got broken and atoms were arranged in random order like amorphous or liquid phase . I did the rdf calculation again for this combination of tstep =0.01 fs and tdamp = 100 but still getting the same rdf.

I will try to change the density to the liquid density and melt the system without relaxing for the box size as per Axel’s recommendation. May be the wrong density is changing the nearest neighbors distance.

Thanks for your help, let me know if you find something else. Will post my findings with fix liquid density.

Regards,
Vishank

Dear Aidan,

I tried 5 tdamp values: 1, 0.1, 10, 100, 500 and tstep of 0.1 and 0.01 fs to
make sure my parameters are correct. For time step 0.1 and 1 fs, I held the

this is all a massive waste of time. you are trying to replace actual
learning and understanding where parameter choices are coming from by
brute force trying out settings.
anything less than a timestep of 0.5fs is going to be rather pointless
and - as aidan already pointed out - 1fs should be working ok for
silicon atoms, even at as high a temperature as 3000K.

system at 3000K for 500k timesteps, i.e. for 50 and 500 ps, respectively,
but the system remained in the crystalline phase when dump file was observed
using Ovito even though the temperature profile was stabilized at 3000K.

now that should give you pause. something is wrong there. you may have
simply created a "flying icecube". that is quite possible, since i
don't see in the script that you posted a velocity command assigning
an initial velocity distribution

However, for tstep=0.01 fs crystal symmetry got broken and atoms were
arranged in random order like amorphous or liquid phase . I did the rdf

that should give you even more pause. if anything, you would see
breaking symmetry by preference at *larger* time steps due to less
accurate numerical integration of the equation of motions.

calculation again for this combination of tstep =0.01 fs and tdamp = 100 but
still getting the same rdf.

so have you validated the potential file that you are using? is it
parameterized for real units?
i just did a quick test with the Si(C) atom type from the
SiCGe.tersoff potential file i mentioned earlier (of course using
metal units) and it reproduces your reference g(r) pretty much spot
on.

a g(r) is not a very sensitive property, so to be far off means that
something is very badly wrong.

I will try to change the density to the liquid density and melt the system
without relaxing for the box size as per Axel's recommendation. May be the
wrong density is changing the nearest neighbors distance.

Thanks for your help, let me know if you find something else. Will post my
findings with fix liquid density.

as aidan already pointed out, the major flaw in your whole effort is
your approach of neglecting to properly learn and understand
fundamentals of MD simulations and as a consequence drawing wildly
wrong conclusions. there is a whole lot of textbook knowledge out
there explaining choices of settings and parameters. what you are
trying to do is a rather simple and straightforward procedure, yet you
have managed to come up with rather strange and sometimes absurd
conclusions and choices.

the procedure is simple:

- if possible, always use an independently tested and validated
potential or force field.
- set up a system and the desired density
- assign initial velocities (making certain, there is no center of
mass momentum) and run MD until the system is reasonably well
equilibrated
- run production for a suitable amount of time and make certain your
data collection is meaningful (not too often, not to infrequent).

this is all stuff that should be well explained in any good text book
on MD. there is no point in worrying about details like thermostat
time constants, if your fundamentals are wrong.

axel.