【分享】尝试lammps中, 分享中...
Lesson 1
Loop inside Lammps input script
想要尝试在lammps的input script里面做循环, 结果因为一个小问题, 捣鼓了我几个小时, 这才完全通过测试. 分享一下.
(1) input script里面的循环块
include tempfile
include tempfile2
variable i loop 10
label loopa
fix 2 all nvt temp ${mytemp} ${mytemp2} 100.0
run 200
unfix 2
next mytemp
next mytemp2
next i
jump SELF loopa
这个的SELF是让程序执行到这里, 跳回自己, 然后从标签 loopa开始执行. 当然loopa是随便取的, 你可以用CHN来做label.
另外, 在某些情况下, 如果要让c++ rewind, 那么最好在执行脚本里面这样写:
lmp -in script
而不是
lmp < script
否则, 你人品不好的时候, 会出问题, 哈哈.
(2) tempfile 和 tempfile2是在同目录下的另外两个文件. tempfile的文件内容是:
variable mytemp index 500.0 700.0 900.0 1100.0 1300.0 1500.0 1700.0 1900.0 2100.0 2300.0
注意, 只有一行! 至于lammps能读多长的一行, 我还没有测试.
tempfile2的文件内容, 也是一行:
variable mytemp2 index 700.0 900.0 1100.0 1300.0 1500.0 1700.0 1900.0 2100.0 2300.0 2500.0
折腾我的问题是, 我原先在这两个tempfile文件中, 数值之间用逗号分隔, 但是事实上, 是需要用空格分隔的.
(3) 如果你需要做的循环不是很多, 那么不需要额外准备tempfile 和 tempfile2文件. 而只是需要在 input script中这样做:
variable mytemp index 500.0 700.0 900.0 1100.0 1300.0 1500.0 1700.0 1900.0 2100.0 2300.0
#variable mytemp2 index 700.0 900.0 1100.0 1300.0 1500.0 1700.0 1900.0 2100.0 2300.0 2500.0
variable i loop 10
label loopa
fix 2 all nvt temp ${mytemp} ${mytemp2} 100.0
run 200
unfix 2
next mytemp
next mytemp2
next i
jump SELF loopa
当然, 这里不需要原先的这两句 include语句了.
总结一下你需要熟悉的命令:
variable, include, jume, next
End of Lesson 1.
ENJOY.
[ Last edited by qphll on 2010-12-3 at 11:54 ]
返回小木虫查看更多
京公网安备 11010802022153号
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?
太好了，正在学习中
I am sharing some other test output with you.
reax FF, Lammps, Temperature Programmed Reax MD
(1) The NVT damping parameter of Lammps
The tests were carried out for 10 different values of that damping coefficients. I find that changing that parameter could not bring forth better control of temperature during the whole damping process. The flucutation is reduced a little bit at some value. But generally I would say there is no pointing of doing this type of parameter optimizaiton.
(2) Heating rate
With a same timestep of 0.25 fs, I tried the different heating rate from 500K to 2500K, namely, 0.01 K/step, 0.05K/step and 1 K/step. From the attached temperature flucutation plots, I observed (a) changing the rate from 0.005 K/step to 0.01 K/step or 0.05 K/step does not change much of the temperature flucutation. Does is mean that we could heat up the system faster? Or we do have very good reason in using the heating rate of 0.005 K/step? (b) at lower rate (0.005 K/step, 0.01 K/step and 0.05 K/step), it is not easy to control the temperature flucutation. But at 1 K/step, compared with the other three rates, the fluctuation looks smooth for [900K, 1800K]. It may comes from fewer statistic points, but is it worthy of trying different heating rates for low temperature, median temperature and high temperature?
0.01 K/step
0.05 K/step
1 K/step
I do not think those tests are going to change anything about the TPRD procedures. They are just my thoughts when I am collecting the outputs.
Q: 如何在跑完reax MD以后, 判断是否有反应发生了呢?
A:
(1) 在input script文件中使用命令 'fix reax/bonds' 输出成键信息, 具体请查询手册;
(2) 使用这里给出的小程序来读取output文件, 将提取出来的信息作图判断.
请根据自己的体系, 自行修改给出的小程序. 我还没来得及改进....
!# DEC.9, 2010
!# QPHLL
!#
!# This is a program to read the output from 'fix reax/bond', TPRD, Lammps
!# The output is saved into file "bonds.connect", where each image is divided
!# into three parts:
!#
!# (1) Head, 7 Lines;
!# (2) Body, No._of_atom Lines;
!# (3) Tail, 1 Line
!#
!# The total number of images is related with the output frequence and number of iterations.
!# In this case, it is "number of iteration+1".
!#
!# Each line in Body part is made up of the following parameters:
!# id, type, nb, id_1, id_2, ... id_nb, mol, bo_1, bo_2, ... bo_nb, abo, nlp, q
!# abo = atomic bond order
!# nlp = number of lone pairs
!# q = atomic charge
!#
!# PLEASE DOUBLE CHECK YOUR OWN LAMMPS INPUT SCRIPT & OUTPUT AND MAKE CORRESPONDING CHSNGES
program main
implicit none
integer I, J, K, L
integer image, natom
integer headline, tailline
integer id, atype, nb, bd1, bd2, bd3, mol
double precision bo1, bo2, bo3, abo, nlp, q
open (unit=10, file='bonds.connect')
open (unit=20, file='N129.txt', status='unknown')
open (unit=21, file='N133.txt', status='unknown')
open (unit=22, file='N137.txt', status='unknown')
open (unit=23, file='N141.txt', status='unknown')
open (unit=24, file='N145.txt', status='unknown')
open (unit=25, file='N149.txt', status='unknown')
open (unit=26, file='N153.txt', status='unknown')
open (unit=27, file='N157.txt', status='unknown')
image = 2000
headline = 7
tailline = 1
natom = 160
do I = 1, image+1
do J = 1, headline
read(10,*)
end do
do K = 1, natom
read(10,*) id, atype, nb
write(*,*) id, atype, nb
if (atype .eq. 4) then
backspace 10
read(10,*) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
write(*,*) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
if (id .eq. 129) then
write(20, 200) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
elseif (id .eq. 133) then
write(21, 200) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
elseif (id .eq. 137) then
write(22, 200) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
elseif (id .eq. 141) then
write(23, 200) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
elseif (id .eq. 145) then
write(24, 200) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
elseif (id .eq. 149) then
write(25, 200) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
elseif (id .eq. 153) then
write(26, 200) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
elseif (id .eq. 157) then
write(27, 200) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
200 format(7I4, 6f14.3)
endif
endif
enddo
do L =1,tailline
read(10,*)
enddo
enddo
end program main
更新一下, 稍微修改了一下前楼的程序. 现在如果在计算中有发生过反应, 那么分析程序会将那部分信息输出到另外的一个文件中.
!# DEC.9, 2010
!# QPHLL
!# www.muchong.com
!#
!# This is a program to read the output from 'fix reax/bond', TPRD, Lammps
!# The output is saved into the file "bonds.connect", where each image is divided
!# into three parts:
!#
!# (1) Head, 7 Lines;
!# (2) Body, No._of_atom Lines;
!# (3) Tail, 1 Line
!#
!# The total number of images is related with the output frequence and number of iterations.
!# In this case, it is "number of iteration+1".
!#
!# Each line in Body part is made up of the following parameters:
!# id, type, nb, id_1, id_2, ... id_nb, mol, bo_1, bo_2, ... bo_nb, abo, nlp, q
!# abo = atomic bond order
!# nlp = number of lone pairs
!# q = atomic charge
!#
!# PLEASE DOUBLE CHECK YOUR OWN LAMMPS INPUT SCRIPT & OUTPUT AND MAKE CORRESPONDING CHSNGES
program main
implicit none
integer I, J, K, L
integer image, natom
integer headline, tailline
integer id, atype, nb, bd1, bd2, bd3, bd4, mol
double precision bo1, bo2, bo3, bo4, abo, nlp, q
open (unit=10, file='bonds.connect')
open (unit=20, file='N129.txt', status='unknown')
open (unit=21, file='N133.txt', status='unknown')
open (unit=22, file='N137.txt', status='unknown')
open (unit=23, file='N141.txt', status='unknown')
open (unit=24, file='N145.txt', status='unknown')
open (unit=25, file='N149.txt', status='unknown')
open (unit=26, file='N153.txt', status='unknown')
open (unit=27, file='N157.txt', status='unknown')
open (unit=30, file='reactionRecord.txt', status='unknown')
!# Make changes accordingly.
image = 2000
headline = 7
tailline = 1
natom = 160
do I = 1, image+1
! Skip the head part
do J = 1, headline
read(10,*)
end do
! Each image has 'natom' lines
do K = 1, natom
! read in the first three number each line to determine:
! (1) what type of atom it is, atype
! the correspondence in Lammps: 1-C, 2-H, 3-O, 4-N, 5-S
! (2) how many bonds it has, nb
! this 'nb' determines the following bond_link information & bond_order paramaters of the same line
read(10,*) id, atype, nb
! TEST
! write(*,*) id, atype, nb
if (atype .eq. 4) then
backspace 10
! Should have some easier way to replace this "IF", I am just toooo lazy.
! Thanks to the fact that the maximum number of bonds is 4. ^-^
!??? is it possible that nb = 0 ??? KEEP THAT IN MIND.
if (nb.eq.0) then
read(10,*) id, atype, nb, mol, abo, nlp, q
if (id .eq. 129) then
write(20, 200) id, atype, nb, mol, abo, nlp, q
elseif (id .eq. 133) then
write(21, 200) id, atype, nb, mol, abo, nlp, q
elseif (id .eq. 137) then
write(22, 200) id, atype, nb, mol, abo, nlp, q
elseif (id .eq. 141) then
write(23, 200) id, atype, nb, mol, abo, nlp, q
elseif (id .eq. 145) then
write(24, 200) id, atype, nb, mol, abo, nlp, q
elseif (id .eq. 149) then
write(25, 200) id, atype, nb, mol, abo, nlp, q
elseif (id .eq. 153) then
write(26, 200) id, atype, nb, mol, abo, nlp, q
elseif (id .eq. 157) then
write(27, 200) id, atype, nb, mol, abo, nlp, q
200 format(4I4, 3f14.3)
endif
! If bd .ne. 3, it measn reaction is happening to Nitrogen atom.
write (30, 300) I, id, atype, nb, mol, abo, nlp, q
300 format(5I4, 3f14.3)
elseif (nb.eq.1) then
read(10,*) id, atype, nb, bd1, mol, bo1, abo, nlp, q
if (id .eq. 129) then
write(20, 201) id, atype, nb, bd1, mol, bo1, abo, nlp, q
elseif (id .eq. 133) then
write(21, 201) id, atype, nb, bd1, mol, bo1, abo, nlp, q
elseif (id .eq. 137) then
write(22, 201) id, atype, nb, bd1, mol, bo1, abo, nlp, q
elseif (id .eq. 141) then
write(23, 201) id, atype, nb, bd1, mol, bo1, abo, nlp, q
elseif (id .eq. 145) then
write(24, 201) id, atype, nb, bd1, mol, bo1, abo, nlp, q
elseif (id .eq. 149) then
write(25, 201) id, atype, nb, bd1, mol, bo1, abo, nlp, q
elseif (id .eq. 153) then
write(26, 201) id, atype, nb, bd1, mol, bo1, abo, nlp, q
elseif (id .eq. 157) then
write(27, 201) id, atype, nb, bd1, mol, bo1, abo, nlp, q
201 format(5I4, 4f14.3)
endif
! If bd .ne. 3, it measn reaction is happening to Nitrogen atom.
write (30, 301) I, id, atype, nb, bd1, mol, bo1, abo, nlp, q
301 format(6I4, 4f14.3)
elseif (nb.eq.2) then
read(10,*) id, atype, nb, bd1, bd2, mol, bo1, bo2, abo, nlp, q
if (id .eq. 129) then
write(20, 202) id, atype, nb, bd1, bd2, mol, bo1, bo2, abo, nlp, q
elseif (id .eq. 133) then
write(21, 202) id, atype, nb, bd1, bd2, mol, bo1, bo2, abo, nlp, q
elseif (id .eq. 137) then
write(22, 202) id, atype, nb, bd1, bd2, mol, bo1, bo2, abo, nlp, q
elseif (id .eq. 141) then
write(23, 202) id, atype, nb, bd1, bd2, mol, bo1, bo2, abo, nlp, q
elseif (id .eq. 145) then
write(24, 202) id, atype, nb, bd1, bd2, mol, bo1, bo2, abo, nlp, q
elseif (id .eq. 149) then
write(25, 202) id, atype, nb, bd1, bd2, mol, bo1, bo2, abo, nlp, q
elseif (id .eq. 153) then
write(26, 202) id, atype, nb, bd1, bd2, mol, bo1, bo2, abo, nlp, q
elseif (id .eq. 157) then
write(27, 202) id, atype, nb, bd1, bd2, mol, bo1, bo2, abo, nlp, q
202 format(6I4, 5f14.3)
endif
! If bd .ne. 3, it measn reaction is happening to Nitrogen atom.
write (30, 302) I, id, atype, nb, bd1, bd2, mol, bo1, bo2, abo, nlp, q
302 format(7I4, 5f14.3)
elseif (nb.eq.3) then
read(10,*) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
if (id .eq. 129) then
write(20, 203) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
elseif (id .eq. 133) then
write(21, 203) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
elseif (id .eq. 137) then
write(22, 203) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
elseif (id .eq. 141) then
write(23, 203) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
elseif (id .eq. 145) then
write(24, 203) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
elseif (id .eq. 149) then
write(25, 203) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
elseif (id .eq. 153) then
write(26, 203) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bd3, abo, nlp, q
elseif (id .eq. 157) then
write(27, 203) id, atype, nb, bd1, bd2, bd3, mol, bo1, bo2, bo3, abo, nlp, q
203 format(7I4, 6f14.3)
endif
elseif (nb.eq.4) then
read(10,*) id, atype, nb, bd1, bd2, bd3, bd4, mol, bo1, bo2, bo3, bo4, abo, nlp, q
if (id .eq. 129) then
write(20, 204) id, atype, nb, bd1, bd2, bd3, bd4, mol, bo1, bo2, bo3, bo4, abo, nlp, q
elseif (id .eq. 133) then
write(21, 204) id, atype, nb, bd1, bd2, bd3, bd4, mol, bo1, bo2, bo3, bo4, abo, nlp, q
elseif (id .eq. 137) then
write(22, 204) id, atype, nb, bd1, bd2, bd3, bd4, mol, bo1, bo2, bo3, bo4, abo, nlp, q
elseif (id .eq. 141) then
write(23, 204) id, atype, nb, bd1, bd2, bd3, bd4, mol, bo1, bo2, bo3, bo4, abo, nlp, q
elseif (id .eq. 145) then
write(24, 204) id, atype, nb, bd1, bd2, bd3, bd4, mol, bo1, bo2, bo3, bo4, abo, nlp, q
elseif (id .eq. 149) then
write(25, 204) id, atype, nb, bd1, bd2, bd3, bd4, mol, bo1, bo2, bo3, bo4, abo, nlp, q
elseif (id .eq. 153) then
write(26, 204) id, atype, nb, bd1, bd2, bd3, bd4, mol, bo1, bo2, bd3, bo4, abo, nlp, q
elseif (id .eq. 157) then
write(27, 204) id, atype, nb, bd1, bd2, bd3, bd4, mol, bo1, bo2, bo3, bo4, abo, nlp, q
204 format(8I4, 7f14.3)
endif
! If bd .ne. 3, it measn reaction is happening to Nitrogen atom.
write (30, 304) I, id, atype, nb, bd1, bd2, bd3, bd4, mol, bo1, bo2, bo3, bo4, abo, nlp, q
304 format(9I4, 7f14.3)
! Corresponding to "if (nb.eq.0) then "
endif
! Corresponding to "if (atype .eq. 4) then"
endif
enddo
do L =1,tailline
read(10,*)
enddo
enddo
end program main，
前面的两个版本, 还是不满意.
Lammps作者给的c代码,也是不满意.
自己写, 想了一天了, 也还是不满意.
忽然发现现在的主要的问题是: 我想从计算的结果中分析什么?!
加油, 加油.
贴出Lammps 的Adian给的c代码, 有需要的, 自行下载.
http://g.zhubajie.com/urllink.php?id=10358303ba0v26t5ccef82ff
楼主真有耐心，一次传这么多，，