|
|
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ghcacj(金币+20):谢谢 2010-12-04 17:09:24
Lesson 2 Reax FF
就不自己翻译成中文了吧, 反正大家都能看明白的. 拼写错误啦, 语法错误啥的, 大家多包涵就是了.
如果核心内容上, 对于原理的理解上有任何错误的, 请指出.
Two tests have been carried out for NH3_Graphene_COOH_OH systems. The same initial.dat and ffielf.reax files are used. Timestep is 0.25 fs and initial velocities are generated according to Maxwell-Boltzmann distribution.
The only difference is the heating up procedure.
Test 1:
(a) Minimize the initial structure;
(b) Run 200000 steps (50 ps) to equilibrate the system, NVT ensemble with the damping parameter at 10.0;
(c) Heat up the system linearly from 500K to 2500K at the rate of 0.005 K/step, NVT ensemble with the damping rate of 10.0.
Test 2
(a) Minimize the initial structure;
(b) Run 200000 steps (50 ps) to equilibrate the system, NVT ensemble with the damping parameter at 1.0;
(c) Divide the temperature range of (500, 2500) into 5 blocks, and heat up the system linearly in each block and continuously between blocks. The rate is still 0.005 K/step, NVT ensemble with the damping rate of 1.0.
Here are the results and figures.
Figure 1 Temperature vs. simulation time, part(b), NVT ensemble.
The temperature damping parameter is 1.0 for test 1 and 10.0 for test 2. Generally, temperature can fluctuate wildly if the damping parameter is too small, while on the other hand, too large damping parameter will take a very long time to equilibrate the whole system. I may do more tests to get the optimized value for this temperature damping parameter.
Figure 1


Figure 1. Temperature vs. simulation time, the NVT ensemble at step (b).


Figure 2. ‘live’ Temperature vs. target Temperature, part (c), NVT ensemble.


Figure 3. Temperature fluctuation vs. target Temperature.
The reason for the blocks in test 2 is that, I was planning to do some NVE between each NVT iteration, hoping that could bring forth the better control over the system temperature. So in the Lammps input script, part (c) is being taken care of by the following commands:
Test 1:
fix 2 all nvt temp 500.0 2500.0 10.0
run 400000
Test 2:
variable mytemp index 500.0, 900.0, 1300.0, 1700.0, 2100.0
variable mytemp2 index 900.0, 1300.0, 1700.0, 2100.0, 2500.0
label loopa
variable i loop 5
fix 2 all nvt temp ${mytemp} ${mytemp2} 1.0
run 80000
unfix 2
next mytemp
next mytemp2
next i
jump in.reaxff loopa
There are two other issues here. One is that I realize the temperature of first step of part (c) is not exactly 500.0K. This is because the NVT of part (b) cannot bring in 500.0K as the temperature of its final output. Does it mean that part (b) is actually not necessary?
The other is that I am not sure I completely understand the philosophy of temperature programmed reaction method.
Here are my thoughts of the whole process:
Minimization is necessary to make sure you don't have any overlaps. Then when you add temperature (give your atoms velocities) your system won't blow up. Once the temperature (atomic velocity) has been assigned, you first should use a thermostat to equilibrate the system. This can be a Berendsen or a Langevin, depending on what you want to do, as these have a random term in them and this helps to equilibrate faster.
But if you want to measure dynamics at a constant temperature, you should use NVE. However, you want to alter the temperature and measure (or follow) the dynamics, so in this I would use a Nose-Hoover thermostat (fix temp/nvt).
That is why I break the whole process into three parts as I describe previously:
(1) Use minimization to minimize the structure (from my test run, it seemed that the initial structure is already pretty close to a local minimum, as the search stopped after 169 steps or so);
(2) Use a thermostat (Berendsen or Langevin) to equilibrate the temperature. You could also use a Nose-Hoover here as well and I think it would be ok. The Berendsen or Langevin thermostat would likely get you there faster.
(3) Use a thermostat to ramp up the temperature from 500 K to 2500 K.
Here I would use a Nose-Hoover thermostat since I have a pretty small system. So if I use a Berendsen or a Langevin thermostat, I would have to use a large damping parameter here and the results likely would not be realistic. If I use a small damping parameter, the results would be better for the dynamics but the temperature control would not be very good. So a Nose-Hoover with a reasonable damping parameter is probably the best. That is what I use in part (c):
fix 2 all nvt temp 500.0 2500.0 10.0
Of course I then need to find that reasonable damping parameter.
If we could manually do several iterations at the very beginning of part (C), let’s see what happens there:
(a) Image 1, t=0, T(target)=500K, we assign the atomic velocities at 500K. With atomic positions and force field, we know the forces on each atom and thus the acceleration. Do one iteration by Verlet algorithm to go to step 2.
(b) Image 2, t= 0.25fs, T(target)=500.005K, from (a) we get the velocities at this step 2, new coordinates, new forces and new accelerations. The real temperature can be calculated at this step and rescaling (of atomic velocities) is applied by Nose-Hoover algorithm to bring the real temperature to the target temperature 500.005K. Use the rescaled velocities to go a step further to image 3.
(c) Keep doing this time integration and controlling over temperature to do the temperature ramp up steadily along the way.
By increasing temperature slowly and steadily, each atom possess more total energy (potential and kinetic), and make it possible to overcome the energy barriers such as bonds, angles, and even Pi-Pi interaction, and thus old-bond breaking and new-bond forming is more and more possible along the way. Also, from this point of view, there is no advantage of dividing the temperature range into several block.
THE MOST SERIOUS PROBLEM is:
From Figure 3, we can see large fluctuation at high temperature. Is that acceptable? If not, how can we bring down the fluctuation? |
|