NPT equilibration of Silicon

Dear Lammps users,

I have been running into this problem for long time now, whenever the equilibration starts the beam starts shaking. When I use the same code for Aluminum it works fine. Any suggestions are welcom. Thanks

------------------------ INITIALIZATION ----------------------------

echo screen
units metal
dimension 3
boundary p p p
atom_style atomic

----------------------- ATOM DEFINITION ----------------------------

variable latparam equal 5.431

variable nL_x equal -65.431
variable L_x equal 6
5.431
#x= nm

variable nL_y equal -95.431
variable L_y equal 9
5.431

y= nm

variable nL_z equal -85.431
variable L_z equal 8
5.431

z= nm

variable L_fixed equal (-5.55.431)
variable L_mobile equal (5.5
5.431)

lattice diamond {latparam} region beam block {nL_x} {L_x} {nL_y} {L_y} {nL_z} ${L_z} units box
region whole block -100 100 -100 100 -100 100 units box
create_box 1 whole

lattice diamond ${latparam} orient x 1 1 0 orient y 0 0 1 orient z 1 -1 0
create_atoms 1 region beam
variable n equal count(all)

------------------------ FORCE FIELDS ------------------------------

pair_style meam
pair_coeff * * /usr/local/lammps-9Dec14/
potentials/library.meam Si /usr/local/lammps-9Dec14/potentials/Ni.meam Si

neighbor 2.0 bin
neigh_modify delay 0 every 2 check yes

---------------------------Minimization-------------------------

#minimize etol ftol maxiter maxeval

Conjugate gradient Algorithm

min_style cg
minimize 1.0e-7 1.0e-7 100000 100000

variable tmp equal “lx”
variable L0 equal {tmp} print "Initial Length, L0: {L0}"

------------------------------Equilibration---------------------

reset_timestep 0
timestep 0.001

velocity all create 300 12345 mom yes rot no
fix 2 all nvt temp 300.0 300.0 1.0

Set thermo output

thermo 1000
thermo_style custom step lx ly lz press pxx pyy pzz pe ke etotal temp
dump 1 all atom 100 dump.eq.lammpstrj

write_restart Compress.equil

Run for at least 10 picosecond (assuming 1 fs timestep)

run 20000
undump 1
unfix 2

variable tmp equal “lx”
variable L0 equal ${tmp}
variable LZ equal “lz”

Strike two for you…

How about putting more effort into PROPERLY explaining what your problem is?

The beam starts shaking…what is that supposed to tell people with no clue which model/system/problem you are trying to simulate?

Carlos

I am sorry, I am trying to simulate a Single crystal Silicon beam bending. So I need to have vacuum around it from all directions. I made the simulation box a little larger and used periodic boundary conditions. I am using NPT during equilibration, the temp. doesn’t fluctuate too much neither the potential energy. But when I visualize the atoms they are vibrating so fast and the beam itself is shaking with large deformations. I alos tried the NVT but also getting similar behavior. the log file results for NVT are attached here. Actually i used to have F F F boundaries and simulate using langevin and NVE, but I want to check if I can get the same physics in case of Nose Hoover. Thanks

Mohamed

Step Lx Ly Lz Press Pxx Pyy Pzz PotEng KinEng TotEng Temp
0 120 140 120 1342.733 1263.0946 1414.0579 1351.0464 -120570.82 1051.6605 -119519.16 300
1000 120 140 120 -668.31786 98.343185 -833.35806 -1269.9387 -120085.22 591.82961 -119493.39 168.82718
2000 120 140 120 -252.10232 -1470.0148 120.36893 593.3389 -120145.04 689.97255 -119455.07 196.82374
3000 120 140 120 464.76278 1631.2318 451.67346 -688.61694 -120151.34 740.53289 -119410.8 211.24675
4000 120 140 120 520.06845 -50.530045 244.60338 1366.132 -120131.55 768.28043 -119363.27 219.1621
5000 120 140 120 -1005.9857 -978.25952 -757.53444 -1282.1632 -120122.26 808.626 -119313.63 230.6712
6000 120 140 120 954.37596 545.91975 1192.6675 1124.5406 -120121.31 858.49015 -119262.82 244.89561
7000 120 140 120 -24.809746 956.83289 -344.30647 -686.95566 -120099.71 887.92154 -119211.78 253.2913
8000 120 140 120 -141.95572 -1392.4993 373.41749 593.21464 -120054.1 892.98846 -119161.12 254.7367
9000 120 140 120 61.02631 1060.8225 -176.61273 -701.13087 -120032 920.66688 -119111.34 262.63234
10000 120 140 120 533.96338 278.07675 345.16988 978.6435 -120016.84 954.14205 -119062.69 272.18157
11000 120 140 120 -891.97121 -947.05369 -619.61541 -1109.2445 -119997.4 980.62516 -119016.78 279.73622
12000 120 140 120 735.50048 407.89046 756.1805 1042.4305 -119977.03 1003.2161 -118973.81 286.18058
13000 120 140 120 112.09825 999.80311 -129.6897 -533.81865 -119953.81 1018.5877 -118935.23 290.56552
14000 120 140 120 -470.37261 -1671.2981 -233.96757 494.1478 -119935.39 1033.8438 -118901.54 294.91754
15000 120 140 120 -234.57018 573.86026 -347.16897 -930.40182 -119937.06 1062.7586 -118874.3 303.16588
16000 120 140 120 640.05591 652.11547 527.84098 740.21127 -119930.1 1075.2258 -118854.88 306.72229
17000 120 140 120 -836.8541 -910.29783 -1005.9458 -594.31864 -119911.48 1065.4473 -118846.03 303.93286
18000 120 140 120 112.83116 -37.644056 377.79204 -1.6544941 -119925.58 1077.247 -118848.33 307.29887
19000 120 140 120 652.16576 1106.2861 372.68384 477.52736 -119943.39 1081.4633 -118861.92 308.50163
20000 120 140 120 -532.39316 -676.13843 -680.12597 -240.91507 -119938.39 1056.4193 -118881.97 301.3575

I am sorry, I am trying to simulate a Single crystal Silicon beam bending.
So I need to have vacuum around it from all directions. I made the
simulation box a little larger and used periodic boundary conditions. I am
using NPT during equilibration, the temp. doesn't fluctuate too much neither
the potential energy. But when I visualize the atoms they are vibrating so
fast and the beam itself is shaking with large deformations. I alos tried

i ran your input and it is more "wobbling" instead of "shaking".

the NVT but also getting similar behavior. the log file results for NVT are
attached here. Actually i used to have F F F boundaries and simulate using
langevin and NVE, but I want to check if I can get the same physics in case
of Nose Hoover. Thanks

the problem is not the thermostat. the problem is in your
expectations. if you cut through Si-Si bonds and have a free surface
with unsaturated dangling bonds, where previously was bulk, how can
you not expect that there will be a significant reconstruction
happening which silicon is well known for. also, how certain are you
that MEAM is a suitable potential for bulk *and* surfaces of Si? i
would have expected you'd be using Tersoff or Stillinger-Weber.

axel.

Actually MEAM was used several times for uniaxial tension, and I want to compare the results. Langevin and NVE are working fine with MEAm for the same structure, but I am still not sure how reliable to use Langavin with NVE! Tahts why I want to use NPT but don’t know if there are parameters that i can change to overcome this wobbling behavior.

Thanks

Actually MEAM was used several times for uniaxial tension, and I want to
compare the results. Langevin and NVE are working fine with MEAm for the
same structure, but I am still not sure how reliable to use Langavin with

langevin will give you the same behavior, if you use it only as a
weakly couple thermostat. if you increase the couple (i.e. increase
friction and random noise) you're dissipating the spontaneously
released potential energy from cutting bonds better, but your dynamics
are affected.

NVE! Tahts why I want to use NPT but don't know if there are parameters that
i can change to overcome this wobbling behavior.

you need to do better/longer equilibration/relaxation of your
structure. using nve+langevin until the system is relaxed and then
switching to nve could be a way. in structures with covalent bonds
equilibration, this can be much trickier than in metals.

axel.

I didn’t understand how can I switch to NVE without controlling the temperature or you mean NVT? and How can I judge if it is weakly coupled or not? I know that these questions are not related to the software now, but it will solve so much confusion. Thanks
Mohamed

I didn’t understand how can I switch to NVE without controlling the temperature or you mean NVT? and How can I judge if it is weakly coupled or not? I know that these questions are not related to the software now, but it will solve so much confusion. Thanks

I meant what I wrote. If you don’t follow or need additional explanation, please start by carefully reading the relevant parts of the documentation and consult your adviser and/or your colleagues. It is exactly because they are not technical questions about the software, that you should make a better effort to resolve the confusion yourself. This is how actual learning happens, not by having someone tell you.

Thanks for your understanding,
Axel