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

丁越

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. 
   
4. fname1=`echo $1`
   
5. echo ' Input indices of the atoms to be constraint (fixed), e.g. 1,5,9-12,14-18: '
   
6. read num
   
7. echo ' Input if_pos(i), e.g. 0 0 0, which means atomic coordinates at x, y, and z direction will be fixed simutaneously:'
   
8. read if_pos
   
9. echo
   
10. 
    
11. pos=`grep -n ATOMIC_POSITIONS ${fname1} |cut -d ":" -f 1 |head -1`
    
12. echo $num|awk -F , '{for(i=1;i<=NF;i++) print $i}' > ${prefix}_tmp1.dat
    
13. 
    
14. #Generating the atoms Index into tmp2.dat
    
15. for i in $(cat ${prefix}_tmp1.dat);do
    
16.         if [[ "$i" =~ ^[0-9]+$ ]];then
    
17.                 echo $i >> ${prefix}_tmp2.dat
    
18.         else
    
19.                 Min=$(echo $i |awk -F \- '{print $1}')
    
20.                 Max=$(echo $i |awk -F \- '{print $2}')
    
21.                 for j in $(seq $Min $Max)
    
22.                 do
    
23.                        echo $j >> ${prefix}_tmp2.dat
    
24.               done               
    
25.         fi
    
26. done
    
27. 
    
28. #fixing atoms in the QE input file
    
29. for i in $(cat ${prefix}_tmp2.dat)
    
30. do
    
31.         line=$(($i + $pos))
    
32.         echo $line >> ${prefix}_tmp3.dat
    
33. done
    
34. 
    
35. for i in $(cat ${prefix}_tmp3.dat)
    
36. do
    
37.         sed -i "${i}s/$/    ${if_pos}/g" ${fname1}
    
38. done
    
39. 
    
40. rm -f ${prefix}_tmp*
    
41. 
    
42. 
    

复制代码

yanleshan

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

yanleshan

老师您好，我使用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

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

丁越

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

可以将脚本第36 ...

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

qczgzly

感谢您的脚本！似乎有个微小问题。固定参数（例如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

复制代码

即可。

卡开发发

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

python可以学学，一些轻量化的处理会轻松很多，ase也就是其中一个生态资源而已。

丁越

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

谢谢卡开发发老师，ASE我还没用过，主要是不会python过段时间再学下

卡开发发

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