In metal units your damping constant is 30, which
means 30 psec, which is 30,000 timesteps. So
your Temp will equilibrate on that timescale.
Probably 100x to 1000x longer than you want.
'
Steve
Hello,
I am interested in studying melting of bulk, 2D and 1D Silicon structures. I am using SW potential based on reading some papers.
For studying bulk Silicon, I have tried two strategies :
-
ramping the temperature of a bulk Si block (1000 atoms with ppp boundary conditions) at a rate of 1K/1ps in steps of 100K and equilibrating at each temperature for 2ns using a series of fix NPT runs. However, even when the simulation ramps to 2000K there are no signs of melting (MSD, volume or coordination do not change)
-
take two blocks of Si (1000 atoms with ppp boundary condition) where one is solid Si and another is liquid (started with solid Si and ramped the temperature upto 3000K) and then study the system near the melting point. These were independently run with NPT. This two phase system is run under NVT and I see something sensible from this simulation.
But, given approach 1 does not give the correct melting point, I was wondering how to setup a simulation of 2D block of Silicon (boundary conditions are no longer ppp) given that NPT in LAMMPS does not work with non-periodic dimensions. If I setup a system with “ppf” boundary conditions (or pff for 1D case) and run NVT, will I get reasonable results for melting ?
can anyone help me understand how to setup simulations to study melting of 1D/2D structures ?
thanks,
Anurag
Hello,
I am interested in studying melting of bulk, 2D and 1D Silicon structures.
I am using SW potential based on reading some papers.For studying bulk Silicon, I have tried two strategies :
1. ramping the temperature of a bulk Si block (1000 atoms with ppp
boundary conditions) at a rate of 1K/1ps in steps of 100K and equilibrating
at each temperature for 2ns using a series of fix NPT runs. However, even
when the simulation ramps to 2000K there are no signs of melting (MSD,
volume or coordination do not change)
no surprise there. melting (and freezing) is a process with a
hysteresis and usually needs some kind of nucleation (e.g. from
internal fluctuations). you very likely have to deal with massive
finite size effects here.
2. take two blocks of Si (1000 atoms with ppp boundary condition) where one
is solid Si and another is liquid (started with solid Si and ramped the
temperature upto 3000K) and then study the system near the melting point.
These were independently run with NPT. This two phase system is run under
NVT and I see something sensible from this simulation.But, given approach 1 does not give the correct melting point, I was
wondering how to setup a simulation of 2D block of Silicon (boundary
conditions are no longer ppp) given that NPT in LAMMPS does not work with
non-periodic dimensions. If I setup a system with "ppf" boundary conditions
what you describe as a 2D system would actually still be a 3D system,
only that it is not periodic in the 3rd dimension. so just set up an input
with ppp boundary and then leave a sufficiently large (at least 2x cutoff)
space between the periodic images and then run npt, but keep the 3rd
dimension at fixed length (fix npt supports that). this is generally referred
to as a slab geometry.
similarly you can set up a 1D periodic rod geometry.
(or pff for 1D case) and run NVT, will I get reasonable results for melting
?can anyone help me understand how to setup simulations to study melting of
1D/2D structures ?
there has to be a _ton_ of literature on simulation and theory/models
of melting and crystallization. most of which will probably tell you that
it is pretty difficult to get it done right.
axel.
Thanks Dr. Kohlmeyer.
I have a follow-up question:
in relation to approach 1, you mentioned about finite size effect due to small box size I am using. does this mean, simulating a much bigger block of Si using approach 1 would give melting temperature close to experimental Tm ?
thanks,
Anurag
Thanks Dr. Kohlmeyer.
I have a follow-up question:in relation to approach 1, you mentioned about finite size effect due to
small box size I am using. does this mean, simulating a much bigger block
of Si using approach 1 would give melting temperature close to experimental
Tm ?
it should, however, the necessary size to see this happen spontaneously,
might be beyond any currently available supercomputing capacity.
there is also the question of time scales, i.e. how long does it take
until melting starts spontaneously?
i have not dealt with that personally, but i know that since quite a
long time ago, different researchers have proposed different
approaches for a "theory of melting" (and probably fought about it).
that is why i suggested to have an extended look at available
literature (search for the phrase "theory of melting")
more commonly know is the inverse problem, i.e. that
it is experimentally possible to have supercooled liquids
and maintain that state for a long time. you can try with a
bottle of water in a freezer, for example. if you don't cool
it too far below the freezing point, it can stay liquid for
quite a while or until you shake it (careful, water expands
on freezing and the bottle can thus explode).
axel.
Hello Prof. Kohlmeyer,
based on your suggestion I tried to setup the 2D/1D structures in LAMMPS but I am seeing some strange errors. these errors are common to both 2D/1D cases and I am not able to understand the origin of these errors. I have copied below the portion of input file I changed (from bulk case) to simulate slab/wire structures. can you please help me understand the source of error ?
regards,
Anurag
PS: I observed that if I use a 5x5x5 size bulk (ppp boundary), LAMMPS generates 1000 atoms. however for 1D case (as shown below) LAMMPS generates 1105 atoms. why ?
when I run LAMMPS in parallel as mpirun -np 4 lammps -in filename
I get the following error
Hello Prof. Kohlmeyer,
based on your suggestion I tried to setup the 2D/1D structures in LAMMPS but
I am seeing some strange errors. these errors are common to both 2D/1D
cases and I am not able to understand the origin of these errors. I have
neither am i. my guess would be that there is something
wrong in the parts of the input that you have omitted.
i have no psychic capabilities and my crystal ball is
broken beyond repair.
if you don't make it easy to reproduce any problems
that you are seeing or that you have questions about,
then your chances to get help will drop massively.
copied below the portion of input file I changed (from bulk case) to
simulate slab/wire structures. can you please help me understand the source
of error ?
not as such.
regards,
AnuragPS: I observed that if I use a 5x5x5 size bulk (ppp boundary), LAMMPS
generates 1000 atoms. however for 1D case (as shown below) LAMMPS generates
1105 atoms. why ?
probably some atoms are exactly on the simulation cell
boundaries and in the 3d case they are removed from
one face because of that, but not in the 1d case.
axel.
Thanks again for your suggestions.
mpirun invocation of lammps works now. I placed the wire in the center of box and specified the processors variable.
I have a question regarding barostat.
I have requested pressure control only along x-direction (axial dimension of wire). When I dump the simulation box dimension using dump custom, I check that the simulation box is contracting in the x-dimension but it is fixed in the y and z.
however, the mean square displacement of atoms increases monotonically and after visualizing the trajectory with VMD, it appears that the solid is reaching a disordered state. is this barostat setting ok ?
thanks and regards,
Anurag
I have a question regarding barostat.
I have requested pressure control only along x-direction (axial dimension of
wire). When I dump the simulation box dimension using dump custom, I check
that the simulation box is contracting in the x-dimension but it is fixed in
the y and z.however, the mean square displacement of atoms increases monotonically and
after visualizing the trajectory with VMD, it appears that the solid
is reaching a disordered state. is this barostat setting ok ?
how the hell should i know??
how many times do you want me to repeat that i cannot
read minds. neither yours nor that of your computer.
axel.
sorry for describing the details.
please find below the input file and the various fix npt options I have tried so far. LAMMPS output is also excerpted.
----------- input --------------
Silicon wire (1D) melting
units metal
boundary p p p
atom_style atomic
processors * 2 2
lattice diamond 5.431
region box block 0 5 0 15 0 15
#leave a vacuum space in y and z
create_box 1 box
region wire block 0 5 5 10 5 10
create_atoms 1 region wire
mass 1 28.0855
timestep 0.002
pair_style sw
pair_coeff * * Si.sw Si
neighbor 0.5 bin
neigh_modify every 5 delay 0 check yes
velocity all create 300 825577 dist gaussian
dump 4 all xyz 1000 wire_si.xyz
compute 3 all pe/atom
compute 4 all ke/atom
compute 5 all coord/atom 3.0
compute 6 all cna/atom 2.5
compute 7 all msd com yes
thermo 1000
thermo_style custom step temp vol etotal pe ke pxx pyy pzz c_7[4]
dump 1 all custom 5000 wire_si.atom id xs ys zs c_3 c_4 c_5 c_6
fix 1 all npt temp 300 300 0.3 x 0.0 0.0 0.3 y NULL NULL 0.3 z NULL NULL 0.3
run 500000
I meant sorry for NOT describing the details earlier.
dear anurag,
please let me know if further information is required. can you please
suggest what is the problem ?
the biggest problem at this point is that you just want
me to tell you what to do so that you don't have to do
any thinking and testing for yourself. this is not how
mailing lists are supposed to work. i don't have the
time and the interest to play 20 questions with people.
i am only interested to help people that want to
understand and learn from their problems so that
they don't need so much help in the future. so to
keep my interest, you have to show me that you are
not only blindly following what i tell you, but that you
are drawing conclusions from your mistakes and
that you learn how to tell what is wrong from right
and not just do something and then have somebody
else give you the thumbs up or down.
half of the answer you are looking for, is in the
description you gave in your mail. the other half is
in the answers i gave you before. re-read them,
draw your conclusions, think about what you want
to achieve and what you are doing. i _strongly_
recommend to carefully monitor your system with
a visualization tool on top of just looking at some
simplified parameter like MSD.
axel.
Hello,
I have a question regarding LAMMPS command : read_data.
Documentation mentions that pair_coeff must be specified after the simulation box is defined using read_data, … while pair_style must be specified before coefficients are set by pair_coeff, read_data, …
I have simple input file which uses read_data to read initial atomic coords. LAMMPS, however, reports error in the way pair_coeff, pair_style commands are specified. any suggestions where am I going wrong ?
thanks,
Anurag
PS: I checked the example file in.meam where pair_style and pair_coeff are specified after read_data but when I try that I again get an error.
units metal
boundary p p p
atom_style atomic
pair_style sw
#pair_coeff * * Si.sw Si
read_data initial.data
timestep 0.001
#pair_style sw
pair_coeff * * Si.sw Si
neighbor 0.5 bin
neigh_modify every 5 delay 0 check yes
velocity all create 1700 825577 dist gaussian
dump 4 all xyz 10000 solid_si.xyz
thermo 1000
thermo_style custom step temp vol etotal pe ke press enthalpy
dump 1 all custom 5000 solid_si.atom id xs ys zs c_3 c_4 c_5 c_6
fix 1 all npt temp 1700 1700 0.3 iso 0.0 0.0 0.2
run 100000
errors reported by LAMMPS: ERROR: Incorrect args for pair coefficients
Hello,
I would like to know if there is a way to extract elastic strain energy from LAMMPS output.
I am considering a two layer system (layer B grown on layer A) where initially atoms are placed at lattice sites corresponding to lattice constant of layer A. The system is then allowed to relax whereby layer B distorts tetragonally. material A and B have cubic lattices.
how to calculate the strain energy stored in this relaxed epilayer (layer B) ?
Thank you.
best regards,
Anurag
How is that different than the total energy of the system or the per-atom
energy summed over some subset of atoms?
Steve
I will explain the source of my confusion.
elastic strain energy for biaxially grown layers (assuming Si as substrate and SiGe as the grown coherent layer) is given using eq. 14 from (D. J. Dunstan, J. of Mater. Sci, vol 8, pp 337, 1997)
as E = 0.5 * (volume summation over stress*strain) where stress and strain are expressed as tensors.
or per unit area using eq. 15 from the same paper as
E = thickness of grown layer (h) * biaxial Young’s modulus (M) * (misfit strain or e_0)^2
If I use SW potential to model Si-Si, Si-Ge and Ge-Ge interaction with LAMMPS, I get lattice parameter for Si_0.3Ge_0.7 as 5.58423 A. Si lattice constant is 5.431A so with this information I can calculate e_0.
again using LAMMPS prescribed method to calculate elastic constants, I get C11= 139.69GPa and C12 = 52.98GPa. which gives me M = 152.482 GPa (using eq. 13 from the paper).
now plugging numbers for a 4.5nm Si_0.3Ge_0.7 layer, I get ~0.549 Joules per sq. meter as the elastic strain energy.
second approach is to simulate a structure with Si_0.3Ge_0.7 layer on top of Si (100) such that it is biaxially strained to have Si lattice constant in the interface plane. when the system is allowed to relax, I see that the SiGe layer distorts tetragonally in the (100) direction. strain tensor extracted using MD data is also reasonable. I take the total energy of the relaxed structure and subtract from it (# of atoms in the SiGe layer * energy per atom of the bulk SiGe) to get “strain energy”. I am not sure if this the right approach.
then to compare, I divide this strain energy by the interface plane area used in the MD simulation. this gives ~2.26 Joules per sq. meter. I do not know how to account for this >4x difference between the two approaches.
any comments/suggestions are welcome. thank you.
regards,
Anurag
E = 0.5 * (volume summation over stress*strain) where stress and strain are expressed as tensors
I don't know what this eq means for an atomistic simulatoin. Sum over what?
What is the volume? All you have are atom's and their per-atom stress
or possibly strain if that means their displacement.
I don't think I can help you further. The LAMMPS doc pages for compute pressure
and compute stress/atom and compute displace/atom are quite clear about
what they calculate. I think you will need to figure out if that is a match
to these papers.
Steve
Hello,
I am looking to study a Si nanowire under uniaxial compression test at different temperatures.
In looking at previous work on LAMMPS mailing list I came across few posts where someone proposed using “fix uniaxial” command however, it seems that this fix is obsolete now. can someone indicate how such a simulation can be setup with the current version of LAMMPS ?
thank you.
sincerely,
Anurag
I don't recall there was ever a fix uniaxial command.
Take a look at fix deform - it can do tenslie and compressive
changes to the LAMMPS simulation box.
Steve
As per LAMMPS documentation on fix deform …"Any parameter varied by this command must refer to a periodic dimension "
so if we have to compare nanowires with different axial length then it would not be possible to do so using fix deform.
thanks for your help.
regards,
Anurag