应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 分子模拟 » LAMMPS/DL_POLY » 【分享】尝试lammps中, 分享中...
51/1返回列表
查看: 5772  |  回复: 15
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
本帖产生 2 个 模拟EPI ，点击这里进行查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

qphll

金虫 (正式写手)

[交流] 【分享】尝试lammps中, 分享中...已有9人参与

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 ]
回复此楼

» 收录本帖的淘帖专辑推荐

资源收集 材料计算模拟实用技巧 分子模拟 MD分子动力学
我学习计算的一些帖子 分子动力学 量化 LAMMPS
关于Lammps

» 本帖已获得的红花（最新10朵）

fly771125 生与死后

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
Life, Love, Laugh.
1楼2010-12-03 11:43:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+10, 模拟EPI+1):谢谢 2010-12-10 10:41:17
更新一下, 稍微修改了一下前楼的程序. 现在如果在计算中有发生过反应, 那么分析程序会将那部分信息输出到另外的一个文件中.



!# 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
赞一下(1人) 回复此楼
Life, Love, Laugh.
6楼2010-12-10 03:29:50
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 16 个回答

qphll

金虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
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?
赞一下(1人) 回复此楼
Life, Love, Laugh.
2楼2010-12-04 15:07:09
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

sg18408926

至尊木虫 (著名写手)
太好了，正在学习中
赞一下 回复此楼
3楼2010-12-04 17:27:08
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+20):谢谢 2010-12-10 10:40:58
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.
赞一下(1人) 回复此楼
Life, Love, Laugh.
4楼2010-12-07 23:19:38
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 16 个回答
普通表情 龙 兔 虎 猫 高级回复(可上传附件)
信息提示
关闭
请填处理意见
关闭