24小时热门版块排行榜    

查看: 1808  |  回复: 3
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

cgzhang_gg

铜虫 (著名写手)

[求助] pair_style coul/cut 的测试

有谁测试过pair_style   coul/cut 的对势,  我通过pair_writes 导出的势,与我fortran按照库伦作用表达式算出的势相差很大 ,定性上都是错的。也就是说,通过pair_writes 导出的 TI_Ti ,Ti_O, O_O 三个对势都是正的,而且与电荷量没有关系。我怀疑是我in.文件里那个地方设置的不对。下面是我的in.文件的内容:
大家帮我看看哪里有不对的地方,或者咱们讨论下,好人有好报!

注:我测试的体系TiO2金红石结构:一个原胞里面有两个Ti和四个O.,晶格常数是4.493A。
设置第一类原子Ti的电荷为2.192e, e为元电荷,一个质子所具有的电量。相应的O的电荷设置为:-1.098e.
单位是按照metal 来确定的。

    # TiO2 coul/cut detect right or wrong
    units metal
    atom_style charge
    boundary p p p

    lattice    custom 4.493    a1   1.0  0.0  0.0    a2   0.0   1.0   0.0    a3   0.0   0.0   0.669&
                                           basis   0.0   0.0   0.0    basis   0.5   0.5   0.5 &
                                            basis  0.303   0.303   0.0   basis   0.697   0.697   0.0 &
                                            basis  0.803   0.197   0.5    basis   0.197    0.803    0.5
    region                 mybox    block   0 1   0 1   0 1
    create_box          2       mybox
    create_atoms      2       box basis  1  1  basis  2  1 &
                                                basis  3  2  basis  4  2  basis  5  2  basis  6   2
    mass                   1  47.867 # 设置第一类原子Ti的质量:47.867
    mass                   2  15.9994 # 设置第二类原子O的质量:15.9994

    set                      type 1    charge   2.196 # 设置第一类原子Ti 的电荷为2.196
    set                      type 2    charge   -1.098 # 设置第二类原子O的电荷为-1.098

    # 1-coul/cut 8.5
    pair_style coul/cut 8.5
    pair_coeff 1 1
    pair_coeff 1 2
    pair_coeff 2 2

     pair_write      1    1     50000    r    0.4     8.0     table.T_T.c    Ti_Ti
     pair_write      1    2     50000    r    0.4     8.0     table.T_O.c   Ti_O
     pair_write      2    2     50000    r    0.4     8.0     table.O_O.c   O_O
  
     quit
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cgzhang_gg

铜虫 (著名写手)

送红花一朵
引用回帖:
3楼: Originally posted by cgzhang_gg at 2013-06-18 09:55:54
经过测试,可以很明确的说,用pair_write 导出库伦对势的时候没有考虑电荷项,只是导出了:C/r 这一项, r 按照所用单位值确定,C是 energy-conversion constant  。如units  metal   r 的单位为 A ,C的值为14.39 ...

本贴已完结,请不用回复了!
4楼2013-06-18 09:57:10
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 4 个回答

cgzhang_gg

铜虫 (著名写手)

我再补充下,我把用pair_write 导出来的对势,T_T,T_O,O_O的三个对势与 1/r 相比较,结果发现比值恰为:14.399645,这正好是energy-conversion constant C 见lammps手册 pair_style coul/cut command 大约872页, 然后我又测试了通过read_data 命令来读取构型以及每个原子的带电量,然后又通过pair_write 导出对势,结果还是不对?我现在怀疑lammps 里面有一个bug, 就是通过pair_write 命令导出对势的时候根本没考虑qi 和qj,至于程序内部计算力,以便更新位置和速度,应该考虑了,要不然会出现很多离奇的错误。我测试的in.文件如下:
# TiO2 TiO2.pair.table detect right or wrong
units          metal
atom_style     charge
#atom_style     atomic
boundary       p p p

read_data      data.TiO2

mass           1 47.867   # Ti
mass           2 15.9994  # O

# 1-coul/cut 8.5
pair_style coul/cut 8.5
#pair_coeff * *
pair_coeff 1 1
pair_coeff 1 2
pair_coeff 2 2

pair_write 1 1 50000 r 0.4 9.0 table.T_T.c Ti_Ti
pair_write 1 2 50000 r 0.4 9.0 table.T_O.c Ti_O
pair_write 2 2 50000 r 0.4 9.0 table.O_O.c  O_O

fix 1 all nve
run 10
quit

data.TiO2 的内容如下:
#TiO2

6       atoms

0       bonds
0       angles
0       dihedrals
0       impropers

2       atom types
0       bond types
0       angle types
0       dihedral types
0       improper types

0 4.493   xlo xhi
0 4.493   ylo yhi
0 3.00582 zlo zhi

Atoms

1 1 2.196 0 0 0
2 1 2.196 2.2465 2.2465 1.50291
3 2 -1.098 1.36138 1.36138 0
4 2 -1.098 3.13162 3.13162 0
5 2 -1.098 3.60788 0.885121 1.50291
6 2 -1.098 0.885121 3.60788 1.50291

在这里为了简单,我只取了一个原胞的原子。即,两个Ti,一个O.  Atoms 每一列的意义如下:
id type q x y z

» 本帖已获得的红花(最新10朵)

2楼2013-06-17 22:49:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cgzhang_gg

铜虫 (著名写手)

送红花一朵
引用回帖:
2楼: Originally posted by cgzhang_gg at 2013-06-17 22:49:16
我再补充下,我把用pair_write 导出来的对势,T_T,T_O,O_O的三个对势与 1/r 相比较,结果发现比值恰为:14.399645,这正好是energy-conversion constant C 见lammps手册 pair_style coul/cut command 大约872页,  ...

经过测试,可以很明确的说,用pair_write 导出库伦对势的时候没有考虑电荷项,只是导出了:C/r 这一项, r 按照所用单位值确定,C是 energy-conversion constant  。如units  metal   r 的单位为 A ,C的值为14.399645 eV.A,   (C/r)的单位为eV 。 但是在程序内部计算时,是考虑电荷的。
    atom_style   charge
。。。。。。。。。。。。。
    set                      type 1    charge   2.196 # 设置第一类原子Ti 的电荷为2.196
    set                      type 2    charge   -1.098 # 设置第二类原子O的电荷为-1.098
。。。。。。。。。。。。。
这样设置来模拟coul作用是完全正确的。

» 本帖已获得的红花(最新10朵)

3楼2013-06-18 09:55:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 建议基金发布提前给出明确的时间点 +10 kulium 2026-08-21 13/650 2026-08-21 21:58 by alongwaytogo
[基金申请] 看来今天不会放榜了? +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
[基金申请] 时间戳又变了 +13 wuchongjun 2026-08-20 19/950 2026-08-21 17:21 by 紫杉醇
[基金申请] filecode,4个jtjc了 +11 ziyangfang 2026-08-19 14/700 2026-08-21 16:51 by applepetty
[基金申请] 让我中一个面上吧! +11 大萍1987 2026-08-20 13/650 2026-08-21 14:48 by 20120902066
[基金申请] 人气不行了 +7 fansofjerry 2026-08-21 7/350 2026-08-21 14:13 by Vivilian
[基金申请] 感觉是下周放榜了 +7 angus9576 2026-08-17 12/600 2026-08-21 13:38 by weiyin
[基金申请] 我面上完蛋了 +7 且听虎啸 2026-08-20 8/400 2026-08-21 12:31 by 酷酷墨镜
[基金申请] 什么时候开奖? +7 CrisMessi 2026-08-18 7/350 2026-08-21 12:15 by weiyin
[基金申请] 今天基金会出结果吗?20260819 +15 kkkl_v 2026-08-19 16/800 2026-08-21 11:22 by 365372687
[基金申请] 应该是下周三26日公布了吧? +4 哈哈蛤? 2026-08-21 4/200 2026-08-21 10:58 by Vivilian
[论文投稿] 投稿咨询 +5 wwm09 2026-08-17 7/350 2026-08-21 10:11 by 期刊论文帮手
[基金申请] 估计是周四 +3 archvillain 2026-08-18 3/150 2026-08-21 01:48 by jnhyjjm
[基金申请] 今天系统多次维护,明天很可能放榜! +10 zju2000 2026-08-16 11/550 2026-08-20 20:02 by cl479861084
[基金申请] 时间戳变了,能看出什么问题? +18 基诺咪客 2026-08-17 23/1150 2026-08-20 17:19 by Godzela
[基金申请] 重要消息,中午系统在维护 +11 yuleib84 2026-08-18 12/600 2026-08-20 11:09 by xskun
[教师之家] 为什么余额宝的年化利率越来越低?主要原因有哪些? +5 瞬息宇宙 2026-08-15 5/250 2026-08-19 20:23 by super2002521
[基金申请] 今天放榜吗? +14 布布和一二 2026-08-19 15/750 2026-08-19 18:07 by gltch
[基金申请] 2027广东省杰青 +4 奶牛小黑 2026-08-15 10/500 2026-08-19 11:02 by wanfengnew
[基金申请] 明天放榜? +5 Shxjjxjkx 2026-08-18 5/250 2026-08-18 18:14 by -大大大大大-
信息提示
请填处理意见