LAMMPS users,
I have been trying to reproduce lattice parameter calculations in the literature. If I start with a unit cell and use lammps replicate command, wouldn’t I just divide lx, ly, lz by the replication factor?
For instance, if I have
replicate 7 7 7
variable length equal “lx/7”
variable length2 equal “ly/7”
variable length3 equal “lz/7”
should be right shouldn’t it? I have noticed that as I increase replicate and divide by replication factor, my lattiice constants are actually getting much smaltler.
I want to reproduce BKS lattice parameters. I can get very close but the paper does not specify a cutoff. How can I know which cutoff to use? If I use cutoff of 5 I get close to desired result. If I use 10 my lattice constant calculations decrease by around 30%.
Can anyone explain why I am seeing this?
Ben
LAMMPS users,
I have been trying to reproduce lattice parameter calculations in the
literature. If I start with a unit cell and use lammps replicate command,
wouldn't I just divide lx, ly, lz by the replication factor?
For instance, if I have
replicate 7 7 7
variable length equal "lx/7"
variable length2 equal "ly/7"
variable length3 equal "lz/7"
should be right shouldn't it? I have noticed that as I increase replicate
and divide by replication factor, my lattiice constants are actually getting
much smaltler.
I want to reproduce BKS lattice parameters. I can get very close but the
paper does not specify a cutoff. How can I know which cutoff to use? If I
use cutoff of 5 I get close to desired result. If I use 10 my lattice
constant calculations decrease by around 30%.
Can anyone explain why I am seeing this?
it is always dangerous to speculate with extremely limited information.
if the behavior changes with the cutoff, your cutoff is too small.
if you get the right result with the wrong cutoff, your potential
parameters are not correct.
you didn't mention how you handle long-range electrostatics.
axel.
Axel,
for long range electrostatics I am using buck/coul/long.
I doubt my parameters are wrong as I am taking them directly from BKS paper. However, I put them below anyways.
The partial charges are +2.4 for Si and -1.2 for O as in the paper.
I am using:
pair_style buck/coul/long 10 10
pair_coeff 1 1 1388.77 .3623188 175.0
pair_coeff 1 2 18003 .2052124 133.5381
pair_coeff 2 2 0 .1 0
After minimizing, I calculate lattice constants. Since I am using VMD alpha-quartz, I should expect to see around 4.913 for a and 5.4 for c but I am getting values in the 2 and 3 range. With low cutoff, I get close to these values but with higher cutoff, I do not.
Another question - What should I use as a cutoff if I am just simulating the unit cell? If I specify 10 angstroms, will the interactions be calculated through periodic boundary conditions multiple times? This is not my main question in this thread - just a side note. I am replicating the unit cell before minimization (I know I could do it after) and then dividing lx, ly, lz by replication factor.
Any help would be appreciated. Thank you for the help Axel.
Ben
Just a side thought… What are your units… metal or real? Are your parameters correct with right units?
My units are metal.
Metal units correspond to the BKS paper as the potential is in the form of eV and Angstroms. And I converted the b term in the paper to rho in LAMMPS buckingham form since they are inverses of one another.
Ben
Are you sure you are using a long range solver for Coulombic interactions?
The default kspace_style is none. Additionally, are you sure your pair coefficient types are properly assign, i.e. which atom is type 1 and type 2?
Carlos
Carlos, I am using pppm to do long-range Coulombics. Additionally, my oxygen is assigned as 1 and Silicon as 2.
I have included my code. I figured that would easier to show exactly what I have done.
units metal
atom_style charge
kspace_style pppm 1.0e-5
kspace_modify gewald 0.29
neighbor 1.0 bin
neigh_modify delay 0 every 5 check yes
boundary p p p
#read_restart data.restart1
read_data SiO2.data
replicate 20 20 20
pair_style buck/coul/long 10 10
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 lx ly lz
thermo_modify norm yes
timestep 0.001
#Minimize the potential
thermo 1
fix 1 all nve
fix 2 all box/relax iso 1.0 vmax 0.001
minimize 1.0e-25 1.0e-25 1000 10000
unfix 1
variable length equal “lx/20”
variable length2 equal “ly/20”
variable length3 equal “lz/20”
print “Lattice constant (Angstoms) = {length};"
print "Lattice constant (Angstoms) = {length2};”
print “Lattice constant (Angstoms) = ${length3};”
Thank you for all the help.
Ben
Benjamin,
to start with, you dont need timestep and fix nve for minimization. Second, your number of steps of minimization is quite low. This is what I suggest. You print out your length, length2, length3 in thermo line (define the variables before thermo). and see how those are changing as you are minimizing your system.
Best regards,
Vikas
Make sure as well that the quotients lx,ly,lz/20 are indeed the lattice constants of the primitive cell. Computing the volume/atom in the primitive cell can help you figure out if you are analyzing things properly.
Carlos
Vikas,
After taking your your suggestions and running it, I get a = 4.8 and c - 5.29.
I know these values are very close to published results - I am still uncertain as to why they are not closer. The BKS paper has 4.9 or so for a and 5.4 for c.
Ben
Did you check if the published results are actually MD at certain temperature and not at 0 K? What I see is that you have lower values. May be if you run MD simulation at room temperature and equilibrate it, the cell parameters will increase and will lead to better agreement.
Vikas,
I am guessing you suggesting a fix npt? Fix nve or nvt will not work because the lattice dimensions wouldnt change due to constant volume. I am running it but it appears that c will be around 5.25 so I do not think it is increasing it.
Ben
Yes, I was suggesting fix npt with independent barostats along x y and z direction.
I am surprised, it is lower than your minimized structure. In generally lattice constants are expected to increase with temperature/ You have to wait until pressure fluctuations become steady. It might be a while. I am not sure what temperature you are running your simulation. did you visualize your system to check it is stable?
Vikas,
Why cant I simply use:
timestep 0.0005
fix 2 all npt temp 300 300 0.05 iso 1 1 0.5
Would that not set a barostat at 1 bar? What is the advantage of setting independent barostats?
Yes, it is surprising - I was assuming it would get larger as well.
Ben
Dear Benjamin,
I think I know the issue.
The unit cell parameters that you are providing us are of trigonal system (one of the angle is 120 degrees). The datafile that you are running your job is orthorhombic. It is messing up the box dimensions. Make your box non-orthagonal properly and rerun your simulations. You should be able to get everything right.
Regards,
Vikas
I will make your life simpler. Use this orthogonal box as your unit cell. This is not the primitive unit cell but it is orthogaonal (2X2 of the primitive and transformed to orthogonal). Use updated box lengths.
REMARK Test
REMARK Created: Tue Oct 08 13:48:06 Eastern Standard Time 2013
CRYST1 8.509 9.826 5.405 90.00 90.00 90.00 P1
ORIGX1 1.000000 0.000000 0.000000 0.00000
ORIGX2 0.000000 1.000000 0.000000 0.00000
ORIGX3 0.000000 0.000000 1.000000 0.00000
SCALE1 0.117516 0.000000 0.000000 0.00000
SCALE2 0.000000 0.101771 0.000000 0.00000
SCALE3 0.000000 0.000000 0.185007 0.00000
ATOM 1 OX1 MOL 2 1.760 0.299 1.157 1.00 0.00 O
ATOM 2 OX1 MOL 2 3.116 8.744 2.959 1.00 0.00 O
ATOM 3 OX1 MOL 2 3.634 0.783 4.761 1.00 0.00 O
ATOM 4 OX1 MOL 2 1.139 1.375 4.248 1.00 0.00 O
ATOM 5 OX1 MOL 2 0.621 3.240 2.446 1.00 0.00 O
ATOM 6 OX1 MOL 2 2.495 2.755 0.644 1.00 0.00 O
ATOM 7 SI1 MOL 2 2.000 8.671 1.802 1.00 0.00 Si
ATOM 8 SI1 MOL 2 0.000 2.310 3.603 1.00 0.00 Si
ATOM 9 SI1 MOL 2 2.255 1.302 0.000 1.00 0.00 Si
ATOM 10 OX1 MOL 2 6.015 7.668 1.157 1.00 0.00 O
ATOM 11 OX1 MOL 2 7.371 6.288 2.959 1.00 0.00 O
ATOM 12 OX1 MOL 2 7.888 8.153 4.761 1.00 0.00 O
ATOM 13 OX1 MOL 2 5.393 8.744 4.248 1.00 0.00 O
ATOM 14 OX1 MOL 2 4.876 0.783 2.446 1.00 0.00 O
ATOM 15 OX1 MOL 2 6.750 0.299 0.644 1.00 0.00 O
ATOM 16 SI1 MOL 2 6.255 6.215 1.802 1.00 0.00 Si
ATOM 17 SI1 MOL 2 4.255 9.679 3.603 1.00 0.00 Si
ATOM 18 SI1 MOL 2 6.509 8.671 0.000 1.00 0.00 Si
ATOM 19 OX1 MOL 2 1.760 5.212 1.157 1.00 0.00 O
ATOM 20 OX1 MOL 2 3.116 3.831 2.959 1.00 0.00 O
ATOM 21 OX1 MOL 2 3.634 5.696 4.761 1.00 0.00 O
ATOM 22 OX1 MOL 2 1.139 6.288 4.248 1.00 0.00 O
ATOM 23 OX1 MOL 2 0.621 8.153 2.446 1.00 0.00 O
ATOM 24 OX1 MOL 2 2.495 7.668 0.644 1.00 0.00 O
ATOM 25 SI1 MOL 2 2.000 3.758 1.802 1.00 0.00 Si
ATOM 26 SI1 MOL 2 0.000 7.223 3.603 1.00 0.00 Si
ATOM 27 SI1 MOL 2 2.255 6.215 0.000 1.00 0.00 Si
ATOM 28 OX1 MOL 2 6.015 2.755 1.157 1.00 0.00 O
ATOM 29 OX1 MOL 2 7.371 1.375 2.959 1.00 0.00 O
ATOM 30 OX1 MOL 2 7.888 3.240 4.761 1.00 0.00 O
ATOM 31 OX1 MOL 2 5.393 3.831 4.248 1.00 0.00 O
ATOM 32 OX1 MOL 2 4.876 5.696 2.446 1.00 0.00 O
ATOM 33 OX1 MOL 2 6.750 5.212 0.644 1.00 0.00 O
ATOM 34 SI1 MOL 2 6.255 1.302 1.802 1.00 0.00 Si
ATOM 35 SI1 MOL 2 4.255 4.766 3.603 1.00 0.00 Si
ATOM 36 SI1 MOL 2 6.509 3.758 0.000 1.00 0.00 Si
Sounds like Vikas found the source of the problem. Relieved on this side to know that Van Beest , Kramer and van Santen indeed knew what they were doing (IRONY).
Carlos