[lammps-users] alpha_quartz(SiO2) system and replicate problem

Dear all:
I want to simulate a system with alpha_quartz(SiO2) and water in it. I use CHARMM force field for SiO2, so
I have to define bonds, angels etc. in the data file. It’s easy to create a data file for unit cell (triclinic), but when I use replicate command to create a larger
system, I found it hard to define bonds, angles that across the periodic boundary, i.e how does the replicate command deal with bonds, angles in a crystal
starting from a unit cell?
I know that when using CHARMM code, they can define “patches” in the topology file to deal with this,
but I don’t know how can lammps do this. Or should I use a force field that don’t need explicit define of bonds? Then, is there any for alpha-quartz (I don’t find ref)?


There should be no issue with replicating a system with bonds
that straddle periodic boundies. Though it probably hasn't been
widely tested for triclinic systems. However, these lines from
the replicate doc page are very important. If you don't do
this, it won't work:

All properties of the atoms are replicated, including their
velocities, which may or may not be desirable. New atom IDs are
assigned to new atoms, as are molecule IDs. Bonds and other topology
interactions are created between pairs of new atoms as well as between
old and new atoms. This is done by using the image flag for each atom
to "unwrap" it out of the periodic box before replicating it. This
means that molecular bonds you specify in the original data file that
span the periodic box should be between two atoms with image flags
that differ by 1. This will allow them to be unwrapped appropriately.


Actually I noticed this, but I’m a little confused. For example, I have a simulation box (unit cell) like this:

1 2-3 4|

1,2,3,4 are atoms with image flag 0 0 0, “-” indicates bond and atoms 1 and 4 are also bonded, i.e. bond section are like this:


1 1 2 3
2 1 1 4

when I replicate it, it will be like this

1 2-3 4 5 6-7 8|… and with bonds 1-4, 5-8

is that right? But I want bonds 4-5 and 1-8. What should I do by defining bond section and atom’s image flag?


In the original system, before you replicate, you
can define a bond 4-5 for 2 atoms at opposite
sides of the box, say in the x dimension.
The image flags for atoms 4 and 5, which you
also specify in the original data file, need to
differ by 1. See the read_data doc page
info on image flags for details.

If you do this, then the replicate command should
just work.


Hey Xiong,
I haven’t tried the CHARMM SiO2 potential (do you have a reference for it? I’m curious about it), but the BKS potential is a widely used silica potential that works very well for alpha-quartz, and it’s a pair-only potential, which does mean that you don’t have to mess around with moving around angles and bonds (but may, depending on which potential you’re comparing it to, come with some reduced accuracy). The paper is:
B. W. H. van Beest; G. J. Kramer; R. A. van Santen; “Force Fields for Silicas and Aluminophosphates Based on Ab Initio Calculations”, Phys. Rev. Lett. 64 1955-1958 (1990).
There are tons of silica potentials out there, many of which work well for alpha-quartz, so if you’ve only found one so far, you may want to go back and search a bit more thoroughly.

Dear Tyler,
Thanks. My system has quartz and water in it, so if I choose BKS potential for quartz, I have to find another potential for the interaction between quartz and water. Here are the refs. for CHARMM SiO2+water potential:

Lopes, P.E.M., et al., Development of an Empirical Force Field for Silica. Application to the Quartz−Water Interface. The Journal of Physical Chemistry B, 2006. 110(6): p. 2782-2792.

Cruz-Chu, E.R., A. Aksimentiev and K. Schulten, Water−Silica Force Field for Simulating Nanodevices. The Journal of Physical Chemistry B, 2006. 110(43): p. 21497-21508.

It’s a little bit complex, but is specific for quartz+water system.


Right, sorry I missed that key fact. Well, one for alpha-quartz and water that has no defined bonds/angles is the SiO2 reaxff:
Adri C. T. van Duin, Alejandro Strachan, Shannon Stewman, Qingsong Zhang, Xin Xu, and, William A. Goddard, III. ReaxFFSiO Reactive Force Field for Silicon and Silicon Oxide Systems 2003 107 (19), 3803-3811
Thanks for the citations, I’ll check those out.

Dear Xiong and Tyler,

You can also use the simple Lee and Rossky potential, if quartz is maintained as rigid solid. These potential have been well studied for water-silica system.

Here is one of the reference, Malani, Ayapp and Murad, J. Phys. Chem. B, 113, 13825, 2009. ( http://pubs.acs.org/doi/full/10.1021/jp902562v)


One more question: Should I write the parameter file ffield.reax according to the reference manually? It’s prone to have errors.
I found a file in lammps /potentials/ffield.reax.budzien which contains elements Si, O, H, C etc, is it suitable to just use this file?


Looks like that one’s from a different paper:
But it doesn’t look too far off. How about we both try it and compare? I’m about to try reax myself.

Yes, sure. I’m trying simple Lee and Rossky potential.