Amorphous silica from cristobalite melting & quenching - morse/coul/long potential

Deal all,

I’m trying to create amorphous silica by melting and quenching cristobalite, with a morse + coul/long hybrid potential.

I use NPT system.

My concern is that I cannot equilibrate my system : the pressure is very high (several thousands of bar instead of 1).

I tried to play with the Tdamp, Pdamp, drag options ==> no effect

I tried to play on the high temperature step duration ==> no effect

I tried to decrease the temperature step from 7000K to 3000K ==> no effect

I tried reneighboring more often ==> no effect

I confess I’m running lack of ideas…

My script is below : does anyone has any comment ? Is there an ugly mistake in that script ?

Best regards,

Alexandre

Deal all,

I'm trying to create amorphous silica by melting and quenching cristobalite,
with a morse + coul/long hybrid potential.

I use NPT system.

My concern is that I cannot equilibrate my system : the pressure is very
high (several thousands of bar instead of 1).

I tried to play with the Tdamp, Pdamp, drag options ==> no effect

I tried to play on the high temperature step duration ==> no effect

I tried to decrease the temperature step from 7000K to 3000K ==> no effect

I tried reneighboring more often ==> no effect

I confess I'm running lack of ideas...

My script is below : does anyone has any comment ? Is there an ugly mistake
in that script ?

yes. your time constant for the cell relaxation is huge. next to
nothing will happen over the duration of your trajectory.

axel.

Pdamp is above the recommended value because with Pdamp=1000*timestep, my system explodes.
Surely a silly reason I must admit...

Going to a recommended value for Pdamp and running O gives me the well known Out of range atoms.
And I cannot figure out why for the moment (forces are small, I think no overlapping atoms), so I believed it was an artefact of the first timestep, maybe I'm wrong about it.
I'll check more in details for that huge pressure. Or maybe minimization with box/relax, even if I dont understand why a beta-cristobalite structure seems so far from equilibrium...

Anyway thank you for your quick answer,
Bye

Alexandre

-----Message d'origine-----

Pdamp is above the recommended value because with Pdamp=1000*timestep, my system explodes.
Surely a silly reason I must admit...

Going to a recommended value for Pdamp and running O gives me the well known Out of range atoms.
And I cannot figure out why for the moment (forces are small, I think no overlapping atoms), so I believed it was an artefact of the first timestep, maybe I'm wrong about it.

have you checked for overlapping atoms due to periodic boundaries?

also, have you used the same force field parameters before?
perhaps you have a typo or mixed up a unit conversion somewhere.

also, try using a shorter time step. perhaps some of the repulsive
parts of the interactions are too steep for your choice.

I'll check more in details for that huge pressure. Or maybe minimization with box/relax, even if I dont understand why a beta-cristobalite structure seems so far from equilibrium...

before doing anything about relaxing the box, you should first see to
get a decent stable simulation with a fixed volume.

axel.

I began with NVT and the pressure was very high (~1e5 bar), so I thought that accommodating the box with NPT could help... wrong guess

Concerning your remarks :
- timestep : I decreased it by factor 10, no change
- concerning the ff : I never used it before, I double-checked the value, I believed it is ok (to be sure I jumped to units real, so that no unit conversion is needed)
- concerning the cristobalite : I used a file taken from that list, written by Paul, and I could run a step at 1K and 1 bar (which is not possible with my beta cristobalite).

so : I have to check for that cristobalite

but : my system still explodes with a step at 7000K with Paul's initial configuration
do you think that a first step at 7000K is too hard for the system ? (since it seems ok at 1K)
do I have to ramp the temperature from 1K to 7000K ?

that very huge pressure bothers me, but I feel that I'm not in the right direction in playing with ramps

Alex

-----Message d'origine-----

Things to check:

Why do you have 16 atoms in the unit cell? Can
you run with a normal diamond lattice with 8 atoms/unit cell?
Can you run NVE at a normal low temperature (300K)
and run a stable system with reasonable pressure? If not,
something is wrong with your model or initial conditions.

Steve

Steve,

I fear I don't understand your first question very well : I have 16 atoms corresponding to oxygen (lattice custom) and 8 for silicon (lattice diamond). Did I misunderstand the way lattices are created in Lammps ?

Concerning the second point, I run NVE with temp/rescale @ 300K for 1ns and the system seems stable (ie does not evolve too much in terms of force) with pressure oscillating between -6000 and 6000 bars...
I have the same oscillations amplitudes (but this times not symmetric around 0) with Tersoff or BKS potentials.
To me that means that the force field parameters are not the cause of the pb.

To go further, I run a ramp from 300 to 7000K with NVE + NH thermostat:
- with Tersoff, I see a transition (E jumps up and pressure jumps down and then is nearly constant) @ 3800K ==> melting, I'm happy (melting temp is a bit high indeed, but...)
- with Morse/coul/long, I see a transition in E @ 3000K (I'm happy) BUT no (or very tiny) pressure decrease + further pressure increase after the transition and *that*, I dont understand...

I did not run Verlet + thermostat + barostat but I suspect a correct decrease in volume for Tersoff and a volume blow with morse/coul/long...
I'll run it to be sure but I'm *quite* pessimistic about it.

So this potential must be the pb, whereas I saw several authors using it with the set of parameters I use (even in Lammps : http://lammps.sandia.gov/abstracts/abstract.2416.html)

Alexandre

-----Message d'origine-----

Alex,
If your potential is a good one, you should recover alpha-quartz as
the lowest energy polymorph of SiO2
at low temp. Have you verified this?
Carlos

I fear I don't understand your first question very well : I have 16 atoms corresponding to oxygen (lattice custom) > and 8 for silicon (lattice diamond). Did I misunderstand the way lattices are created in Lammps ?

I don't know. I think of a diamond lattice as a cubic lattice with 8
basis atoms. But there are many
unit cells with basis atoms that can be used. So long as you are sure
you have the right
geometry and density, you should be OK.

I also don't know why you are having problems with the Morse
potential. It may help to monitor
the thermo output at a high frequency of output to see if you are getting large
temperature or volume excursions. That often indicates bad dynamics,
e.g. too big a
timestep.

Steve

A quick message for the list to say that I managed to quench my Si - O mixture (either beta cristobalite or random) to amorphous SiO2 with pretty good RDF and angles.
It seems that it was mainly a question of melting temperature : decreasing it to 4000K works fine with NPT dynamics.
Alternatively, melting at 7000K and quenching to moderate temperature with NVT dynamics, then quenching to ambient with NPT dynamics gives the same final structure.

Thank you again for your help

Alexandre

PS : if anyone wants the inputs, just email me

-----Message d'origine-----