Si Nanowire

Hello LAMMPS users,

I am in the process of thermalizing a Si nanowire. The lattice is a “diamond” lattice. I use periodic boundary conditions in all dimensions. The axis is aligned along x and due to the periodicity; it is thus “an infinitely long wire” (and it is surrounded by vacuum; I make the box large enough to zero image effects). The wire is only 2 nm long with about the same diameter (so very tiny, around 220 atoms). I am using a Tersoff potential as the force field.

The goal is to get this wire first in an NVT ensemble and then in an NVE ensemble.
My approach is: Thermalize -> NVT -> NVE

From reading other posts on this forum I found this method of thermalization (I have been working on this for a couple of months now so I have tried a lot of permutations of this, I just can’t seem to get it right):
To get the wire at 300 K:

  1. Minimize the structure first:
    min_style sd
  2. Zero linear and angular velocities (so the structure should not pick up any of these later on?)
    velocity all zero linear
    velocity all zero angular
  3. Use a Langevin thermostat starting at T=20K (some low T) going to T=300K in 200000 steps. The damping I use is 0.06 (I explain why later). I use “zero yes” to make sure the wire does not drift.
    fix NVE all nve
    fix LANGEVIN all langevin 20 300 0.06 31415926 zero yes
  4. Use an NPT ensemble to zero the pressure along the x direction (I don’t put a barostat on the y and z directions since this is vacuum) with a 0.01 T damping and 10 times larger P damping. I tried varying these damping parameters and found no noticeable effect.
    fix NPT all npt temp 300 300 0.01 x 0.0 0.0 0.1
  5. Use a Langevin thermostat for 70 million steps (when I observe the energy evolution versus step it really seems to take this long! “long time” is relative, but comparing to papers in the literature it seems like an extremely long time) (step size = 0.5 fs). (again with “zero yes”.
    fix NVE all nve
    fix LANGEVIN all langevin 300 300 0.06 31415926 zero yes
    At this point I assume we are in an NVT ensemble (I also tried following this Langevin step with an actual “nvt” ensemble for millions of steps, the following result is the same).

Now I consider the structure is thermalized and in an NVT ensemble. Then, I want it in an NVE ensemble.

  1. Do a pure NVE run for 8 million steps.
    fix NVE all nve

Before my actual questions: Why the 0.06 damping in Langevin?
This number seems very low but of course “low” is relative and depends on the system. I would like the lowest possible damping as my understanding is that this will get me to equilibrium fastest?
I tried values around 0.06 as well. When I observe the total energy it seems to behave like an exponential decay (this should be fine), though decaying very slowly meaning: It takes on the order of 70 million steps(!) before this decay reaches a plateau. There are some issues (of physical origin) with the surface atoms in such a wire: they have many dangling bonds, but even with this I am not sure 70 million steps is okay. But, my understanding is that larger damping will affect the structure less (less viscosity between the “background fluid” we imagine created by the thermostat and the atoms in the structure). However, and this is why I don’t do it, it also takes longer to thermalize as far as I understand. I would like to spend much less than 70 million steps.

In the NVE ensemble I want to compute the heat flux along x.

What is the problem?
Well, 70 million thermalization steps seems like too much. And from comparing my result to what I know it should give I am getting something close, but still too far, from what is expected. This makes me worry about the process I use.

Now my questions are:

  1. Do you have any suggestions as to how I figure out whether my Langevin damping is okay (I tried a bunch of values and the structure seems to stick together and the energy evolves in this exponential, but slow, fashion for values from 0.01 to 10; I do see a large dependence on the damping parameter in general: the property I am computing in the NVE ensemble can jump an order of magnitude when changing the damping coefficient although the thermalizations look very similar)?

  2. In the NVE ensemble my structure picks up a large angular momentum. I don’t like not being in control of the structure since I expect it to stand still (overall). I don’t want to use “momentum zero” to get rid of it since this changes the physics of the velocities which I need for the heat flux calculation, but maybe I can use this during thermalization and then stop doing it during the NVE run? I would just assume it should not pick up such a rotation since it starts with zero angular momentum? Does anybody know how to get rid of this? I can think of reducing the time step but I am already pretty low compared to other values in the literature (0.55 and 0.8, mine is 0.5) for very similar (if not the same) problems (I can try anyway I know I cannot reason based on other peoples work like that; however I am not sure the rotation is of main concern physically).

  3. Does anybody have familiarity with this problem and can direct me to a (maybe) better method not requiring this many steps or just give some pointers as to how I should proceed?

From a paper I found I know that this is suppose to work (although I am the first to do it in LAMMPS; the paper uses another MD program).

I appreciate your time and suggestions.

Thank you,

Hi Jesper,
I have worked previously on similar stuff. However, I did not use langevin thermostat for equilibration. I just used NVT followed by the NPT. My nanowires were quite long then yours but if you want to dig more into that, here is the link to the paper. Please have a look at supporting information as well.!divAbstract

Best Regards,

Sorry, forgot to mention, after NPT, I did run NVE to make sure energies are conserved and there is no angular momentum.

Thank you for your response Vikas.

  • Do you mind telling me how many steps you used for the NVT and NPT runs (order of magnitude)?
  • Did you minimize before the run?

Of course, as the wire becomes larger I (eventually) approach the bulk situation which I have no problems working in and reproducing; as you point out yourself your Si wire seems very large compared to mine.

How did your total energy evolve during the NVT run, did you see that “exponential decay” behavior, see attached figure (ignore labels etc. it is just the trend of the graph I am showing)?
That is actually for 90 million steps (just to make a very long run) which I think is way to much!
I will try your thermalization approach.

Again, thank you,

I forgot to ask:

  • How did you get your system to 300 K?
    Did you just use “nvt” with start and stop temperature of 300K or did you slowly heat it up like I am trying to do?
  • Did you use “velocity create” or something similar before the nvt call to start the atoms with some velocities?

You help is much appreciated,

Dear Jesper,

Out of top of my head, this is what I did

  1. Minimize the structure

  2. Velocity create (~600 K) and start running NVE first. This lead to initial velocity of 300 K (I think half of the energy went in potential energy; Check it to cofirm) which started increasing and saturated at ~340 K (after ~25 ns). I know it is long but I had ~80000 atoms (lots of unrelaxed surface area in the beginning) and I wanted to make sure that temperature does saturate before I do anything else.

  3. Then, ramped down the temperature from 340 to 300 K in NPT ensemble over 1 ns.

  4. Then, ran NPT for ~ 3-4 ns at 300 K (To adjust the nanowire length and make sure there are no residual stresses left in nanowire axial direction due to periodicity)

  5. Finally, I ran NVE to make sure everything is equilibrated (tempeature is not increasing, potnetial energy is equilibrated and fluctuating around some mean value)

Here are my graphs (in the supporting information from the paper). Have a look at supplementary information A.

Let me know if anything is not clear. For you, since the nanowire is not long, thing should equilibrate much quicker.

c2nr30602f.pdf (712 KB)

Dear Jesper,
Out of top of my head, this is what I did
1. Minimize the structure
2. Velocity create (~600 K) and start running NVE first. This lead to
initial velocity of 300 K (I think half of the energy went in potential
energy; Check it to cofirm) which started increasing and saturated at ~340 K
(after ~25 ns). I know it is long but I had ~80000 atoms (lots of unrelaxed
surface area in the beginning) and I wanted to make sure that temperature
does saturate before I do anything else.

this is the step, where it would be beneficial to use fix langevin,
since that will accelerate the process of equipartitioning of kinetic
energy. in my tests, it typically works the most effective to
initially start with a comparatively short relaxation time for fix
langevin and then two thirds down the line make it much longer, to
also have a smooth transitioning to nvt or npt ensemble. using plain
fix nve has the big risk that any low frequency soundwaves that may
result from the initial high potential energy structure quickly
relaxing do may still be present and taint any analysis. fix langevin
will quickly dissipate them (so would multiple quenches with
minimization, but that is more tedious to set up).

the slow ramping of temperatures is really only needed, if you have a
very delicate structure, e.g. biopolymers, membranes and other complex
systems that require a complex relaxation and equilibration workflow.