Hello LAMMPS users,
I have a question regarding ground states, I am working with MEAM potentials, I am using a version of the example script to get the elastic constats, here it makes a minimimization as it follows:
I am assuming that with this kind of minimization I am obtaining the ground state of the structure, my main goal is to see a phase transformation for which I am using a script following some steps found in a paper where they see this transformation, the initial temperature is 500K and is decreased by 0.01K/step reaching 100k and then heated again to 500k, the script highlights goes as follows:
velocity all create 500 16723
Minimization (aforementioned)
fix 1 all nvt temp 500.0 100.0 0.1
run 40000
fix 2 all nvt temp 500.0 100.0 0.1
run 40000
Until the minimization everything is going ok but as soon as I start cooling down the structure is really distorted, but more important the energy obtained by the minimization is higher than the one showed at high temperature with the NVT ensemble that, in theory, it should match the high temperature energy.
I have tried starting from the low temperature structure, also I have added some timesteps at the limit temperatures to equilibrate the energy, I have used different box sizes and even tried using NPT ensemble but the results are similar.
Can someone “enlighten” me why does this happens? I have tried some potentials and with all of them I obtain similar behaviors regarding to the energy.
Thank you for taking the time to read and for your help,
Hello LAMMPS users,
I have a question regarding ground states, I am working with MEAM
potentials, I am using a version of the example script to get the elastic
constats, here it makes a minimimization as it follows:
I am assuming that with this kind of minimization I am obtaining the ground
state of the structure,
for that assumption to be true, you either need to have a very simple
structure or a *lot* of luck.
my main goal is to see a phase transformation for
which I am using a script following some steps found in a paper where they
see this transformation, the initial temperature is 500K and is decreased by
0.01K/step reaching 100k and then heated again to 500k, the script
highlights goes as follows:
velocity all create 500 16723
Minimization (aforementioned)
fix 1 all nvt temp 500.0 100.0 0.1
run 40000
fix 2 all nvt temp 500.0 100.0 0.1
run 40000
Until the minimization everything is going ok but as soon as I start cooling
down the structure is really distorted, but more important the energy
obtained by the minimization is higher than the one showed at high
temperature with the NVT ensemble that, in theory, it should match the high
temperature energy.
I have tried starting from the low temperature structure, also I have added
some timesteps at the limit temperatures to equilibrate the energy, I have
used different box sizes and even tried using NPT ensemble but the results
are similar.
Can someone "enlighten" me why does this happens? I have tried some
potentials and with all of them I obtain similar behaviors regarding to the
energy.
there is not enough information here to draw any conclusions and give
any recommendation. the simulation commands you quote, don't quite
match the simulation protocol you describe. for your 500K structure,
you not only need to assign an initial temperature, but also
equilibrate to that temperature. the minimization does the inverse.
then you have two cooling down temperature ramps. this raises the
suspicion, that there are more parts of your input that are
inconsistent or incorrect.
before anybody can help you for real, you need to show a full input
deck that is actually consistent with what you claim it would do.
Well regarding to my assumption of the ground state of the structure, for the high temperature I am working with a binary BCC structure and a monoclinic structure at low temperature, I have measured the elastic constants and obtained the parameters after the minimization for both structures, I believe BCC is quite a simple structure thats why I start with a high temperature in my script.
Here is my script, first I create at 500k, then I equilibrate, then I cool down, equilibrate, reheat and equilibrate again, with this I expect to see a fully phase transformation
#------------NVT part by steps ----------------------------------------
######## 1st step stabilizing @500K
fix 2 all npt temp 500 500 0.1 iso 0 0 1 drag 1
thermo 2000
thermo_style custom step pe temp lx ly lz press pxx pyy pzz c_eatoms cellalpha cellbeta cellgamma
dump 2 all atom 1000 T1.lammpstrj
dump 7 all custom 20000 *.T1.equil id type x y z c_csym
run 40000
unfix 2
undump 2
undump 7
######## 2nd step cooling down from 500K to 100K
fix 3 all npt temp 500 100 0.1 iso 0 0 1 drag 1
thermo 2000
thermo_style custom step pe temp lx ly lz press pxx pyy pzz c_eatoms cellalpha cellbeta cellgamma
dump 3 all atom 1000 Cooling.lammpstrj
dump 8 all custom 100000 *.cool.equil id type x y z c_csym
run 40000
unfix 3
undump 3
undump 8
######## 3rd step stabilizing @100K
fix 4 all npt temp 100 100 0.1 iso 0 0 1 drag 1
thermo 2000
thermo_style custom step pe temp lx ly lz press pxx pyy pzz c_eatoms cellalpha cellbeta cellgamma
dump 4 all atom 1000 T3.lammpstrj
dump 9 all custom 20000 *.T3.equil id type x y z c_csym
run 40000
unfix 4
undump 4
undump 9
######## 4rd step heating from 100K to 500K
fix 5 all npt temp 100 500 0.1 iso 0 0 1 drag 1
thermo 2000
thermo_style custom step pe temp lx ly lz press pxx pyy pzz c_eatoms cellalpha cellbeta cellgamma
dump 5 all atom 1000 Heating.lammpstrj
dump 10 all custom 100000 *.heat.equil id type x y z c_csym
run 40000
unfix 5
undump 5
undump 10
######## 5th step stabilizing @ 500K
fix 6 all npt temp 500 500 0.1 iso 0 0 1 drag 1
thermo 2000
thermo_style custom step pe temp lx ly lz press pxx pyy pzz c_eatoms cellalpha cellbeta cellgamma
dump 6 all atom 1000 T5.lammpstrj
dump 11 all custom 20000 *.T5.equil id type x y z c_csym
run 40000
unfix 6
undump 6
1. If I understand correctly, you are trying to get the crystal to
switch from BCC to monoclinic.
2. Solid-solid phase transformations are difficult to reproduce,
because it is very easy to get stuck in all kinds of kineticly trapped
configurations that are not at thermodynamic equilibrium. So there is
no guarantee that a simulation will achieve this goal, even if it is
free of obvious errors.
3. However, the above script is guaranteed not to achieve this goal,
because you used the keyword "iso". To see a shift to monoclinic, you
will need to use "tri"
4. You might be able to force the phase transition by applying a
suitable strain while running nvt dynamics. See fix deform.
Thank you very much for your answer,
I haven’t realized about the “tri” option since I was starting from B2, thank you for your recommendation.
I still have 2 questions,
As Axel mentioned to get the ground structure threw a minimization its a matter of having a simple structure or a lot of luck, what would be a good way of doing it? a combination of NVE energy equilibration and NPT pressure equilibration would work?
To get the elastic constants of the “heated” cubic structure is it a good way to first equilibrate the structure at the temperature and then read the elastic constants in the same way described in the example included by LAMMPS (perhaps without the minimization parts)?
Thank you very much for your answer,
I haven't realized about the "tri" option since I was starting from B2,
thank you for your recommendation.
I still have 2 questions,
1. As Axel mentioned to get the ground structure threw a minimization its
a matter of having a simple structure or a *lot* of luck, what would be a
good way of doing it? a combination of NVE energy equilibration and NPT
pressure equilibration would work?
There is no guarantee that any simulation will generate the lowest energy
crystal all by itself. It is a lot easier to invert the problem. Build
known crystal structures and relax them with your potential at P = 0, T = 0
using fix box/relax, and see which one has the lowest potential energy.
2. To get the elastic constants of the "heated" cubic structure is it a
good way to first equilibrate the structure at the temperature and then
read the elastic constants in the same way described in the example
included by LAMMPS (perhaps without the minimization parts)?
First start by getting elastic constants at T = 0, P = 0 using the LAMMPS
example in examples/ELASTIC. Then try using the T > 0 LAMMPS example in
examples/ELASTIC_T. The reason is that there are many more things that can
go wrong for T > 0.