I am having some trouble with Tersoff potential when coupling berendsen barostat and NVT ensemble. When I minimize at 0K and 0bars, I get a density of 2.42 g/cc. When I run NVT with barostat at 0bars, density decreases as it should. Can anyone physically explain why when minimizing to 1atm or 1.01325 bars, my density is about 2.42 g/cc still, but upon running NVT with barostat, the density increases up to about 2.6 g/cc?
Ben
#Equilibration for quartz using the Tersoff Potential
units metal
atom_style full
#Data format is read in below
read_data …/…/quartz.data
####################Potential#############################
pair_style tersoff
pair_style tersoff
pair_coeff * * …/…/…/…/potentials/SiO.tersoff O Si
##############Minimization###############################
fix 1 all box/relax x 1.01325 y 1.01325 z 1.01325 vmax 0.001
minimize 1.0e-50 1.0e-50 100000 1000000
unfix 1
Axel has explained this to you before - it is likely that the data
structure is metastable. A potential energy surface may have multiple
minima, and minimization and geometry optimization may not take the
metastable structure away from its original local minimum. Once
temperature is supplied the atoms may have a better chance evolving to
a lower energy minimum (local or global), and this may be why the
density increased a little bit after equilibrating at 300K and 1 atm.
Ray
I understand, but why did it do the same at 0atm?
Ben
Possibly due to similar reason. There are some problems with your script.
Ray
I understand, but why did it do the same at 0atm?
Ben
Axel has explained this to you before - it is likely that the data
structure is metastable. A potential energy surface may have multiple
minima, and minimization and geometry optimization may not take the
metastable structure away from its original local minimum. Once
temperature is supplied the atoms may have a better chance evolving to
a lower energy minimum (local or global), and this may be why the
density increased a little bit after equilibrating at 300K and 1 atm.
Ray
> I am having some trouble with Tersoff potential when coupling berendsen
> barostat and NVT ensemble. When I minimize at 0K and 0bars, I get a
> density
> of 2.42 g/cc. When I run NVT with barostat at 0bars, density decreases
> as it
> should. Can anyone physically explain why when minimizing to 1atm or
> 1.01325
> bars, my density is about 2.42 g/cc still, but upon running NVT with
> barostat, the density increases up to about 2.6 g/cc?
>
> Ben
>
> #Equilibration for quartz using the Tersoff Potential
> units metal
> atom_style full
You don't need full for Tersoff.
>
> #Data format is read in below
> read_data ../../quartz.data
>
> ####################Potential#############################
>
> pair_style tersoff
> pair_style tersoff
Duplicated entry.
> pair_coeff * * ../../../../potentials/SiO.tersoff O Si
>
> ##############Minimization###############################
>
> fix 1 all box/relax x 1.01325 y 1.01325 z 1.01325 vmax 0.001
> minimize 1.0e-50 1.0e-50 100000 1000000
Minimizing to e-50 is meaningless.
> unfix 1
>
> #########################################################
>
>
> fix 2 all nvt temp 300 300 0.1
> fix 3 all press/berendsen x 1.01325 1.01325 100 y 1.01325 0 100 z
> 1.01325
> 1.01325 100
pdamp of 100 ps is definitely too large.
Also the default timestep of 1 fs may be too large for Tersoff.
I understand that I do not need full for Tersoff but does it hurt anything to have the charges in the data file? They wont play a role since I ignore Columbic forces.
Ben
Not just charge, you also don't need molecule_ID. If something is not
needed, why keep it there?
I am having some trouble with Tersoff potential when coupling berendsen
barostat and NVT ensemble. When I minimize at 0K and 0bars, I get a density
of 2.42 g/cc. When I run NVT with barostat at 0bars, density decreases as it
should. Can anyone physically explain why when minimizing to 1atm or 1.01325
bars, my density is about 2.42 g/cc still, but upon running NVT with
barostat, the density increases up to about 2.6 g/cc?
we have been over this already and established that a) your initial
structure was metastable (minimizing usually does not break symmetries
for analytical potentials) and b) that after running for a bit your
structure changes significantly (i.e. you have a phase change) during
MD and then you get - obviously - different properties.
i fail to see, why this is so difficult to understand and why you need
to keep wasting your time with pointless simulations that do nothing
but re-establish those facts.
one more thing i am noticing - and that goes along the of basic
understanding of what simulations and MD do - is that you are
mistaking a model with the reality. models can have all kinds of
limitations and specific parameterizations are good at representing
parts of a phase diagram of a material and others are less so and all
of them have limitations. so what you are looking at are less the
properties of the material, but rather the properties of the
potential. if those properties match the material, is a measure for
the quality of the potential, but also a measure for how meaningful
your simulation setups are. in this specific case, you *first* need to
make a choice what kind of phase (or structure) you want to model
exactly and then you validate your potentials/model in two ways
(against the reference structures provided in the paper describing the
potential and against your desired structure) to determine correctness
of the implementation and transferability. your rather random
simulation attempts and almost obsessive desire to just change
parameters until you get the numbers that you desire to see are not
helpful in that regard.
if you want to get a feel for what good or bad a model can be, you
should rather do a DFT calculation (and be surprised by how far those
are away from the perfect experimental data depending on the
functional used). with tersoff, you should also consider the
complexity of the model and how difficult it is to balance all those
parameters to the different components of the model, especially for
structure like SiO2. in fact, i am not fully convinced that plain
tersoff is a good model for this in the first place, since it is based
on some assumptions that are not really true for a material like SIO2.
there are other manybody potentials that do capture this better.
please also note, that you can have error cancellation, particularly
through highly symmetric input structures. another thing that you seem
to be completely ignoring. finally, you may need to consider the
effect of cutoffs and how those affect structure, specifically if the
number of atoms that are nearest neighbors changes due to a change in
the structure.
all of this have nothing to do with LAMMPS, but are really general
considerations related to doing simulations of a material overall.
trying to find correlation to certain input flags in LAMMPS will
likely not provide much insight. if anything, it will be confusing,
since the correlations are different.
axel.
Axel, I know I am lacking basic MD knowledge. I never really understood this phase change. I heat to 300K over 10ps and no phase change. Then I equilbrate at 300K for 10ps and then a phase change. How can there be a phase change in this situation? I am not heating anything at the point of phase change, just trying to look at properties at 300K. I am doing these simulations because I obviously lack this knowledge and need to learn it. If I went on to more complex models without understanding this one, I would not get very far.
Ben
Axel, I know I am lacking basic MD knowledge. I never really understood this
why don't you talk to your adviser? or to colleagues? those are the
ones that can see what you do and discuss with you. by trying things
all on your own you are massively wasting your time.
there are times when it is good exercise to do something on your own,
but you also need to be smarter about when to ask people for help. as
a rule of the thumb, you need to get basic/fundamental knowledge and
experience with local help, and go to external resources for specialty
problems or technical issues.
phase change. I heat to 300K over 10ps and no phase change. Then I
equilbrate at 300K for 10ps and then a phase change. How can there be a
phase change in this situation?
2 x 10ps is nothingness in terms of MD time. a blink of an eye.
as was pointed out, you likely start from a geometry that is - with
respect to your potential(!) - metastable. note that phase changes are
activated processes and - particularly for very small systems - you
can have a comparatively large hysteresis.
I am not heating anything at the point of
phase change, just trying to look at properties at 300K. I am doing these
simulations because I obviously lack this knowledge and need to learn it. If
I went on to more complex models without understanding this one, I would not
get very far.
yes, but you *have* to understand that this is not a LAMMPS problem.
you need to attend classes, courses, and get tutoring in person. no
mailing list can ever replace this.
axel.