[lammps-users] Problem with pair-style table and BKS potential for a-Quartz silica

Dear Users,

I am having a problem when trying to implement the BKS potential for a-quartz using the Pair-style table command.
My code used to run fine when implementing the BKS potential using the Buck/coul/long pair style.

Now I want to add an additional 24-6 LJ potential on top of the BKS potential so I want to implement these two using a table (No built in 24-6 LJ potential in lammps so cannot use Hybrid overlay command).

The problem is when I use the Pair style table, I lose atoms:

  1. When not specifying RSQ parameters in the table I lose atoms after about 160 steps
  2. When specifying RSQ, I lose almost all atoms at the first step.

Temperature and kinetic energy start increasing from the initial 450K (although I run NVT 300K) to 10^12 K before losing atoms.

Pressure is negative at the first step.

  1. When not specifying RSQ in table, pressure is --9035954.5 (metal units) and steady till the last three steps before losing atoms where it explodes up to +10^13 (metal units)
  2. When specifying RSQ, I lose atoms after the first step and pressure is -10^11 (metal units).

I used the same parameters to generate my potential and force table in excel as in the Buck/coul/long potential. The same ones that from the original Ban best/Kramer/Van Santen (1990) publication. And I plotted them to double check they looked all good.

Also I tried to refine my table more and more (I am currently using 10000 data point for each pair, ranging from 0.0001A to 10A plus Kspace PPPM 10^-4.
I plotted the potentials using Pair write and they look very good. No problem.

Once again, This code used to be working good and giving good thermal conductivity values when using Pair style buck/coul/long. Here the only thing I did was to replace it by a table.

Here is the input script followed by the first few lines of the Potential table file.

Does anybody know what it is that I am doing wrong. Help would be greatly appreciated.

Thank you very much in advance.

Sincerely,

INPUT SCRIPT

#Initialization############################################################

units metal

dimension 3

boundary p p p

atom_style charge

Atom definition##########################################################

lattice custom 5.4054 &

a1 0.9095 0.0000 0.0000 &

a2 -0.4547 0.7876 0.0000 &

a3 0.0000 0.0000 1.0000 &

basis 0.4697 0.0000 0.0000 basis 0.0000 0.4697 0.6667 basis 0.5303 0.5303 0.3333 basis 0.4135 0.2669 0.1191 basis 0.2669 0.4135 0.5475 basis 0.7331 0.1466 0.7857 basis 0.5865 0.8534 0.2142 basis 0.8534 0.5865 0.4524 basis 0.1466 0.7331 0.8809

region simbox block 0 4 0 4 0 4 units lattice

create_box 2 simbox

create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 &

basis 4 2 basis 5 2 basis 6 2 &

basis 7 2 basis 8 2 basis 9 2

mass 1 28.0855

mass 2 15.9994

group siliconatoms type 1

group oxygenatoms type 2

set group siliconatoms charge 2.4

set group oxygenatoms charge -1.2

Atoms interactions settings##################################

Si type 1, O type 2

#pair_style buck/coul/long 8 8

#pair_coeff 1 1 0.0 0.10 0.0

#pair_coeff 1 2 18003.7572 0.205205 133.5381

#pair_coeff 2 2 1388.7730 0.362319 175.0

pair_style table linear 10000

pair_coeff 1 1 BKSLJ.table SiSi

pair_coeff 1 2 BKSLJ.table SiO

pair_coeff 2 2 BKSLJ.table OO

kspace_style pppm 1.0e-4

neighbor 2.0 bin

neigh_modify every 1 delay 0 check yes

pair_write 1 2 1000 r 0.01 5.0 table.txt BKS_Si_O 2.4 -1.2

pair_write 2 2 1000 r 0.01 5.0 table.txt BKS_O_O -1.2 -1.2

pair_write 1 1 1000 r 0.01 5.0 table.txt BKS_Si_Si 2.4 2.4

Equilibrate temperature in NVT

fix 1 all nvt 300.0 300.0 100.0

velocity all create 450.0 5674845

timestep 0.00055

thermo_style custom step temp press &

ke pe etotal vol lx ly lz epair

thermo 1

dump 1 all atom 1000 dump1.qua

run 100000

NOW HERE IS AN EXAMPLE OF THE POTENTIAL TABLE FILE

table for BKS only Potential for LAMMPS

Si-O potential first, O-O potential second, SiSI third

4 columns with: index, distance, Potential, Force

SiO
N 10000 RSQ 0.0001 9.9991

1 0.0001 -1.33538E+26 -8.01229E+30
2 0.0011 -7.53788E+19 -4.11157E+23
3 0.0021 -1.557E+18 -4.44858E+21
4 0.0031 -1.50465E+17 -2.91222E+20
5 0.0041 -2.81127E+16 -4.11405E+19
6 0.0051 -7.58899E+15 -8.92823E+1

9997 9.9961 -4.148849365 -0.415113756
9998 9.9971 -4.148434293 -0.415030674
9999 9.9981 -4.148019304 -0.414947616
10000 9.9991 -4.147604398 -0.414864583

OO
N 10000 RSQ 0.0001 9.9991

1 0.0001 -1.75E+26 -1.05E+31
2 0.0011 -9.87829E+19 -5.38816E+23
3 0.0021 -2.04043E+18 -5.82981E+21

9998 9.9971 2.073974958 0.207369985
9999 9.9981 2.073767609 0.207328557
10000 9.9991 2.073560301 0.207287142

SiSi
N 10000 RSQ 0.0001 9.9991

1 0.0001 0 0
2 0.0011 0 0
3 0.0021 0 0
4 0.0031 0 0

Dear Users,

dear thomas,

Now I want to add an additional 24-6 LJ potential on top of the BKS
potential so I want to implement these two using a table (No built in
24-6 LJ potential in lammps so cannot use Hybrid overlay command).

how about using a table only for the 24-6 LJ
(quickly) writing a 24-6 LJ style.

since there are now a large number of different
exponent LJ styles around, perhaps this would be
the time to implement a generic m-n-LJ style that
would make most of the others obsolete?

The problem is when I use the Pair style table, I lose atoms:
1. When not specifying RSQ parameters in the table I lose atoms after
about 160 steps
2. When specifying RSQ, I lose almost all atoms at the first step.

hmmm... 24-6 is a very steep potential on the repulsive side.
have you tried reducing the time step to handle this better?

second thing. are you sure about the sign of your table?
i ran into this myself recently, and it took a while to
realize that the prefactor in the formula had a sign that
was opposite from what i expected.

cheers,
   axel.

Dear Thomas, Axel and all,

Just to recall the context of this issue: It is necessary to add a
short range strong repulsion to the BKS model of silica, in order to
avoid unphysical fusion events between Si and O at high temperatures.
There has been approximately one way to do it per paper published on
the subject.

For instance, following a paper on the phase diagram of the BKS
potential (http://scholar.google.fr/scholar?q=Phase+diagram+of+silica+from+computer+simulation),
I recently wrote a 30-6 LJ potential. Please find enclosed the source
files, and the corresponding parameters for lammps (in metal units):

pair_style hybrid/overlay lj306/cut 10.0 buck/coul/long 10.0 8.0
pair_coeff * * none
pair_coeff 8 8 buck/coul/long 1389.0 0.3623 175.0
pair_coeff 8 14 buck/coul/long 18000.0 0.2052 133.5
pair_coeff 14 14 buck/coul/long 0.0 0.3333 0.0
pair_coeff 8 8 lj306/cut 0.001051 1.779
pair_coeff 8 14 lj306/cut 0.003098 1.314
pair_coeff 14 14 lj306/cut 0.0 1.5

As for implementing a generic m-n LJ potential, I would be willing to
do it if you let me some time. Yet there is the problem that it won't
be possible to make specific optimizations to compute the powers of r
as it is done in the standard and 9-6 LJ implementations...

Best regards,
Laurent

pair_lj306_cut.cpp (20.1 KB)

pair_lj306_cut.h (1.52 KB)

Dear Axel, Dear All,

Thank you for your answer,

Now I want to add an additional 24-6 LJ potential on top of the BKS
potential so I want to implement these two using a table (No built in
24-6 LJ potential in lammps so cannot use Hybrid overlay command).

how about using a table only for the 24-6 LJ
(quickly) writing a 24-6 LJ style.

I wanted to implement the BKS potential myself too because I am worried the one obtained using buck/coul/long and the BKS parameters (BKS 1990) is not exactly the same because of maybe a difference in the hidden factor C in front of the coulombic term (“energy conversion factor”).

since there are now a large number of different
exponent LJ styles around, perhaps this would be
the time to implement a generic m-n-LJ style that
would make most of the others obsolete?

I am running LAMMPS on some Cluster where I cannot recompile it so I cannot really build my own new pair_style but yes having this would be very useful.

The problem is when I use the Pair style table, I lose atoms:

  1. When not specifying RSQ parameters in the table I lose atoms after
    about 160 steps
  2. When specifying RSQ, I lose almost all atoms at the first step.

hmmm… 24-6 is a very steep potential on the repulsive side.
have you tried reducing the time step to handle this better?

Sorry I was maybe not clear enough on this but first I tried with the BKS+24-6LJ potential in my table. Then because I was having this problem of losing atoms, I removed the 24-6LJ potential to just keep the BKS potential in my table. And still the exact same problem remained so it is not linked to the LJ potential but to an error I must be making when implementing the pair-style table.

second thing. are you sure about the sign of your table?
i ran into this myself recently, and it took a while to
realize that the prefactor in the formula had a sign that
was opposite from what i expected.

When I plot the potential and force using Pair-Write… I get totally relevant curves.

Subsidiary question: Do I still need to keep defining each atoms charges in the input script? Since basically all the information is already contained into the table?

I have been stuck on this for multiple days. Any further insight would be greatly appreciated.

Thank you very much in advance.

Sincerely,

Thomas Coquil
Research Assistant
UCLA Mechanical and Aerospace Engineering Department
420 Westwood Plaza
43-132 ENGINEERING IV
Los Angeles CA, 90095
Tel: 1 626 429 2110
[email protected]…24…

Dear Laurent, Axel and all,

Thank you for your help and for providing your own source files. I appreciate.

I am aware of this issue with the BKS potential at high temperature but I am simply trying to run it at 300K here and I do not know why temperature starts increasing like crazy at some point.

As I was trying to say in my previous email, I tried first with a 24-6LJ potential on top of the BKS potential. Then I tried with the BKS potential alone. Both time using the pair_style table to implement them. And I have been having the same problem both times… at 300K with nvt or nve.

Any ideas where it could come from?

Thank you so much

Sincerely,

Thomas Coquil
Research Assistant
UCLA Mechanical and Aerospace Engineering Department
420 Westwood Plaza
43-132 ENGINEERING IV
Los Angeles CA, 90095
Tel: 1 626 429 2110
[email protected]…24…

As for implementing a generic m-n LJ potential, I would be willing to
do it if you let me some time. Yet there is the problem that it won't
be possible to make specific optimizations to compute the powers of r
as it is done in the standard and 9-6 LJ implementations...

most of these optimizations could be done at compile time using
templated functions. in the simplest case one you just program an
if/elseif/elseif/elseif/ block each time specific combinations of exponents
are required and then the compiler would eliminate the unnecessary
branches for the give M-N exponent combination.

or even simpler is to do it with a switch statement like our cg/cmm
potential does, which allows to select 12-4,9-6, or 12-6 depending
on the given pair. extending that for other exponent combinations
would be fairly simple.

cheers,
     axel.

I wanted to implement the BKS potential myself too because I am worried the
one obtained using buck/coul/long and the BKS parameters (BKS 1990) is not
exactly the same because of maybe a difference in the hidden factor C in
front of the coulombic term ("energy conversion factor").

lammps has that included. check out the use of force->qqrd2e
in all coul containing pair styles.

[...]

I am running LAMMPS on some Cluster where I cannot recompile it so I cannot
really build my own new pair_style but yes having this would be very useful

who is running a cluster where you cannot compile your
own application? that is like owning a ferrari in an area
with only dirt roads....

[...]

Subsidiary question: Do I still need to keep defining each atoms charges in
the input script? Since basically all the information is already contained
into the table?

charges are ignored if you use pair style table.

I have been stuck on this for multiple days. Any further insight would be
greatly appreciated.

the best way to get help is to provide the means to exactly
reproduce what is giving your the problems.

cheers,
    axel.

Dear Laurent, Axel and all,
Thank you for your help and for providing your own source files. I
appreciate.
I am aware of this issue with the BKS potential at high temperature but I am
simply trying to run it at 300K here and I do not know why temperature
starts increasing like crazy at some point.
As I was trying to say in my previous email, I tried first with a 24-6LJ
potential on top of the BKS potential. Then I tried with the BKS potential
alone. Both time using the pair_style table to implement them. And I have
been having the same problem both times... at 300K with nvt or nve.
Any ideas where it could come from?

so your tables are not good enough for how you run your input. there are
a lot of possible reasons, but without exactly seeing what you are doing,
it is near impossible to tell.

cheers,
     axel.

As for implementing a generic m-n LJ potential, I would be willing to
do it if you let me some time. Yet there is the problem that it won't
be possible to make specific optimizations to compute the powers of r
as it is done in the standard and 9-6 LJ implementations...

most of these optimizations could be done at compile time using
templated functions. in the simplest case one you just program an
if/elseif/elseif/elseif/ block each time specific combinations of exponents
are required and then the compiler would eliminate the unnecessary
branches for the give M-N exponent combination.

or even simpler is to do it with a switch statement like our cg/cmm
potential does, which allows to select 12-4,9-6, or 12-6 depending
on the given pair. extending that for other exponent combinations
would be fairly simple.

Alas this kind of things goes well beyond my programing skills, I'm
afraid I'll have to let someone else do it...

Best regards,
Laurent

Dear Axel and all,

the best way to get help is to provide the means to exactly
reproduce what is giving your the problems.

so your tables are not good enough for how you run your input. there are
a lot of possible reasons, but without exactly seeing what you are doing,
it is near impossible to tell.

Thank you very much, I am providing here my input script and my potential/force table for your convenience. It’s a very simple script so the problem must just be a little detail I am not getting right.

charges are ignored if you use pair style table.

I just realized I might actually be making a mistake here. I do not have much experience with LAMMPS indeed. But since I am now using a table instead of the buck/coul/long pair style… Do I need to add an additional coul/long pair style with Hybrid command so that long range attractive coulombic interactions can be accounted for using the Kspace solver?

lammps has that included. check out the use of force->qqrd2e
in all coul containing pair styles.

I agree, but still something must be a little different because when plotting the potentials from LAMMPS buck/coul/long and using the original parameters from scratch… the two curves have an offset. Force curves (derivatives) look very similar however so it’s fine.

I am running LAMMPS on some Cluster where I cannot recompile it so I cannot
really build my own new pair_style but yes having this would be very useful

who is running a cluster where you cannot compile your
own application? that is like owning a ferrari in an area
with only dirt roads…

I agree, I just thought it would be simpler to try by myself with a table first though.

Thank you so much for all this help.

Best regards,

Thomas

in.quaLJ2 (2.28 KB)

BKSLJ.table (874 KB)

Dear Axel and all,

dear thomas,

Thank you very much, I am providing here my input script and my
potential/force table for your convenience. It's a very simple script
so the problem must just be a little detail I am not getting right.

the BKSLJ.table does not include any charges.
also it is of R style but advertises RSQ.

if i use:

pair_write 1 1 1000 rsq 0.2 10.0 table.txt SiSi 2.4 2.4
pair_write 1 2 1000 rsq 0.2 10.0 table.txt SiO 2.4 -1.2
pair_write 2 2 1000 rsq 0.2 10.0 table.txt OO -1.2 -1.2

with the original non-tabulated pair style, i can generate
a table that reproduces the analytical result fairly close
using:

pair_style table spline 1000
pair_coeff 1 1 table.txt SiSi 8.0
pair_coeff 1 2 table.txt SiO 8.0
pair_coeff 2 2 table.txt OO 8.0

original:
Step Temp Press KinEng PotEng TotEng Volume Lx Ly Lz E_pair
       0 450 -23524.037 50.198177 -16762.602
-16712.404 10860.435 29.496187 17.029172 21.6216
-16762.602
      20 315.05909 10329.846 35.145315 -16747.719
-16712.574 10860.435 29.496187 17.029172 21.6216
-16747.719
      40 304.88494 61856.368 34.010374 -16746.659
-16712.649 10860.435 29.496187 17.029172 21.6216
-16746.659
      60 281.11586 71712.592 31.358897 -16744.002
-16712.644 10860.435 29.496187 17.029172 21.6216
-16744.002

table:
Step Temp Press KinEng PotEng TotEng Volume Lx Ly Lz E_pair
       0 450 -23524.905 50.198177 -16762.628
-16712.429 10860.435 29.496187 17.029172 21.6216
-16762.628
      20 315.06005 10328.972 35.145423 -16747.738
-16712.593 10860.435 29.496187 17.029172 21.6216
-16747.738
      40 304.88546 61855.418 34.010432 -16746.677
-16712.666 10860.435 29.496187 17.029172 21.6216
-16746.677
      60 281.11618 71711.975 31.358933 -16744.019
-16712.66 10860.435 29.496187 17.029172 21.6216
-16744.019

charges are ignored if you use pair style table.

I just realized I might actually be making a mistake here. I do not
have much experience with LAMMPS indeed. But since I am now using a
table instead of the buck/coul/long pair style.... Do I need to add an
additional coul/long pair style with Hybrid command so that long range
attractive coulombic interactions can be accounted for using the
Kspace solver?

you only need to add a coul/long via hybrid overlay,
if you don't include the charges _with_ the long-range
dependent damping into the table. the above table output
includes the charges in the tables accordingly.

for the sake of simplicity, it is probably easiest to go
that route, as this would allow you to produce a proper
table manually without having to mess with erfc().

cheers,
    axel.

Dear Axel, dear all,

Just one more question please:

You said yesterday that I had two options:

  1. Generating the table from the original non tabulated Buck/coul/long potential (like you did)
  2. Generating the table myself without the coulombic term and then add a coul/long term on top of it with hybrid/overlay to use the pppm kspace solver.
  • You said the first route would easier because not having to mess with and erfc-style screening function…
    What I do not get is… If using the first method with the coulombic term already in the table… How are the long range coulombic interaction accounted for when using the table as the potential…? Since this table has a finite (pretty short) cutoff…! right?
    The kspace solver does not work with the table anymore here… since charges are already included into the table?

  • If I use method 2: What is the error I am making without using and erfc function?
    I found your previous thread on this: LAMMPS mail list thread #1996:

please note that the point of ewald summation is that you screen
out the long range coulomb interactions in real space and the 
short range coulomb in k-space and thus the screening has to have
the exact same shape in both cases (using an erfc-style screening 
function equivalent to adding/subtracting a gaussian).

But I do not have much experience with this and I do not really see how to implement this and make the screenings have same shape…

Thank you very much for your insight

Best regards,

Thomas Coquil
Research Assistant
UCLA Mechanical and Aerospace Engineering Department
420 Westwood Plaza
43-132 ENGINEERING IV
Los Angeles CA, 90095
Tel: 1 626 429 2110
[email protected]…24…

Dear Axel, dear all,
Just one more question please:
You said yesterday that I had two options:
1. Generating the table from the original non tabulated Buck/coul/long
potential (like you did)

that is not what i meant. i meant generating the table from scratch with all
short range terms included. i.e. the buckingham, the lennard-jones and the
screened coulomb. for proper long-range electrostatics, you would have
to use pppm with that and in that case you have to make sure that you
set the proper real space cutoff that matches what is used for the screening.

2. Generating the table myself without the coulombic term and then add a
coul/long term on top of it with hybrid/overlay to use the pppm kspace
solver.
- You said the first route would easier because not having to mess with and
erfc-style screening function...

there is a third option: you could tabulate the 24-6 lennard jones and
overlay this with the existing buck/coul/long pair style.

What I do not get is... If using the first method with the coulombic term
already in the table... How are the long range coulombic interaction
accounted for when using the table as the potential...? Since this table has
a finite (pretty short) cutoff....! right?

exactly. you can still use the kspace solver, but you have to match the
cutoffs with what is used in the screening to indicate to the pppm code
whay settings to use. this is what makes it a bit messy. as you saw
from the example i quoted, it does work. only my example was missing
the 24-6 lj term that you wanted to add in the first place.

The kspace solver does not work with the table anymore here... since charges
are already included into the table?

kspace works with tables, if you provide a cutoff (and the same cutoff for
all tables).

- If I use method 2: What is the error I am making without using and erfc
function?

in the second case, the coulomb is properly treated with the coul long.

you can test this by using hybrid/overlay on plain pair_style buck
with coul/long and compare to buck/coul/long. for your tables, you'd
then just replace the analytical buck with buckingham plus 24-6 lj.

I found your previous thread on this: LAMMPS mail list thread #1996:

please note that the point of ewald summation is that you screen
out the long range coulomb interactions in real space and the
short range coulomb in k-space and thus the screening has to have
the exact same shape in both cases (using an erfc-style screening
function equivalent to adding/subtracting a gaussian).

But I do not have much experience with this and I do not really see how to
implement this and make the screenings have same shape....

that is why i suggested to use coul/long with hybrid overlay (as counterpart
to pppm) to treat the short range part of coulomb and putting the
non-coulomb part together in one table.

cheers,
     axel.

Paul - don't know if you've seen this thread in the LAMMPS mail
list, but someone is having problems using pair table for short-range
BKS in silica - a problem you solved before. Maybe you
want to check the messages and weigh in.

Steve

Thomas,

As Axel mentioned, you can use PPPM with pair-style table, but it is a bit tricky. When you run the simulation with pair-style table, you have to make sure you are using the same g_ewald parameter that was used when the tables were created, presumably using a pair_write command. This can be done with a command like:

kspace_modify gewald 0.29

There are also some other little subtleties when using pair-style table that you need to be careful about in order to get it right. See:

http://lammps.sandia.gov/doc/pair_table.html

As Steve mentioned, I've solved this problem before, similarly to what you're doing, but with a different steep potential to avoid the "fusion" problem. I'll attach some input from my old BKS quartz work that might offer some clues.

Paul

quench.in (664 Bytes)

silica.tabulated (121 KB)

Dear Paul, dear Steve,

Thank you very much for your help. I appreciate. I think I found a way to do it that I like using a table parametrizing the Buckingam part of the BKS (using Beest/Kramed/Santen parameters) + my 24-6 LJ and by adding the coul/long term using hybrid overlay and the pppm solver.
That should do it.

I am indeed looking forward to quenching a-quartz to make amorphous silica.

Paul, is that the script (that you forwarded) you actually used for quenching? I have never performed a quench before. Are there other things I should be careful with when doing it?

Thank you so much.

Sincerely
Thomas Coquil
Research Assistant
UCLA Mechanical and Aerospace Engineering Department
420 Westwood Plaza
43-132 ENGINEERING IV
Los Angeles CA, 90095
Tel: 1 626 429 2110
[email protected]…24…

Thomas,

You’re welcome. I hope you are successful in getting this going.

The script I sent was used for quenching a-quartz to make amorphous silica, so if you follow the procedure, it should work the same for you. I’m also attaching the corresponding data file. If you update the input script I sent to work with the latest version of LAMMPS, you should be able to run the same simulation I ran and get a feel for the procedure I used.

Best wishes,

Paul

quench.data (109 KB)