Hello all. I have problem heating silicon with Tersoff potential. Even if i set velocity as low as 1 after several timesteps temeperature is around 1000 and system falls apart. Same happens if i heat it using nvt, changing timestep or dumping parameter doesnt help. Here are the codes:
#Velocity create:
units metal
atom_style atomic
boundary p p p
lattice diamond 5.43
region box block -10.0 10 -10.0 10 -10 10.0 units lattice
processors * * 1
create_box 2 box
region layer block -2 2 -2 2 -2 2
create_atoms 1 region layer
mass 1 28.0
mass 2 28.0
pair_style tersoff
pair_coeff * * Si.tersoff Si Si
minimize 1.0e-5 1.0e-8 10000 10000
velocity all create 10 4928459
fix 1 all nve
thermo 100
timestep 0.0005
dump 2 all xyz 100 out2.xyz
run 10000
#nvt heating:
units metal
atom_style atomic
boundary p p p
lattice diamond 5.43
region box block -10.0 10 -10.0 10 -10 10.0 units lattice
processors * * 1
create_box 2 box
region layer block -2 2 -2 2 -2 2
create_atoms 1 region layer
mass 1 28.0
mass 2 28.0
pair_style tersoff
pair_coeff * * Si.tersoff Si Si
minimize 1.0e-5 1.0e-8 10000 10000
reset_timestep 0
thermo 100
fix 1 all nvt temp 3000 3000.0 0.02
timestep 0.0002
dump 1 all xyz 100 out2.xyz
run 100000
Your structure has 6 not equilibrated surfaces, so it is unstable. You can, for example, add a short run with fix viscous before setting velocities.
The system falls apart in nvt heating, because it is set to 3000K.
Michal
i ran velocity down to 1e-14 with fix viscous and system still blows when i heat it.
i ran velocity down to 1e-14 with fix viscous and system still blows when i heat it.
i ran velocity down to 1e-14 with fix viscous and system still blows when i
heat it.
why should it not? what is the boiling point for your Si parameter set?
axel.
i mean, even if i heat it to as low as 10K it blows up. Here is the code:
#Velocity create:
units metal
atom_style atomic
boundary p p p
lattice diamond 5.43
region box block -10.0 10 -10.0 10 -10 10.0 units lattice
processors * * 1
create_box 2 box
region layer block -2 2 -2 2 -2 2
create_atoms 1 region layer
mass 1 28.0
mass 2 28.0
pair_style tersoff
pair_coeff * * Si.tersoff Si Si
minimize 1.0e-5 1.0e-8 10000 10000
fix 1 all nve
thermo 100
fix vsc all viscous 1
run 5000
unfix vsc
velocity all create 10 4928459
timestep 0.0005
dump 2 all xyz 100 out2.xyz
run 10000
i mean, even if i heat it to as low as 10K it blows up. Here is the code:
i see no "blowing up" at all. atoms stay in their location.
either you are not telling the whole story, or you have a modified
potential or you have a miscompiled or buggy executable.
I havent modified the potential and I’m actually running it on windows.
I think this is the whole story. My code ends up like this in case i didnt copy it correctly the last time
unfix vsc
velocity all create 10 4928459
timestep 0.0005
dump 2 all xyz 100 out2.xyz
run 10000
I havent modified the potential and I'm actually running it on windows.
I think this is the whole story. My code ends up like this in case i didnt
copy it correctly the last time
you didn't.
unfix vsc
velocity all create 10 4928459
timestep 0.0005
dump 2 all xyz 100 out2.xyz
run 10000
so if it look at the trajectory, the explanation is obvious and the one
that has already been given by michal way, way back in this discussion.
you have a very small object with unstable surfaces (lots of dangling
bonds). now the system will reconstruct and that will generate energy. that
energy is pushing specifically the corners of your cube so far that it is
enough to separate them from the rest.
if you want your system to be better behaved, you have to either generate
an initial structure with lower energy or use a dissipative thermostat with
a rather slow temperature ramp to slowly relax your structure, e.g.
something along the lines of this (you'll have to tweak the total number of
timesteps and the thermostat parameters to get the desired behavior):
thermo 100
velocity all create 10 4928459
timestep 0.0005
fix f0 all nve
fix f1 all temp/csvr 10.0 1000.0 0.01 12353
dump 2 all xyz 100 out2.xyz
run 100000
Thank you for your help.
By the way, I didnt say i havent got Michal’s argument.