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

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 分子模拟 » LAMMPS/DL_POLY » 【分享】尝试lammps中, 分享中...
161/2返回列表12下一页
查看: 5769  |  回复: 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的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

qphll

金虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+10):谢谢 2010-12-10 10:41:06
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
赞一下(2人) 回复此楼
Life, Love, Laugh.
5楼2010-12-10 01:47:12
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

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的回帖

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的回帖

qphll

金虫 (正式写手)
★ ★ ★ ★ ★
zh1987hs(金币+5):谢谢 2010-12-11 23:02:52
内容已删除
回复此楼
Life, Love, Laugh.
7楼2010-12-11 10:32:02
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hbcsucy

木虫 (小有名气)

AAA
楼主真有耐心，一次传这么多，，
赞一下 回复此楼
day day
8楼2010-12-13 16:11:36
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+10, 模拟EPI+1): 谢谢 2011-02-24 13:30:03
好久没有做什么更新, 这两个月来在Lammps, 尤其是reaxff上面花了不少时间, 写了不少小程序, 现在基本到一个尾声. 很多程序写的时候一鼓作气, 没有很好地加入备注, 今天花了一整天的时间注释了一些程序的用途. 有可能会慢慢发出来.

这是才写下的一段话, 是对于lammps+reax中最重要两个输出文件的处理意见. 程序本身对别人没有什么帮助, 但是这段备注我想和大家分享.


! The "REAXFRAG" code needs two files:
! (1) coordinate table <---- "molfra.configure" <----- from Lammps "molfra" dump file, prepared by this code
! (2) bond order table <---- "bonds.connect"    <---- directly from Lammps 'fix' dump
!
! The corresponding commamds in "in.reaxff" file look like:
!
! "dump    1  all custom 10 molfra id type q x y z vx vy vz mass"
! "dump_modify 1 sort id"   <---- This "sort" option is very important.
!
! "fix 2 all reax/bonds 10 bonds.connect"
!
! Some further explanations to the "molfra" and "bonds.connect" files.
!
! Each frame of "molfra" is made up of the head part, 9 line and body part, N line (N=number of atoms of the system).
! So an output of "m" MD iterations will produce (m+1)*(N+9) lines in "molfra".
!
! Each frame of "bonds.connect" file has head part, 7 line, body part, N line (N=number of atoms of the system)
! and tail part, 1 line, which contains a "#" sign as the ending remark.
! So an output of "m" MD iterations will produce (m+1)*(N+8) lines in "molfra".
!
! "m+1" comes from the fact that Lammps dose output the t=0 image.
!
!!!!!!!!!!!!
! IT IS VERY IMPORTANT that you use same output frequency for all the dumps in Lammps.
! You also need to 'sort' the output so that atoms in different dump files are in good match.
! That explains why I prefer to seprate the minimization part from MD runs.
!!!!!!!!!!!!
!
! The other input is the number of frames in 'molfra', which is the number of iterations
! you request in "in.reaxff" file, such as the following:
!
! "run 1200000"
!
!
! Required input:
! "molfra" from Lammps & number of frames in the molfra file
!
! Please always double check to make sure that is what you want to have.
! @ QPHLL, Feb.23,2010
赞一下(1人) 回复此楼
Life, Love, Laugh.
9楼2011-02-24 07:13:34
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+20): 谢谢 2011-02-24 13:30:14
这是另外一个script, 就是用来格式转化. 没有什么, 但是往往需要很多这样的script在各个软件之间作为桥梁.

有谁需要, 尽管拿去吧.


! This is to write the output XYZ file from Lammps output format into:
! (1) the normal XYZ format which is readable by most progrms, such as imol
! (2) the "initial.dat" file which is the required input when you use "read_data"
!     in the "in.reaxff" file
!
! You can use this script to 'look' at each image along the way of your simulation
! or convert at any point of Lammps XYZ file into a new "initial.dat" file.
! I find it useful since I seperate minimization from MD production, and I use this
! script if there is a need for rerun a job.
!
! Be aware of the two facts:
! (1) the IDs of atoms are defined in "initial.dat" file of Lammps
! You have to compare the definition of atom with "in.reaxff" file to make sure
! that you've got the right atom types;
!
! (2) the box size is hard-coded here as:
!
! lx = 39.352194
! ly = 42.600000
! lz = 60.000000
! halflz = lz/2.0
!
! Pay attention to the PBC issue and wrap your coordinates back in one smooth image.
!
! You need to change accordingly or treat it as an interactive input.
!
! I have a different version(lmp2lmp_ver.f90), where basically I match-up
! the atom types between "initial.dat" and "in.reaxff".

! The other version(batchlmp2lmp.f90) is capable of batch conversion:
! (1) generated an index of files that you want to convert
! This can be done by this command at your target folder
! # find *.xyz > list
! (2) the format of that index file is as the following:
! Line one: Number of files (N)
! Line two to N: the name of file
!
! Required input: name of XYZ file and number of atom type
!
! Please always double check to make sure that is what you want to have.
! @ QPHLL, Feb.23,2010
!

program lmp2lmp
implicit none
integer natom, ntype,I, id
integer zero, one,two,three,four,five
double precision x, y, z,lx,ly,lz, halflz
character*1 C, H, O, N, S, atype*2, ctmp1*2, ctmp2*2,ctmp3*2
character*100 filename, temp, temp2

C='C'
H='H'
O='O'
N='N'
S='S'

zero=0
one=1
two=2
three=3
four=4
five=5

lx = 39.352194
ly = 42.600000
lz = 60.000000
halflz = lz/2.0


write (6,*)'Output from LAMMPS MINIMIZATION (probably structure.*.xyz) ?'
read (5,'(a40)')filename

temp='initial.dat'//filename
temp2 ='formatted'//filename

open(unit=30,file=filename,status='old')
open(unit=20, file=temp, status='unknown')
open(unit=40, file=temp2, status='unknown')


read(30,*) natom
read(30,*)
write (6,*)'How many types of atoms ?'
read (5,*) ntype
write(40,*) natom
write(40,*)

    write(20,*)
    write(20,*)
    write(20,*) natom, "atoms"
    write(20,*)
    write(20,*) ntype,"atom types"
    write(20,*)
    write(20,*)   0.00000,lx, "xlo xhi"
    write(20,*)   0.00000,ly, "ylo yhi"
    write(20,*)   0.00000,lz, "zlo zhi"
    write(20,*)
    write(20,*)
    write(20,*) "Atoms"
    write(20,*)
   
   

do I = 1, natom

read(30,*)id, x, y, z

if(z.lt.halflz) then
write(20,100)I,id,zero,x, y, z
else
write(20,100)I,id,zero,x, y, z-lz
endif

if(z.lt.halflz) then

if (id.eq.1) then
write(40,*)C, x, y, z
elseif(id.eq.2) then
write(40,*)H, x, y, z
elseif(id.eq.3) then
write(40,*)O, x, y, z
elseif(id.eq.4) then
write(40,*)N, x, y, z
elseif(id.eq.5) then
write(40,*)S, x, y, z
else
write(6,*)'Wrong atom ID!'
stop
endif

else

if (id.eq.1) then
write(40,*)C, x, y, z-lz
elseif(id.eq.2) then
write(40,*)H, x, y, z-lz
elseif(id.eq.3) then
write(40,*)O, x, y, z-lz
elseif(id.eq.4) then
write(40,*)N, x, y, z-lz
elseif(id.eq.5) then
write(40,*)S, x, y, z-lz
else
write(6,*)'Wrong atom ID!'
stop
endif

endif


100 format(I6,1x,I3,1x,I3,1x,3f12.6)

enddo

end program lmp2lmp
赞一下(1人) 回复此楼
Life, Love, Laugh.
10楼2011-02-24 07:16:55
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 qphll 的主题更新
161/2返回列表12下一页
普通表情 龙 兔 虎 猫 高级回复(可上传附件)
信息提示
关闭
请填处理意见
关闭