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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 第一性原理 » Vasp&MedeA » 从outcar中提取原子坐标的脚本
52/1返回列表
查看: 5927  |  回复: 14
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
本帖产生 2 个 1ST强帖 ，点击这里进行查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

mika

金虫 (正式写手)

纳米宝贝

[求助] 从outcar中提取原子坐标的脚本

想从outcar里面提取每一步自洽后指定原子的xyz坐标，请问有没有相关的脚本啊，谢谢
回复此楼

» 猜你喜欢
1楼 2012-05-02 20:13:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

showtime0703

新虫 (初入文坛)
怎么样根据OUTCAR中的坐标在MS中建立模型呢
赞一下 回复此楼
14楼2014-07-13 00:01:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 15 个回答

gfei2008

木虫 (小有名气)

【答案】应助回帖

★ ★
感谢参与，应助指数 +1
WDD880227: 金币+2, 感谢您的详细解答~~很有帮助 2012-05-03 08:24:54
有个简单的指令你试试：grep "POSITION" -n OUTCAR
这个n表示你的原子数目，你把得到的数据存到一个文件下面，然后用ultraedit 去掉一些不用的提示就可以得到你想要的数据。
例如：grep "POSITION" -n OUTCAR>data
然后得到的坐标和受力在这个data文件中，然后用ultraedit处理一下就可以啦。
赞一下(1人) 回复此楼
2楼2012-05-02 23:37:23
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dxf_dxf_dxf

银虫 (小有名气)

【答案】应助回帖

★ ★ ★ ★ ★
感谢参与，应助指数 +1
fzx2008: 金币+5, 1ST强帖+1, 谢谢指教，详细 2012-05-03 13:04:02
你好！我用官方的例子来解释。
首先，vasp进行几何优化时（IBRION=2 or 3，ISIF=2，3，5等）时，才会移动原子。vasp会根据ISIF的值改变晶格常数或移动原子位置，求出体系能量最小或受力最小对应的几何结构。当设置了IBRION和ISIF参数之后，vasp实际会进行三层循环（见mannual09版7.1节Fig3）：
最内层循环，是用来矩阵对角化，求解薛定谔方程。
vasp会先给出一个初始解，带入薛定谔方程。
如果没有达到收敛条件，则修改后再次带入方程
中层循环，是用来求出自洽的电荷密度。
vasp会先猜测一个初始的电荷密度，由此构造薛定谔方程，求解得到波函数，再根据新的波函数算出（输出的）电荷密度，和（输入的）初始电荷密度对比。
如果差别没有达到收敛条件，则二者混合后作为新的输入电荷密度，重复上一步。
在mannual6.7节NWRITE-tag中，这一步对应于“electronic step”
在OUTCAR中，对应类似与“----------- Iteration    3(   2)  ----------------”中的“2”
外层循环，是用来求出能量最低、受力最小的几何结构。
vasp会根据初始POSCAR文件求出体系的总能和受力情况。根据算法改变晶格常数或离子位置，得到新的结构，再次求出总能和受力。
如果差别或绝对值未达到收敛条件，则继续改变晶格常数或离子位置。
在mannual6.7节NWRITE-tag中，这一步对应于“ion step”
在OUTCAR中，对应类似与“----------- Iteration    3(   2)  ----------------”中的“3”

如果在INCAR文件中没有设置NWRITE参数，则在OUTCAR文件中可以看到每一次ion step结束后vasp会给出优化情况：

vasp.5.2.2 15Apr09 complex
executed on             LinuxIFC date 2012.05.03  12:04:41
running on    8 nodes
distr:  one band on    1 nodes,    8 groups
--------------------------------------------------------------------------------------------------------
INCAR:
   ...
--------------------------------------------------------------------------------------------------------
...
...
----------------------------------------- Iteration    1(   1)  ---------------------------------------
    POTLOK:  cpu time    0.01: real time    0.01
    ...
----------------------------------------- Iteration    1(   2)  ---------------------------------------
     ...
--------------------------------------------------------------------------------------------------------
...
------------------------aborting loop because EDIFF is reached----------------------------------------
...
VOLUME and BASIS-vectors are now :
-----------------------------------------------------------------------------
...
direct lattice vectors reciprocal lattice vectors 0.000000000 2.750000000 2.750000000 -0.181818182 0.181818182 0.181818182 2.750000000 0.000000000 2.750000000 0.181818182 -0.181818182 0.181818182 2.750000000 2.750000000 0.000000000 0.181818182 0.181818182 -0.181818182
...
POSITION                                       TOTAL-FORCE (eV/Angst)
-----------------------------------------------------------------------------------
      4.81250      4.81250      4.81250         0.000000      0.000000      0.000000
      0.68750      0.68750      0.68750         0.000000      0.000000      0.000000
-----------------------------------------------------------------------------------
...
----------------------------------------- Iteration    2(   1)  ---------------------------------------
  ...
------------------------ aborting loop because EDIFF is reached----------------------------------------
...
VOLUME and BASIS-vectors are now :
-----------------------------------------------------------------------------
...

direct lattice vectors                 reciprocal lattice vectors
     0.000000000  2.721754370  2.721754370    -0.183705042  0.183705042  0.183705042
     2.721754370  0.000000000  2.721754370     0.183705042 -0.183705042  0.183705042
     2.721754370  2.721754370  0.000000000     0.183705042  0.183705042 -0.183705042
...
POSITION                                       TOTAL-FORCE (eV/Angst)
-----------------------------------------------------------------------------------
      4.76307      4.76307      4.76307         0.000000      0.000000      0.000000
      0.68044      0.68044      0.68044         0.000000      0.000000      0.000000
-----------------------------------------------------------------------------------

...
----------------------------------------- Iteration    3(   1)  ---------------------------------------
...
------------------------ aborting loop because EDIFF is reached ----------------------------------------
...
VOLUME and BASIS-vectors are now :
-----------------------------------------------------------------------------
...

  direct lattice vectors                 reciprocal lattice vectors
     0.000000000  2.732453619  2.732453619    -0.182985723  0.182985723  0.182985723
     2.732453619  0.000000000  2.732453619     0.182985723 -0.182985723  0.182985723
     2.732453619  2.732453619  0.000000000     0.182985723  0.182985723 -0.182985723
...
POSITION                                       TOTAL-FORCE (eV/Angst)
-----------------------------------------------------------------------------------
      4.78179      4.78179      4.78179         0.000000      0.000000      0.000000
      0.68311      0.68311      0.68311         0.000000      0.000000      0.000000
-----------------------------------------------------------------------------------
...
--------------------------------------------------------------------------------------------------------
    POTLOK:  cpu time    0.01: real time    0.01
--------------------------------------------------------------------------------------------------------

--------------------------------------------------------------------------------------------------------
reached required accuracy - stopping structural energy minimisation
...

由上可知，整个优化过程，外层循环执行了3次，共产生了四个结构：
初始结构
引用回帖:
原胞基矢
0.000000000 2.750000000 2.750000000
2.750000000 0.000000000 2.750000000
2.750000000 2.750000000 0.000000000
原子位置（相对于原胞基矢）
0.87500000 0.87500000 0.87500000
0.12500000 0.12500000 0.12500000
原子位置（相对于笛卡尔系）
4.81250000 4.81250000 4.81250000
0.68750000 0.68750000 0.68750000

第一次外层循环结束后的结构
引用回帖:
原胞基矢
     0.000000000  2.750000000  2.750000000
     2.750000000  0.000000000  2.750000000
     2.750000000  2.750000000  0.000000000
原子位置（相对于笛卡尔系）
      4.81250      4.81250      4.81250
      0.68750      0.68750      0.68750

第二次外层循环结束后的结构
引用回帖:
原胞基矢
     0.000000000  2.721754370  2.721754370
     2.721754370  0.000000000  2.721754370
     2.721754370  2.721754370  0.000000000
原子位置（相对于笛卡尔系）
     4.76307      4.76307      4.76307
      0.68044      0.68044      0.68044

第三次外层循环结束后的结构
引用回帖:
原胞基矢
     0.000000000  2.732453619  2.732453619
     2.732453619  0.000000000  2.732453619
     2.732453619  2.732453619  0.000000000
原子位置（相对于笛卡尔系）
      4.78179      4.78179      4.78179
      0.68311      0.68311      0.68311
赞一下(4人) 回复此楼

» 本帖附件资源列表

  • 欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
    本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：xiaomuchong@tal.com
  • 附件 1 : 2_4_diamondSi_vol_rex.tar.bz2
    • 2012-05-03 12:15:48, 1.39 M
3楼2012-05-03 12:59:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dxf_dxf_dxf

银虫 (小有名气)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
fzx2008: 1ST强帖+1, 谢谢指教 2012-05-03 13:34:14
uuv2010: 金币+10, 多谢分享~ 2012-05-04 18:36:54
由此给出思路及相应脚本：
  首先，要知道体系有多少个原子。这个可以从POSCAR中读取。
  然后，要告诉脚本，你关心的是哪个原子。我们把这个原子在POSCAR中的序号写到一个叫“in”的文件里面，让脚本去读就好了！
  再然后，利用grep从OUTCAR中读取我们关心的信息，并加以处理。
  依据本例，体系共有2个原子，假设我们关心POSCAR中的第二个原子。
  首先，让grep去找关键字“POSITION”，我们关心的数字在这个关键字之后的1（“---"一行）+2（2个原子）行。所以用
CODE:
grep -A3 "POSITION" OUTCAR

命令，把冗长的OUTCAR简化成如下的模样：
引用回帖:
POSITION                                       TOTAL-FORCE (eV/Angst)
-----------------------------------------------------------------------------------
      4.81250      4.81250      4.81250         0.000000      0.000000      0.000000
      0.68750      0.68750      0.68750         0.000000      0.000000      0.000000
--
POSITION                                       TOTAL-FORCE (eV/Angst)
-----------------------------------------------------------------------------------
      4.76307      4.76307      4.76307         0.000000      0.000000      0.000000
      0.68044      0.68044      0.68044         0.000000      0.000000      0.000000
--
POSITION                                       TOTAL-FORCE (eV/Angst)
-----------------------------------------------------------------------------------
      4.78179      4.78179      4.78179         0.000000      0.000000      0.000000
      0.68311      0.68311      0.68311         0.000000      0.000000      0.000000

这样还不够。我们要把”POSITON“行、”-------“行、分隔符”--“删掉：
这3种类型的行的共同特点是，如果以空格（默认的）为分隔符，则最多只有2项，第三项为空。于是用awk语句的$3>0来过滤掉它们，然后用print $1“ ”$2" "$3打印有用的前三个数字（我们不关心受力）
CODE:
grep -A3 "POSITION" OUTCAR | awk '$3>0 {print $1“  ”$2"  "$3}'

于是结果变成这个样子：
引用回帖:
4.81250  4.81250  4.81250
0.68750  0.68750  0.68750
4.76307  4.76307  4.76307
0.68044  0.68044  0.68044
4.78179  4.78179  4.78179
0.68311  0.68311  0.68311

这样还不够。我们需要的只是我们关心的那一个原子。怎么把它调出来呢？
由于我们保留了所有的原子的位置，因此到目前为止，输出的行数应该=体系原子数x外层循环次数。我们只需要找到行数对体系原子数取余数后结果为我们关心的原子对应的序数的行就行了：
NR%2==0(第2个原子，也就是余数为2，即余数为0)
CODE:
grep -A3 "POSITION" OUTCAR | awk '$3>0 {print $1"  "$2"  "$3}' | awk 'NR%2==0 {print $0}'

结果如下：
引用回帖:
0.68750  0.68750  0.68750
0.68044  0.68044  0.68044
0.68311  0.68311  0.68311

大功告成！
所以语句如下：
引用回帖:
grep -A3 "POSITION" OUTCAR | awk '$3>0 {print $1"  "$2"  "$3}' | awk 'NR%2==0 {print $0}'

当然，你需要明白“A3”中3的含义，还有"NR%2==0"中2和0的含义。否则生搬硬套会得到错误的结果~
赞一下(5人) 回复此楼
4楼2012-05-03 13:26:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 15 个回答
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 为什么资助数各大高校都创新高，自己申请怎么就这么难 +9 Kittylucky 2026-08-27 9/450 2026-08-28 01:08 by jurkat.1640
[硕博家园] 售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急 +3 XLGfUIdOmEAO 2026-08-27 3/150 2026-08-28 00:49 by 3RAFyYXYSwGS
[基金申请] 哪位高人中了，把查询到的截图贴出来让我看看，让我长长见识 +5 yuleib84 2026-08-26 6/300 2026-08-28 00:02 by yudaoqian88
[基金申请] 基金不中，共勉 +10 eulota 2026-08-26 10/500 2026-08-27 17:23 by lqllinqiaoli
[基金申请] 基金未中，这种答复是模板吗？ +5 zhaosm1982 2026-08-27 6/300 2026-08-27 16:00 by lfy8008
[基金申请] 怎么看青基中了没有啊 +5 叶九微 2026-08-26 5/250 2026-08-27 10:35 by l_zh2008
[基金申请] 2026年的国家社科基金项目通讯评审的新规则与新动向、新挑战 +7 process2012 2026-08-23 10/500 2026-08-26 19:23 by hmhminy
[基金申请] 系统查不到 +10 董八千 2026-08-26 10/500 2026-08-26 16:30 by Equinoxhua
[基金申请] 怎么查啊 +5 huang1991js 2026-08-26 5/250 2026-08-26 16:09 by Equinoxhua
[基金申请] 2026年8月25日国自然放榜前突然收到列入评审专家邮件，有关系吗？ +25 木水思豆 2026-08-25 28/1400 2026-08-26 14:53 by draco1987
[基金申请] 为什么国自然不能直接公布 +4 bjdxyxy 2026-08-26 4/200 2026-08-26 13:12 by qingmu1201
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 国合可查了 +3 paperzjh 2026-08-26 3/150 2026-08-26 10:41 by LemmonTr
[职场人生] 学生撞见辅导员送外卖，第二天全班都沉默了 +3 鲸鱼融金_浙江_� 2026-08-22 3/150 2026-08-26 09:01 by zzuzxg
[基金申请] 牛来！米来！面来！ +8 beefly 2026-08-26 8/400 2026-08-26 08:37 by xuzhipiao
[基金申请] 在坚冰还盖着北海的时候，我看到了怒放的梅花。 (金币+10) +6 ziyangfang 2026-08-25 9/450 2026-08-25 20:26 by huagongfeihu
[基金申请] 明天应该可查了！？ +6 chengyan1220 2026-08-23 6/300 2026-08-25 19:45 by zfd97
[基金申请] 人气不行了 +11 fansofjerry 2026-08-21 11/550 2026-08-25 11:04 by 孤独的英雄6
[教师之家] 跳槽后在研项目怎么办？ +5 简单化xn 2026-08-22 10/500 2026-08-23 12:38 by 简单化xn
信息提示
关闭
请填处理意见
关闭