计算化学公社

标题: QE中批量固定原子的小脚本QEfix.sh [打印本页]

作者
Author: 丁越 时间: 2022-2-8 16:31
标题: QE中批量固定原子的小脚本QEfix.sh
本帖最后由丁越于 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意思是脚本文件接受的第一个参数。

#!/bin/bash
#usage: QEfix.sh [qe input file]
fname1=`echo $1`
echo ' Input indices of the atoms to be constraint (fixed), e.g. 1,5,9-12,14-18: '
read num
echo ' Input if_pos(i), e.g. 0 0 0, which means atomic coordinates at x, y, and z direction will be fixed simutaneously:'
read if_pos
echo
pos=`grep -n ATOMIC_POSITIONS ${fname1} |cut -d ":" -f 1 |head -1`
echo $num|awk -F , '{for(i=1;i<=NF;i++) print $i}' > ${prefix}_tmp1.dat
#Generating the atoms Index into tmp2.dat
for i in $(cat ${prefix}_tmp1.dat);do
if [[ "$i" =~ ^[0-9]+$ ]];then
echo $i >> ${prefix}_tmp2.dat
else
Min=$(echo $i |awk -F \- '{print $1}')
Max=$(echo $i |awk -F \- '{print $2}')
for j in $(seq $Min $Max)
do
echo $j >> ${prefix}_tmp2.dat
done
fi
done
#fixing atoms in the QE input file
for i in $(cat ${prefix}_tmp2.dat)
do
line=$(($i + $pos))
echo $line >> ${prefix}_tmp3.dat
done
for i in $(cat ${prefix}_tmp3.dat)
do
sed -i "${i}s/$/ ${if_pos}/g" ${fname1}
done
rm -f ${prefix}_tmp*

复制代码

作者
Author: 卡开发发 时间: 2022-2-8 19:30
本帖最后由卡开发发于 2022-2-8 19:51 编辑

其实利用ase gui可以直接在ase自带的图形界面直接批量选择原子进行固定：
（1）先启动ase

$ 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的输入方式在工程角度看稍有一些缺陷，很多东西也不是太完善。

作者
Author: 丁越 时间: 2022-2-8 20:24

卡开发发发表于 2022-2-8 19:30
其实利用ase gui可以直接在ase自带的图形界面直接批量选择原子进行固定：
（1）先启动ase（2）然后界面上 ...

谢谢卡开发发老师，ASE我还没用过，主要是不会python

过段时间再学下

作者
Author: 卡开发发 时间: 2022-2-8 20:41

丁越发表于 2022-2-8 20:24
谢谢卡开发发老师，ASE我还没用过，主要是不会python过段时间再学下

python可以学学，一些轻量化的处理会轻松很多

，ase也就是其中一个生态资源而已。

作者
Author: qczgzly 时间: 2022-8-5 09:43
感谢您的脚本！似乎有个微小问题。固定参数（例如1 1 1）与原子坐标z值相连，没有空格。

N 6.690430 15.690897 3.6265201 1 1

复制代码

可以将脚本第36行作如下调整

#original
sed -i "${i}s/$/${if_pos}/g" $1
#modified
sed -i "${i}s/$/ ${if_pos}/g" $1

复制代码

即可。

作者
Author: 丁越 时间: 2022-8-6 00:11

qczgzly 发表于 2022-8-5 09:43
感谢您的脚本！似乎有个微小问题。固定参数（例如1 1 1）与原子坐标z值相连，没有空格。

可以将脚本第36 ...

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

作者
Author: yanleshan 时间: 2024-7-2 12:08
老师您好，我使用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

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

作者
Author: yanleshan 时间: 2024-7-2 13:18

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不换行，而是在同一行跟在坐标的末尾

欢迎光临计算化学公社 (http://bbs.keinsci.com/)