Adding Malavasi repulsive core to BKS

Lammps users,

I am trying to implement Malavasi’s potential where a repulsive core is added to BKS potential, as in the paper: http://www.sciencedirect.com/science/article/pii/S0022309305007994

When I model 5,184 atoms of quartz at 300K with just BKS potential, system seems stable and silica is modeled well.

When I add the repulsive core, D/r^12 as given in the paper, the temperature of the system increases to ~5000K and slowly goes down. However, the density also increases whereas, one would think that a repulsive core would cause the density to decrease.

My question is this: What would a soaring temperature in an NPT ensemble that is supposed to stay at 300 K indicate? The system is getting about twice as dense as it should, and like I said, I would expect it to become less dense, if anything.

But at STP, the repulsive core should not come into play anyways.

I know that Dr. Crozier developed a repulsive core tabulated file, but I’m trying to implement the one in the paper. Please give me suggestions on what could be causing this issue.

The general setup of the potential file is such:

O-Si
N 15000 R 0.00100000 15.00000000

1 0.001 1.8E+37 2.16E+41
2 0.002 4.39453E+33 2.63672E+37
3 0.003 3.38702E+31 1.35481E+35
4 0.004 1.07288E+30 3.21865E+33
5 0.005 7.3728E+28 1.76947E+32
6 0.006 8.26909E+27 1.65382E+31
7 0.007 1.30046E+27 2.22936E+30
8 0.008 2.61934E+26 3.92902E+29
9 0.009 6.37327E+25 8.49769E+2

Lammps users,

I am trying to implement Malavasi's potential where a repulsive core is
added to BKS potential, as in the paper:
http://www.sciencedirect.com/science/article/pii/S0022309305007994

When I model 5,184 atoms of quartz at 300K with just BKS potential, system
seems stable and silica is modeled well.

When I add the repulsive core, D/r^12 as given in the paper, the temperature
of the system increases to ~5000K and slowly goes down. However, the density
also increases whereas, one would think that a repulsive core would cause
the density to decrease.

are you sure that you have the correct sign on the force resulting
from that term?

My question is this: What would a soaring temperature in an NPT ensemble
that is supposed to stay at 300 K indicate? The system is getting about

there is not much detail in your description, but in general, any time
to significantly change a potential, your potential energy would go up
and thus you should see a spike in the kinetic energy (that would then
be removed when using a thermostat). this is just common sense, since
you are kicking a system that is in equilibrium out of it. so you will
need to minimize and re-equilibrate.

twice as dense as it should, and like I said, I would expect it to become
less dense, if anything.

But at STP, the repulsive core should not come into play anyways.

I know that Dr. Crozier developed a repulsive core tabulated file, but I'm
trying to implement the one in the paper. Please give me suggestions on what
could be causing this issue.

if the density is too high, then the potential is too attractive. no
big mystery there.

axel.

Dr. Kohlmeyer,

Thank you for your quick response.

I am confident that the sign is correct. I took the negative derivative of potential to find the force. Also, the term is a repulsive term, not an attractive term. When I plot the potential, there are no attractive characteristics about it. It’s just a steep potential that goes to INF at small r. I will continue to search minimization techniques, because the default just goes to a much higher density than expected as well.