Lammps users,
I am having some trouble modeling silica via Buckingham. I was able to model it with tersoff just fine but I am trying to compare results of both potentials.
I took the parameters from the BKS paper but when I run the system, my temperature never settles and will eventually explode and atoms are lost. I feel like I have followed the appropriate protocol for determining Tdamp and Pdamp. Does anyone see a red flag in my code that could be causing my system not to settle to desired conditions. I also know that unless my system is sufficiently big, then I will get fluctuations. Even when I replicate the temperature does not seem to settle down.
Any advice would be appreciated.
Ben
units metal
atom_style charge
atom_style charge
kspace_style pppm 1.0e-5
neighbor 1.0 bin
neigh_modify delay 0 every 1 check yes
#read_restart data.restart1
read_data SiO2.data.txt
replicate 2 2 2
pair_style buck/coul/long 10.0 12.0
pair_coeff 1 1 1388.77 .3623188 175.0
pair_coeff 1 2 18003 .2052124 133.5381
pair_coeff 2 2 0 .1 0
thermo_style custom step temp pe ke etotal press vol
thermo_modify norm yes
timestep 0.0001
#Minimize the potential
thermo 1
fix 1 all nve
fix 2 all box/relax iso 1.0 vmax 0.001
minimize 1.0e-12 1.0e-12 1000 10000
unfix 1
velocity all create 600.0 1234567 rot yes mom yes dist gaussian
timestep 0.0002
fix 21 all npt temp 300.0 300.0 0.02 iso 0.0 0.0 0.05
run 100000
Your script has a minor problem of using "fix 1 all nve" and "fix 2
all box/relax iso 1.0 vmax 0.001" together for minimize and again
box/relax with npt for the run. But the effect should be negligible.
Are you sure your structure is good? I can run a 1200
alpha-cristobalite and a 540 atom alpha-quartz without any problem.
Ray
I got the structure from VMD inorganic builder and used topotools to convert to LAMMPS.
How quickly does your structure converge to 300K with my script?
Ben
Within the first 1000 steps. With 1200 atoms, P oscillates about
plus/minus 3000 bars.
Ray
That’s interesting. I am not sure where the problem is. I suppose it could be my data file but I doubt VMD’s inorganic builder for silica would be wrong. Somehow my implementation of it has gone astray. I will continue experimenting to see if I can find the bug.
Ben
Ray,
Are you getting a potential when normalized of about -19.3eV?
Ben
That's interesting. I am not sure where the problem is. I suppose it could
be my data file but I doubt VMD's inorganic builder for silica would be
wrong. Somehow my implementation of it has gone astray. I will continue
experimenting to see if I can find the bug.
what about bonded interactions?
remember if you add bonds, it will lead to exclusions, even if you
don't define a bond style.
axel
I have not added any bonds. Everything is modeled simply with the buck/coul/long potential. I think it may be working now. I used vmd and replicated the alpha quartz structure in LAMMPS. It equilibrated at 300K (with fluctuations due to the small system) after about 1000-1500 timesteps. It also has a normalized potential energy of about -19.3eV/atom which closely matches Paul Croizer’s -19.1eV from his tabulated BKS potential.
Ben
With your script, I get -19.44 for quartz and -19.37 for cristobalite.
Ray
Interesting, I am getting -19.37 for quartz (VMD’s alpha quartz).
Ben