计算化学公社

 找回密码 Forget password
 注册 Register
Views: 5247|回复 Reply: 7
打印 Print 上一主题 Last thread 下一主题 Next thread

[Quantum ESPRESSO] QE中批量固定原子的小脚本QEfix.sh

[复制链接 Copy URL]

432

帖子

11

威望

3422

eV
积分
4074

Level 6 (一方通行)

本帖最后由 丁越 于 2022-8-6 00:19 编辑

QE中批量固定原子的小脚本QEfix.sh


      虽然在pwgui中提供了固定原子的选项,但是一个一个地固定原子坐标显得颇为麻烦,遂写了个小脚本可以很方便的固定原子坐标。脚本的使用很简单,例如将rutile的结构载入gview中,然后选择所要固定的原子后,在tools -> atom selection中就可以复制所要固定原子的序号了。接下来,通过pwgui生成如结构优化输入文件rutile.relax.in(不需要固定原子),然后执行 ./QEfix.sh rutile.relax.in 后,脚本首先会提示需要固定原子的序号,粘贴回车;紧接着会提示如何固定,如输入 0 0 0就表示将x、y、z三个方向同时固定。
      接下来说说脚本的思路:原子的序号样式是这样的 2,3,4-10,因此,首先我们需要将这些序号拆解成一个一个的单个序号形式(2,3,4,5....10),然后写个for循环在批量去在每个需要固定的原子后面添加如0 0 0。拆解原子序号时,我们先用awk以逗号为分割符将其拆开,然后存进第一个临时文件tmp1.dat;接下来,我们需要判断tmp1.dat中数据假如是单个数字的,那就直接把它存进临时文件tmp2.dat中,假如是以4-10这样的数据,我们设置Min=4,Max=10,然后利用for循环把Min~Max中的数据缀加到tmp2中,此时,整个原子序号就被拆解成了单个序号存进了tmp2.dat中。接下来,由于输入文件中第一个原子所在行数不是1,那么我们还需要将tmp2.dat数据都加上行数得到tmp3.dat; 最后,再次利用for循环按照tmp3.dat中原子所在行数去添加如0 0 0。脚本中,[ "$i" =~ ^[0-9]+$ ]是比较常用的判断字符串包含关系的语法。sed -i "${i}s/$/${if_pos}/g" $1中$表示行尾,s表示替换,g为global,$1意思是脚本文件接受的第一个参数。

  1. #!/bin/bash
  2. #usage: QEfix.sh [qe input file]

  3. fname1=`echo $1`
  4. echo ' Input indices of the atoms to be constraint (fixed), e.g. 1,5,9-12,14-18: '
  5. read num
  6. echo ' Input if_pos(i), e.g. 0 0 0, which means atomic coordinates at x, y, and z direction will be fixed simutaneously:'
  7. read if_pos
  8. echo

  9. pos=`grep -n ATOMIC_POSITIONS ${fname1} |cut -d ":" -f 1 |head -1`
  10. echo $num|awk -F , '{for(i=1;i<=NF;i++) print $i}' > ${prefix}_tmp1.dat

  11. #Generating the atoms Index into tmp2.dat
  12. for i in $(cat ${prefix}_tmp1.dat);do
  13.         if [[ "$i" =~ ^[0-9]+$ ]];then
  14.                 echo $i >> ${prefix}_tmp2.dat
  15.         else
  16.                 Min=$(echo $i |awk -F \- '{print $1}')
  17.                 Max=$(echo $i |awk -F \- '{print $2}')
  18.                 for j in $(seq $Min $Max)
  19.                 do
  20.                        echo $j >> ${prefix}_tmp2.dat
  21.               done               
  22.         fi
  23. done

  24. #fixing atoms in the QE input file
  25. for i in $(cat ${prefix}_tmp2.dat)
  26. do
  27.         line=$(($i + $pos))
  28.         echo $line >> ${prefix}_tmp3.dat
  29. done

  30. for i in $(cat ${prefix}_tmp3.dat)
  31. do
  32.         sed -i "${i}s/$/    ${if_pos}/g" ${fname1}
  33. done

  34. rm -f ${prefix}_tmp*


复制代码






QEfix.sh

1021 Bytes, 下载次数 Times of downloads: 13

评分 Rate

参与人数
Participants 2
威望 +1 eV +1 收起 理由
Reason
1049714804 + 1
sobereva + 1

查看全部评分 View all ratings

自由发挥,野蛮生长

3621

帖子

3

威望

1万

eV
积分
18429

Level 6 (一方通行)

第一原理惨品小作坊

2#
发表于 Post on 2022-2-8 19:30:12 | 只看该作者 Only view this author
本帖最后由 卡开发发 于 2022-2-8 19:51 编辑

其实利用ase gui可以直接在ase自带的图形界面直接批量选择原子进行固定:
(1)先启动ase
  1. $ ase gui CeO2.xsd
复制代码
(2)然后界面上工具栏的Tools->Constraints调出固定的菜单。
(3)此时可以通过鼠标框选原子(被选中的原子的轮廓会加粗).
(4)再到Constraints呼出的菜单点击Fix按钮,此时选中的原子上会显示X记号表面这些原子被固定了。
(5)保存,文件名设置为*.pwi就行,ase能识别为qe的in文件(小心不要用poscar之类的命名,否则会被误认为vasp的结构)。

其实磁矩也可以使用ase的gui设置,大家可以自己试试。最后,即便导出的结果不能满意,其实也可以按照ase的json格式存储,再简单写一些程序进一步处理。事实上在上述基础上之前我用PyQT设计出一版额外增加QE参数的图形界面,目前不方便公开,因为一方面是和别人合作的项目,另一方面就是QE的输入方式在工程角度看稍有一些缺陷,很多东西也不是太完善。


评分 Rate

参与人数
Participants 1
eV +3 收起 理由
Reason
hebrewsnabla + 3

查看全部评分 View all ratings

日常打哑谜&&探寻更多可能。
原理问题不公开讨论,非商业性质讨论欢迎私聊。
本周忙

432

帖子

11

威望

3422

eV
积分
4074

Level 6 (一方通行)

3#
 楼主 Author| 发表于 Post on 2022-2-8 20:24:26 | 只看该作者 Only view this author
卡开发发 发表于 2022-2-8 19:30
其实利用ase gui可以直接在ase自带的图形界面直接批量选择原子进行固定:
(1)先启动ase(2)然后界面上 ...

谢谢卡开发发老师,ASE我还没用过,主要是不会python过段时间再学下
自由发挥,野蛮生长

3621

帖子

3

威望

1万

eV
积分
18429

Level 6 (一方通行)

第一原理惨品小作坊

4#
发表于 Post on 2022-2-8 20:41:14 | 只看该作者 Only view this author
丁越 发表于 2022-2-8 20:24
谢谢卡开发发老师,ASE我还没用过,主要是不会python过段时间再学下

python可以学学,一些轻量化的处理会轻松很多,ase也就是其中一个生态资源而已。
日常打哑谜&&探寻更多可能。
原理问题不公开讨论,非商业性质讨论欢迎私聊。
本周忙

111

帖子

0

威望

4239

eV
积分
4350

Level 6 (一方通行)

5#
发表于 Post on 2022-8-5 09:43:02 | 只看该作者 Only view this author
感谢您的脚本!似乎有个微小问题。固定参数(例如1 1 1)与原子坐标z值相连,没有空格。
  1. N     6.690430   15.690897    3.6265201 1 1
复制代码

可以将脚本第36行作如下调整
  1.        #original
  2.         sed -i "${i}s/$/${if_pos}/g" $1
  3.        #modified
  4.         sed -i "${i}s/$/ ${if_pos}/g" $1
复制代码

即可。
gcq#413: "I know poetry is not dead, nor genius lost; nor has Mammon gained power over either, to bind or slay; they will both assert their existence, their presence, their liberty and strength again one day." (Jane Eyre in Jane Eyre by Charlotte Bronte)

432

帖子

11

威望

3422

eV
积分
4074

Level 6 (一方通行)

6#
 楼主 Author| 发表于 Post on 2022-8-6 00:11:34 | 只看该作者 Only view this author
qczgzly 发表于 2022-8-5 09:43
感谢您的脚本!似乎有个微小问题。固定参数(例如1 1 1)与原子坐标z值相连,没有空格。

可以将脚本第36 ...

谢谢提醒,我刚刚再写QEtoolkit2.0版本的时候也发现了这个问题,估计明后天就发出来了,用它生成QE的输入文件将会非常方便。
自由发挥,野蛮生长

5

帖子

0

威望

53

eV
积分
58

Level 2 能力者

7#
发表于 Post on 2024-7-2 12:08:53 | 只看该作者 Only view this author
老师您好,我使用QEfix.sh时报错如下:

grep: 101-353.relax.in: Is a directory
./QEfix.sh: line 31: 1 + : syntax error: operand expected (error token is "+ ")
cat: _tmp3.dat: No such file or directory

请问如何解决,本人纯小白实在想不到办法,盼解答

5

帖子

0

威望

53

eV
积分
58

Level 2 能力者

8#
发表于 Post on 2024-7-2 13:18:24 | 只看该作者 Only view this author
yanleshan 发表于 2024-7-2 12:08
老师您好,我使用QEfix.sh时报错如下:

grep: 101-353.relax.in: Is a directory

找到错误原因了,这是由于我使用的GUI是BURAI 3.1而非pwGUI,BURAI生成输入文件为一个文件夹,.sh读取不了,需要单拎出来文件夹里的geom.in文件作为脚本的输入文件,随之而来的问题是,运行QEfix生成的新geom.in文件,每个原子坐标后跟的0 0 0处在换行后的下一行,故BURAI读取不了正确的坐标

我目前只会手动把0 0 0合并到同一行,希望老师可以教教我如何修改脚本使得输出时0 0 0不换行,而是在同一行跟在坐标的末尾

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2024-11-24 06:07 , Processed in 0.194003 second(s), 25 queries , Gzip On.

快速回复 返回顶部 返回列表 Return to list