octopus 快速上手

参考
简单运行
inp文件
结构优化
- 关于续算的一些测试
- 尝试解决报错循环报错
  - 报错iteration is not making progress towards solution
  - 续算和交替优化
性质计算
- 输入
- Output参数
- 输出文件
  - static
- TDOutput
搭建说明网页
示例
- 周期性示例
- td 计算示例
报错

octopus-8.4快速上手记录
短期内应该不会用octopus了，记录免忘记
不得不说octopus官方的教程写的真是太好了

参考

简单运行

octopus有原子坐标既可以运行
如，新建inp填入

%Coordinates
'H'  |0 | 0 | 0
%

然后

mpirun -np 核数 octopus

就可以了，octopsu会提供赝势等等参数

运算之后，生成exec restart static
static里面包含计算结果和输出文件

inp文件

文件格式

可以在inp中自定义变量，变量名区分大小写，并进行相关数学计算Mathematical expressions
计算参数名不区分大小写，不区分顺序
系统变量也不区分大小写,如Output = dos 等价于Output = DOS
注释#开头
可以用 include 文件名 把其他文件放到inp文件中
输入文件只能是inp,不能改变名字
说明书中的参数解释:
- DOSComputePDOS ;参数名
- Section: Output ;分类,没用
- Type: logical ;类型:float,logical,integer,string,block等, logical类型赋值为0/1
- Default: false ; 默认参数
不同版本的参数会出现不同,本文主要基于Octopus-8/10.4
逻辑变量赋值为yes/no(或1/0),不是true/T之类的

计算参数

MaximumIter=10 最大自洽步数,用于快速debug代码

赝势

自己下载的赝势,UPF的赝势虽然可以用,但是计算量特别大

%Species
'Ag' | species_pseudo | file |  './Ag_ONCV_PBE_fr.upf'
'C' | species_pseudo | file |  './C_ONCV_PBE_fr.upf'
'H' | species_pseudo | file |  './H_ONCV_PBE_fr.upf'
%
#XML (QSO, UPF, and PSML, PSP8) pseudopotential support is an experimental feature.
#需要设置ExperimentalFeatures=yes
ExperimentalFeatures=yes

octopus目录的赝势

如sg15则读取octopus-10.4/share/octopus/pseudopotentials/quantum-simulation.org/sg15/目录的赝势,
由于是upf格式,还要加上ExperimentalFeatures=yes

standard的赝势只有有限的几个H, Li, C, N, O, Na, Si, S, Ti, Se, Cd.

%Species
'Ag' | species_pseudo | set | sg15
#set后面支持的参数,说明书中的PseudopotentialSet由可用列表
%
ExperimentalFeatures=yes
#如果有很多元素则可以这样设置所有都默认是sg15
PseudopotentialSet=sg15

泛函

有的赝势,octopus无法识别泛函的类型,这时会用默认的LDA进行计算,这时需要指定泛函

**   Unknown XC functional for species 'Ag'

如

#exchange:x   correlation:c
XCFunctional=gga_x_pbe+gga_c_pbe

原子质量

默认根据元素名识别相对原子质量,不用指定也可以. 如果元素名随意设置,会报错Cannot determine the element for species CC

vdw修正

VDWCorrection=vdw_d3
ExperimentalFeatures=yes

控制参数

`BoxShape` 计算网格形状

BoxShape

默认是minimum,原子越远离原点,网格点越多

使用sphere时,如果原子/电子超出设置的球形边界,不会被计算,会出现误差

This variable decides the shape of the simulation box. The default is minimum for finite systems and parallelepiped for periodic systems. Note that some incompatibilities apply:

Spherical or minimum mesh is not allowed for periodic systems. Cylindrical mesh is not allowed for systems that are periodic in more than one dimension. box_image is only allowed in 2D.

Options:

sphere: The simulation box will be a sphere of radius Radius. (In 2D, this is a circle.), 半径radius = 15.0*angstrom
cylinder: The simulation box will be a cylinder with radius Radius and height (in the x-direction) of 2 Xlength.
minimum: The simulation box will be constructed by adding spheres created around each atom (or user-defined potential), of radius Radius.
parallelepiped: The simulation box will be a parallelepiped whose dimensions are taken from the variable Lsize.
box_image: The simulation box will be defined through an image, specified with BoxShapeImage. White (RGB = 255,255,255) means that the point is contained in the simulation box, while any other color means that the point is out. The image will be scaled to fit Lsize, while its resolution will define the default Spacing. The actual box may be slightly larger than Lsize to ensure one grid point = one pixel for default Spacing.
hypercube: (experimental) The simulation box will be a hypercube or hyperparallelepiped. This is equivalent to the parallelepiped box but it can work with an arbitrary number of dimensions.
user_defined: The shape of the simulation box will be read from the variable BoxShapeUsDef.

BoxShape = sphere
spacing = 0.19*angstrom
radius = 15.0*angstrom

`TDTimeStep`时间步长

TDTimeStep默认单位是$3.674932539796e-02 \hbar/eV$, $1 \hbar/eV = 0.658211951 fs$

TDTimeStep = 0.002 # 4.83776903362502e-05 fs, *826827 ~ 40fs

`KPointsUseTimeReversal`时间反演对称性

动力学参数

是否移动原子,只有打开MoveIons,原子位置后面的yes和no才会生效

MoveIons = yes
%Coordinates
'C' | angstrom*   -0.0142124300 | angstrom*   -1.1141926050 | angstrom*    4.7610890850  | no
'H' | angstrom*   -0.0135859800 | angstrom*   -1.1175805350 | angstrom*    3.6428618450  | yes
'H' | angstrom*   -0.9162065400 | angstrom*   -1.6344720150 | angstrom*    5.1074087650  | yes
'H' | angstrom*    0.8873471600 | angstrom*   -1.6333341650 | angstrom*    5.1089549050  | yes
'H' | angstrom*   -0.0152745800 | angstrom*   -0.0720769150 | angstrom*    5.1038265250  | yes
%

续算

RestartWriteInterval=50 #每隔多少步保存一次文件,用于续算
FromScratch=no #no读入保存文件,进行续算. yes从头算

终止程序

在运行目录,建一个stop文件

** Warning:
**   Clean STOP

  27658      55.316000 -20002.140099         1         3.865
Info: Writing output to output_iter/td.0027658

结构优化

先用其他程序(VASP)优化可以减少很多时间
cg_,steep等算法,优化到极小值点时程序会崩溃掉,这时只能换其他算法如fire
steep,cg_到极小值后换cg_,steep还是极小值,只能换fire 交替优化时GOMaxIter=50就可以, 提高精度继续重新开始
先使用运动速度快的,如GOMethod=fire配GOFireMass=0.01.GOFireMass越小越大,0.01比默认的1快很多倍
fire算法是通过跑MD的方式,也可以设置一个大的GOMaxIter=500,GOFireMass=0.01,苛刻的优化参数GOTolerance = 0.01*eV/angstrom,让程序跑一段时间
fire不存在到达极小值无法优化奔溃掉
GOObjective的取值,影响不大

关于续算的一些测试

使用steep算法,续算后力的变化也比较有限
steep续算到无法续算后的结构,steep,bfgs算法输入次结构,也是直接无法计算
使用fire继续算这个结构,倒是可以又深层次优化一下

CalculationMode = go
GOMethod = cg_bfgs #先快速steep,最后用cg_bfgs,然后逐渐提高各项精度
GOTolerance = 0.01*eV/angstrom #Default: 0.001 H/b (0.051 eV/A)
#yes移动原子,no固定原子
%Coordinates
'Ag' | angstrom*   -0.0013139800 | angstrom*   -2.7930235850 | angstrom*   -2.7300387650 | no
'Ag' | angstrom*   -0.0013139800 | angstrom*   -2.7930235850 | angstrom*   1.7300387650 | yes
%

更多固定方式

#0动,1固定
%GOConstrains
  'C' | 1 | 0 | 0
  'O' |  1 | 0 | 0
%

输出

work-geom.xyz 每一步的结构
last.xyz 最后一步结构.当设置续算FromScratch = no时,会读入last.xyz进行续算
geom/ 日志,力,位移. 每进行一次最小化循环MINIMIZATION ITER,就输出一个geom/go.nnnn.xyz文件
min.xyz 优化结束时的输出

跑完了,不一定是收敛了,可能是内存溢出程序终止了,结束输出

iter   18 : etot -2.19005958E+02 : abs_dens 3.71E-08 : force  4.73E-07 : etime     0.0s
ETA: .......1......2.......3......4......5.......6......7.......8......9......0

iter   19 : etot -2.19005957E+02 : abs_dens 2.30E-08 : force  4.92E-07 : etime     0.0s
#这里的force 是前后两步force的差别,用来判断scf是否收敛的.同static/convergence
             Info: Writing states. 2022/09/14 at 15:16:35


        Info: Finished writing states. 2022/09/14 at 15:16:35

Info: SCF converged in   19 iterations



++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++ MINIMIZATION ITER #:    1 ++++++++++++++++++++++
  Energy    = -219.0059565865 eV
  Max force =    0.0486830383 eV/A
  Max dr    =    0.0013009175 A
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



Writing final coordinates to min.xyz
Info: Finished writing information to 'restart/gs'.

             Calculation ended on 2022/09/14 at 15:16:35

                          Walltime:  09.651s

Octopus emitted 5 warnings.

尝试解决报错循环报错

报错`iteration is not making progress towards solution`

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++ MINIMIZATION ITER #:    7 ++++++++++++++++++++++
  Energy    = -219.0067014605 eV
  Max force =    0.0153721653 eV/A
  Max dr    =    0.0000000000 A
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++




**************************** FATAL ERROR *****************************
*** Fatal Error (description follows)
*--------------------------------------------------------------------
* From node =    0
*--------------------------------------------------------------------
* Error occurred during the GSL minimization procedure:
* iteration is not making progress towards solution
**********************************************************************

Abort(999) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 999) - process 0

多核运行时根本看不到报错信息啊….,不同服务器都报错tag=1620299

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++ MINIMIZATION ITER #:    7 ++++++++++++++++++++++
  Energy    = -219.0066996657 eV
  Max force =    0.0153715180 eV/A
  Max dr    =    0.0000000000 A
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



Abort(336193796) on node 7 (rank 7 in comm 0): Fatal error in PMPI_Send_init: Invalid tag, error stack:
PMPI_Send_init(152): MPI_Send_init(buf=0x1c55854, count=1, MPI_INTEGER, dest=0, tag=1620299, MPI_COMM_WORLD, request=0x7fffb210d1f8) failed
PMPI_Send_init(103): Invalid tag, value is 1620299

看起来就是,两个构型之间的force一样,原子相对位移为0,Max dr = 0.0000000000 A优化不动了,官方给的解决方案就是用go.xxxx.xyz重算

This means that, although the minimization algorithm managed to get very close (~ 0.014 eV/Å), it was unable to improve the solutions up to the desired accuracy. Unfortunately there are not many things that can be done in this case. The most simple one is just to restart the calculation starting from the last geometry (use the last go.XXXX.xyz file) and see if the code is able to improve further the geometry.

续算和交替优化

有时候系统的续算,一重算就错,自己写的续算更稳

除了构型的输入参数input.inp
初始构型structure.xyz

function xyz2oct()
{
OUTPUTFILE=inp

if [ -d geom ] #存在目录,说明算过了,开始续算
then
  #geom会被octopus清空
  #检查是否存在正常的xyz文件,否则退出
  INPUTFILE=geom/go.0001.xyz
  if [ ! -f $INPUTFILE  ]; then echo "Error: No xyz in geom";exit;fi
  if [ $(wc -l $INPUTFILE | awk '{ print $1 }') -lt 3 ]; then echo $INPUTFILE error;exit;fi
  if [ $(wc -l $INPUTFILE | awk '{ print $1 }') -lt $(($(head -1 $INPUTFILE)+2)) ] ; then echo $INPUTFILE error;exit;fi
  #
  #选择最新正常的xyz文件
  INPUTFILE=$(ls geom/go*.xyz | sort | tail -1)
  if [ $(wc -l $INPUTFILE | awk '{ print $1 }') -lt 3 ]; then echo $INPUTFILE error;INPUTFILE=$(ls geom/go*.xyz | sort | tail -2 | head -1);fi
  if [ $(wc -l $INPUTFILE | awk '{ print $1 }') -lt 3 ]; then echo $INPUTFILE error;INPUTFILE=$(ls geom/go*.xyz | sort | tail -2 | head -1);fi
else #第一次开始计算,使用structure作为结构
  INPUTFILE=structure.xyz
fi
#最后确定xyz文件是否正常
if [ -f  $INPUTFILE ]
then
echo "-------------------------------"
ls -l $INPUTFILE
echo "build input from "$INPUTFILE
cat $INPUTFILE
echo "-------------------------------"
#检查结构,失败直接退出
if [ $(wc -l $INPUTFILE | awk '{ print $1 }') -lt 3 ]; then echo $INPUTFILE error;exit;fi
if [ $(wc -l $INPUTFILE | awk '{ print $1 }') -lt $(($(head -1 $INPUTFILE)+2)) ] ; then echo $INPUTFILE error;exit;fi
#生成输入文件
cp input.inp $OUTPUTFILE
startline=3
x=0;y=0;z=0
#x=$(tail -n +$startline $INPUTFILE | awk 'BEGIN {max = -9999; min=9999} {if ($2+0 > max) max=$2; if ($2+0 < min) min=$2 } END {printf "%.10f" , 0.5*(max+min) }')
#y=$(tail -n +$startline $INPUTFILE | awk 'BEGIN {max = -9999; min=9999} {if ($3+0 > max) max=$3; if ($3+0 < min) min=$3 } END {printf "%.10f" , 0.5*(max+min) }')
#z=$(tail -n +$startline $INPUTFILE | awk 'BEGIN {max = -9999; min=9999} {if ($4+0 > max) max=$4; if ($4+0 < min) min=$4 } END {printf "%.10f" , 0.5*(max+min) }')
echo "%Coordinates" >> $OUTPUTFILE
tail -n +$startline $INPUTFILE | awk '{printf "'\''%s'\'' | angstrom* %15.10f | angstrom* %15.10f | angstrom* %15.10f \n",$1,$2-'$x',$3-'$y',$4-'$z'}'  >> $OUTPUTFILE
#tail -n +$startline $INPUTFILE | awk '{printf "'\''%s'\'' |  %15.10f |  %15.10f |  %15.10f \n",$1,$2-'$x',$3-'$y',$4-'$z'}'  >> $OUTPUTFILE
echo "%" >> $OUTPUTFILE
else
echo "Can't find "$INPUTFILE "breakdown"
exit
fi
}

rm -rf stop1
xyz2oct
mpirun octopus | tee result

cp work-geom.xyz work-geom.bak.xyz
cp result result.bak
echo geom.log > geom.log 
#cg算法找不到极小值会报错不输出min.xyz
#fire虽然会输出min.xyz但是可能达到了最大步数,这时换成cg算法看可不可以
while  (($(grep 'min.xyz' result | wc -l ) < 1)) || [ $(grep 'GOMethod' inp | awk '{ print $3 }') == fire ] ;
do
echo CNQ $(date) >> result
echo "---------------" >> geom.log
for i in $(ls geom); do echo $i >> geom.log; cat geom/$i >> geom.log ; done

#收集上一次计算用的算法
GOMethod=fire
if [ $(grep 'GOMethod' inp | awk '{ print $3 }') == fire ];     then GOMethod=cg_bfgs ; fi
if [ $(grep 'GOMethod' inp | awk '{ print $3 }') == steep ];    then GOMethod=fire ; fi
if [ $(grep 'GOMethod' inp | awk '{ print $3 }') == cg_bfgs ];  then GOMethod=fire ; fi
if [ $(grep 'GOMethod' inp | awk '{ print $3 }') == cg_bfgs2 ]; then GOMethod=fire ; fi
if [ $(grep 'GOMethod' inp | awk '{ print $3 }') == cg_pr ];    then GOMethod=fire ; fi
if [ $(grep 'GOMethod' inp | awk '{ print $3 }') == cg_fr ];    then GOMethod=fire ; fi

xyz2oct
#>[可选]准备更换演化算法
echo "Change GOMethod = $GOMethod"; sed -i "s/GOMethod/GOMethod = $GOMethod #/g" inp
mpirun  octopus > result #报错和其他日志信息在slurm的输出
#可以通过修改参数,然后touch stop换算法、参数续算
cat work-geom.xyz >> work-geom.bak.xyz
cat result >> result.bak

if [ -f stop1 ]; then exit;fi
done

echo "Relax done"

如果最后还没有优化出来,则把geom.log中最后一步的结构保存到geom/go.0001.xyz再输入继续跑,spacing等其他计算参数也决定了能收敛的最高精度利用这种方式,优化的CH4键长和vasp误差0.001ang,还可以

#算之前清空result,文件存在与否决定是否是完成了第一次计算
if [ -f result ]; then mv result result.bak;fi
rm -rf stop1 stop

function xyz2oct()
{
OUTPUTFILE=inp

if [ -f result ] #存在result文件,说明算过了,开始续算
then
  #有时候cg算法失败了,找不到最优解,也输出不了xyz,此时换其他cg然后不改变输入文件中的结构参数
  #下面都采用return.终止则替换为exit
  #
  INPUTFILE=geom/go.0001.xyz
  if [ ! -f $INPUTFILE  ]; then echo "Error: No xyz in geom";return;fi
  if [ $(wc -l $INPUTFILE | awk '{ print $1 }') -lt 3 ]; then echo $INPUTFILE error;return;fi
  if [ $(wc -l $INPUTFILE | awk '{ print $1 }') -lt $(($(head -1 $INPUTFILE)+2)) ] ; then echo $INPUTFILE error;return;fi
  #
  #选择最新正常的xyz文件
  INPUTFILE=$(ls geom/go*.xyz | sort | tail -1)
  if [ $(wc -l $INPUTFILE | awk '{ print $1 }') -lt 3 ]; then echo $INPUTFILE error;INPUTFILE=$(ls geom/go*.xyz | sort | tail -2 | head -1);fi
  if [ $(wc -l $INPUTFILE | awk '{ print $1 }') -lt 3 ]; then echo $INPUTFILE error;INPUTFILE=$(ls geom/go*.xyz | sort | tail -2 | head -1);fi
else #第一次开始计算,使用structure作为结构
  INPUTFILE=structure.xyz
fi
#最后确定xyz文件是否正常
if [ -f  $INPUTFILE ]
then
echo "-------------------------------"
ls -l $INPUTFILE
echo "build input from "$INPUTFILE
cat $INPUTFILE
echo "-------------------------------"
#检查结构,失败直接退出
if [ $(wc -l $INPUTFILE | awk '{ print $1 }') -lt 3 ]; then echo $INPUTFILE error;exit;fi
if [ $(wc -l $INPUTFILE | awk '{ print $1 }') -lt $(($(head -1 $INPUTFILE)+2)) ] ; then echo $INPUTFILE error;exit;fi
#生成输入文件
cp input.inp $OUTPUTFILE
startline=3
x=0;y=0;z=0
#x=$(tail -n +$startline $INPUTFILE | awk 'BEGIN {max = -9999; min=9999} {if ($2+0 > max) max=$2; if ($2+0 < min) min=$2 } END {printf "%.10f" , 0.5*(max+min) }')
#y=$(tail -n +$startline $INPUTFILE | awk 'BEGIN {max = -9999; min=9999} {if ($3+0 > max) max=$3; if ($3+0 < min) min=$3 } END {printf "%.10f" , 0.5*(max+min) }')
#z=$(tail -n +$startline $INPUTFILE | awk 'BEGIN {max = -9999; min=9999} {if ($4+0 > max) max=$4; if ($4+0 < min) min=$4 } END {printf "%.10f" , 0.5*(max+min) }')
echo "%Coordinates" >> $OUTPUTFILE
tail -n +$startline $INPUTFILE | awk '{printf "'\''%s'\'' | angstrom* %15.10f | angstrom* %15.10f | angstrom* %15.10f \n",$1,$2-'$x',$3-'$y',$4-'$z'}'  >> $OUTPUTFILE
#tail -n +$startline $INPUTFILE | awk '{printf "'\''%s'\'' |  %15.10f |  %15.10f |  %15.10f \n",$1,$2-'$x',$3-'$y',$4-'$z'}'  >> $OUTPUTFILE
echo "%" >> $OUTPUTFILE
else
echo "Can't find "$INPUTFILE "breakdown"
exit
fi
}

xyz2oct
mpirun octopus | tee result

cp work-geom.xyz work-geom.bak.xyz
cp result result.bak
echo geom.log > geom.log 
#cg算法找不到极小值会报错不输出min.xyz
#fire虽然会输出min.xyz但是可能达到了最大步数,这时换成cg算法看可不可以
while  (($(grep 'min.xyz' result | wc -l ) < 1)) || [ $(grep 'GOMethod' inp | awk '{ print $3 }') == fire ] ;
do
echo CNQ $(date) >> result
echo "---------------" >> geom.log
for i in $(ls geom); do echo $i >> geom.log; cat geom/$i >> geom.log ; done

#收集上一次计算用的算法
GOMethod=steep
if [ $(grep 'GOMethod' inp | awk '{ print $3 }') == fire ];     then GOMethod=steep ; fi
if [ $(grep 'GOMethod' inp | awk '{ print $3 }') == steep ];    then GOMethod=cg_fr ; fi
if [ $(grep 'GOMethod' inp | awk '{ print $3 }') == cg_fr ];  then GOMethod=cg_bfgs ; fi
if [ $(grep 'GOMethod' inp | awk '{ print $3 }') == cg_bfgs ]; then GOMethod=cg_bfgs2 ; fi
if [ $(grep 'GOMethod' inp | awk '{ print $3 }') == cg_bfgs2 ];    then GOMethod=cg_pr ; fi
#注释此行则全用类cg算法
#if [ $(grep 'GOMethod' inp | awk '{ print $3 }') == cg_fr ];    then GOMethod=fire ; fi
if [ $(grep 'GOMethod' inp | awk '{ print $3 }') == cg_fr ];    then GOMethod=steep ; fi

xyz2oct
#>[可选]准备更换演化算法
echo "Change GOMethod = $GOMethod"; sed -i "s/GOMethod/GOMethod = $GOMethod #/g" inp
mpirun  octopus > result #报错和其他日志信息在slurm的输出
#可以通过修改参数,然后touch stop换算法、参数续算
cat work-geom.xyz >> work-geom.bak.xyz
cat result >> result.bak
if [ -f stop1 ]; then exit;fi
done

echo "Relax done"

性质计算

输入

#基态计算
CalculationMode = gs
#单位
UnitsOutput = ev_angstrom
#以原点的球
BoxShape = sphere
#dr
spacing = 0.19*angstrom
#球的半径,注意直径是2倍
radius = 10*angstrom
#计算的空态,在需要td投影时可以增加空态
ExtraStates=20
#输出信息
Output =dos+ density
#输出格式
OutputFormat = xcrysden
%Coordinates
'H' | angstrom*    0.0000000000 | angstrom*    0.0000000000 | angstrom*    0.0000000000
%

Output参数

可以不断向Output上增加值控制输出

#输出信息
Output =  density
Output = Output + dos

密度

输出density.xsf

Output = density
OutputFormat = xcrysden

dos & pdos

Output = Output + dos
DOSComputePDOS=yes

会输出dos,每一个能级的dos,pdos,多k点时,输出的是每条能带的dos(共nst个)

(python37) cndaqiang@mommint:~/work/octopus/H$ ls static/
convergence   dos-0003.dat  dos-0007.dat  dos-0011.dat  dos-0015.dat  dos-0019.dat  info
density.xsf   dos-0004.dat  dos-0008.dat  dos-0012.dat  dos-0016.dat  dos-0020.dat  pdos-at001-H1s.dat
dos-0001.dat  dos-0005.dat  dos-0009.dat  dos-0013.dat  dos-0017.dat  dos-0021.dat  total-dos.dat
dos-0002.dat  dos-0006.dat  dos-0010.dat  dos-0014.dat  dos-0018.dat  forces        total-dos-efermi.dat

画dos

波函数

wfs_fourier: (Experimental) Outputs wavefunctions in Fourier space. This is only implemented for the ETSF file format output. The file will be called wfs-pw-etsf.nc.
wfs: Outputs wavefunctions. Which wavefunctions are to be printed is specified by the variable OutputWfsNumber – see below. The output file is called wf-, or lr_wf- in linear response.
wfs_sqmod: Outputs modulus squared of the wavefunctions. The output file is called sqm-wf-. For linear response, the filename is sqm_lr_wf-.

示例

Output = Output + wfs #带符号
Output = Output + wfs_sqmod
OutputWfsNumber="1-5" #默认全部
OutputFormat = xcrysden
OutputInterval=5 #输出间隔，含时不要都输出，非常占内存

输出

(python37) cndaqiang@mommint:~/work/octopus/H2O$ ls static/
convergence      sqm-wf-k001-st0001.xsf  sqm-wf-k002-st0001.xsf  wf-k001-st0001.xsf  wf-k002-st0001.xsf
density-sp1.xsf  sqm-wf-k001-st0002.xsf  sqm-wf-k002-st0002.xsf  wf-k001-st0002.xsf  wf-k002-st0002.xsf
density-sp2.xsf  sqm-wf-k001-st0003.xsf  sqm-wf-k002-st0003.xsf  wf-k001-st0003.xsf  wf-k002-st0003.xsf
forces           sqm-wf-k001-st0004.xsf  sqm-wf-k002-st0004.xsf  wf-k001-st0004.xsf  wf-k002-st0004.xsf
info             sqm-wf-k001-st0005.xsf  sqm-wf-k002-st0005.xsf  wf-k001-st0005.xsf  wf-k002-st0005.xsf

如果是没有自旋分辨的单k点

(python37) cndaqiang@mommint:~/work/octopus/H2O$ !ls
ls static/
convergence  info               sqm-wf-st0003.xsf  wf-st0001.xsf  wf-st0004.xsf
density.xsf  sqm-wf-st0001.xsf  sqm-wf-st0004.xsf  wf-st0002.xsf  wf-st0005.xsf
forces       sqm-wf-st0002.xsf  sqm-wf-st0005.xsf  wf-st0003.xsf

在含时的情况，会输出各个波函数随时间的变化

(python37) cndaqiang@mommint:~/work/octopus/H2O$ ls output_iter/td.0000025
density.xsf         wf-st0002.real.xsf  wf-st0004.real.xsf  wf-st0006.real.xsf  wf-st0008.real.xsf
wf-st0001.imag.xsf  wf-st0003.imag.xsf  wf-st0005.imag.xsf  wf-st0007.imag.xsf  wf-st0009.imag.xsf
wf-st0001.real.xsf  wf-st0003.real.xsf  wf-st0005.real.xsf  wf-st0007.real.xsf  wf-st0009.real.xsf
wf-st0002.imag.xsf  wf-st0004.imag.xsf  wf-st0006.imag.xsf  wf-st0008.imag.xsf

输出文件

restart/gs/states总轨道数nst,k点数nik(spin=2时k点加倍),波函数分量dim(spin=4时,二分量波函数dim=2) 对于3个电子的体系,空轨道是20,则

SpinComponents=spin_polarized:每个自旋nst=占据(1.5)+未占据(20)=22,两种自选态在两套k点上计算,nik=2

(python37) cndaqiang@mommint:~/work/octopus/H$ cat restart/gs/states
nst=                        22
dim=                         1
nik=                         2

SpinComponents=spinors:每种自旋态只占据一个电子,nst=占据(3)+未占据(20)=23,一套k点上计算,nik=1,波函数二分量dim=2

(python37) cndaqiang@mommint:~/work/octopus/H$ cat restart/gs/states
nst=                        23
dim=                         2
nik=                         1

static

static/convergence scf的收敛信息
static/info 网格,本征值,能量,Dipole,收敛性,受力等诸多信息

TDOutput

TDOutput默认情况是输出coordinates energy laser multipoles temperature文件在td.general目录的,如果设置了具体输出TDOutput =td_occup内容,那么其他文件不会输出了,我们要手动添加输出

TDOutput=multipoles+laser+energy #+td_occup

multipoles: Outputs the (electric) multipole moments of the density to the file td.general/multipoles. This is required to, e.g., calculate optical absorption spectra of finite systems. The maximum value of l can be set with the variable TDMultipoleLmax.

angular: Outputs the orbital angular momentum of the system to td.general/angular, which can be used to calculate circular dichroism.

gauge_field: If set, outputs the vector gauge field corresponding to a spatially uniform (but time-dependent) external electrical potential. This is only useful in a time-dependent periodic run. On by default if GaugeVectorField is set.

temperature: If set, the ionic temperature at each step is printed. On by default if MoveIons = yes.

ftchd: Write Fourier transform of the electron density to the file ftchd.X, where X depends on the kick (e.g. with sin-shaped perturbation X=sin). This is needed for calculating the dynamic structure factor. In the case that the kick mode is qbessel, the written quantity is integral over density, multiplied by spherical Bessel function times real spherical harmonic. On by default if TDMomentumTransfer is set.

dipole_velocity: When set, outputs the dipole velocity, calculated from the Ehrenfest theorem, in the file td.general/velocity. This file can then be processed by the utility oct-harmonic-spectrum in order to obtain the harmonic spectrum.

eigenvalues: Write the KS eigenvalues.

ionization_channels: Write the multiple-ionization channels using the KS orbital densities as proposed in C. Ullrich, Journal of Molecular Structure: THEOCHEM 501, 315 (2000).

total_current: Output the total current (average of the current density over the cell).

partial_charges: Bader and Hirshfeld partial charges. The output file is called ‘td.general/partial_charges’.

td_kpoint_occup: Project propagated Kohn-Sham states to the states at t=0 given in the directory restart_proj (see %RestartOptions). This is an alternative to the option td_occup, with a formating more suitable for k-points and works only in k- and/or state parallelization

td_floquet: Compute non-interacting Floquet bandstructure according to further options: TDFloquetFrequency, TDFloquetSample, TDFloquetDimension. This is done only once per td-run at t=0. works only in k- and/or state parallelization

spin: (Experimental) Outputs the expectation value of the spin, which can be used to calculate magnetic circular dichroism.

n_excited_el: Output the number of excited electrons, based on the projections of the time evolved wave-functions on the ground-state wave-functions. The output interval of this quantity is controled by the variable TDOutputComputeInterval

coordinates_sep: Writes geometries in a separate file.

velocities_sep: Writes velocities in a separate file.

forces_sep: Writes forces in a separate file.

total_heat_current: Output the total heat current (average of the heat current density over the cell).

total_magnetization: Writes the total magnetization, where the total magnetization is calculated at the momentum defined by TDMomentumTransfer. This is used to extract the magnon frequency in case of a magnon kick.

populations: (?好像是多体波函数激发态的投影)(Experimental) Outputs the projection of the time-dependent Kohn-Sham Slater determinant onto the ground state (or approximations to the excited states) to the file td.general/populations. Note that the calculation of populations is expensive in memory and computer time, so it should only be used if it is really needed. See TDExcitedStatesToProject.

geometry: If set (and if the atoms are allowed to move), outputs the coordinates, velocities, and forces of the atoms to the the file td.general/coordinates. On by default if MoveIons = yes.

dipole_acceleration: When set, outputs the acceleration of the electronic dipole, calculated from the Ehrenfest theorem, in the file td.general/acceleration. This file can then be processed by the utility oct-harmonic-spectrum in order to obtain the harmonic spectrum.

laser: If set, outputs the laser field to the file td.general/laser. On by default if TDExternalFields is set.

energy: If set, octopus outputs the different components of the energy to the file td.general/energy. Will be zero except for every TDEnergyUpdateIter iterations.

td_occup: (向0时刻的ks投影)(Experimental) If set, outputs the projections of the time-dependent Kohn-Sham wavefunctions onto the static (zero-time) wavefunctions to the file td.general/projections.XXX. Only use this option if you really need it, as it might be computationally expensive. See TDProjStateStart. The output interval of this quantity is controled by the variable TDOutputComputeInterval In case of states parallelization, all the ground-state states are stored by each task.还需要设置ExperimentalFeatures=yes,间隔由TDOutputComputeInterval=50决定

local_mag_moments: If set, outputs the local magnetic moments, integrated in sphere centered around each atom. The radius of the sphere can be set with LocalMagneticMomentsSphereRadius.

搭建说明网页

安装apache&php,复制doc目录到网站目录后,主机hosts为cndaqiang时浏览器打开http://cndaqiang/octopus/doc/html/vars.php?page
注意chrome和edge必须有?page,不然打开的是空白网页,safari不加也可以打开,原因暂时不明

或者直接打开相应的html文件,就是丑了点

示例

CalculationMode = gs #计算方式 gs:基态(默认) td:时间依赖 unocc:非自洽计算  
UnitsOutput =  eV_Angstrom  #默认atomic原子单位制(能量是Hartree)

Nitrogen_mass=14.0 #可以在inp文件中定义变量，支持+-*/sqrt()等多种数学运算
	#更多数学运算https://octopus-code.org/wiki/Manual:Input_file#Mathematical_expressions


%Species  #元素描述block
 'N' | species_pseudo | set | standard | lmax | 1 | lloc | 0 | mass | Nitrogen_mass
#原素名 | 性质1 | 值 | 性质2 | 值
#赝势 species_pseudo | set默认就是用系统的 | standard
#     species_pseudo | file | `赝势文件名，支持 .psf(SIESTA),.upf,.fhi,.cpi,.xml,.hgh`  
# 更多参考https://octopus-code.org/doc/8.4/html/vars.php?section=System&name=Species
#lmax: The maximum angular-momentum channel that will be used for the pseudopotential.
#lloc: The angular-momentum channel of the pseudopotential to be considered local.
%


XYZCoordinates = 'N.xyz'
ExtraStates = 1

%Occupations
 2|1|1|1
%

BoxShape = sphere

Radius = 5.0*angstrom
#Spacing = 0.18*angstrom  #时空间网格大小

#参数也可以通过环境变量设置
#设置linux变量
#```
#export OCT_PARSE_ENV=1 #开启从环境读入功能
#export OCT_Spacing=0.18*angstrom  #以OCT_开头设置参数
#```

Output = wfs + density + elf + potential
#默认只产生info，显示能量，特征值，等，其他信息，设置输出
#wfs wavefunction ; elf : the electron localization function  ;potential : the Kohn-Sham potential
#dos : 态密度

OutputFormat = cube + xcrysden + dx + axis_x + plane_z
#输出文件格式，用于给不同的软件(vislt,gunplot , xcrysden,UCSF(density.dx))打开
#https://octopus-code.org/wiki/Tutorial:Visualization

#octopus目录中有很多小程序可以用

#支持 include 文件

周期性示例

CalculationMode = gs

PeriodicDimensions = 3
Spacing = 0.5 

%LatticeVectors
0.0 | 0.5 | 0.5
0.5 | 0.0 | 0.5
0.5 | 0.5 | 0.0
%

a=10.18
%LatticeParameters
 a | a | a
%

%ReducedCoordinates
 "Si" | 0.0 | 0.0 | 0.0 
 "Si" | 1/4 | 1/4 | 1/4 
%

nk=2

%KPointsGrid
  nk |  nk |  nk  #第一行：倒空间每个方向k点数量，na | nb | nc 总数na*nb*nc Monk-Pack网格
 0.5 | 0.5 | 0.5  #之后的为K点的移位
 0.5 | 0.0 | 0.0
 0.0 | 0.5 | 0.0
 0.0 | 0.0 | 0.5
%

KPointsUseSymmetries = yes

ExtraStates = 1
Output = dos

td 计算示例

CalculationMode = td #gs
UnitsOutput = eV_Angstrom

Radius = 3.5*angstrom
Spacing = 0.18*angstrom

CH = 1.097*angstrom
%Coordinates
  "C" |           0 |          0 |           0
  "H" |  CH/sqrt(3) | CH/sqrt(3) |  CH/sqrt(3)
  "H" | -CH/sqrt(3) |-CH/sqrt(3) |  CH/sqrt(3)
  "H" |  CH/sqrt(3) |-CH/sqrt(3) | -CH/sqrt(3)
  "H" | -CH/sqrt(3) | CH/sqrt(3) | -CH/sqrt(3)
%

Tf = 1/eV
dt = 0.002/eV

TDPropagator = aetrs
TDMaxSteps = Tf/dt  #总步数
TDTimeStep = dt     #步长

FromScratch = yes
 
amplitude = 1*eV/angstrom
omega = 18.0*eV
tau0 = 0.5/eV
t0 = tau0

%TDExternalFields
  electric_field | 1 | 0 | 0 | omega | "envelope_cos"
#  电场/磁场/矢量势/标量势/ f(x,y,z) | nx | ny | nz 场的极化方向(nx,ny,nz) | omega | g(t) 字符型，定义在%TDFunctions中 |[ 相位phase(t) ]
# 外场 = f(x,y,z)*cos(w*t+phase(t))*g(t)
# f(x,y,z) : electric field, magnetic field, vector_potential
# g(t)与phase(t) 都是字符串，字符串内容随意，定义式在%TDFunctions中
#phase(t) 不设置时默认为0
%
 
%TDFunctions
  "envelope_cos" | tdf_cosinoidal | amplitude | tau0 | t0
# %TDExternalFields中g(t),phase(t) | 函数名:tdf_cw,tdf_gaussian,tdf_cosinoidal,...来自文件，自定义 | 相关参数
# https://octopus-code.org/doc/8.4/html/vars.php?section=Time-Dependent&name=TDFunctions

%

TDOutput = laser + multipoles
#TD输出
#laser: If set, outputs the laser field to the file td.general/laser. On by default if TDExternalFields is set.
#multipoles: Outputs the (electric) multipole moments of the density to the file td.general/multipoles. This is required to, e.g., calculate optical absorption spectra of finite systems. The maximum value of l can be set with the variable TDMultipoleLmax.

td计算输文件td.general里面包含TDOutput设定的文件

报错

使用中文的tab键

Input: [PseudopotentialSet = standard]
Reading Coordinates from Coordinates block

**************************** FATAL ERROR *****************************
*** Fatal Error (description follows)
*--------------------------------------------------------------------
* From node =    0
*--------------------------------------------------------------------
* Error in block Coordinates line # 1
**********************************************************************

核数不合适

并行原因: 改成2，4，5，10个核都没有问题，6个核就报错
核数不能超过总轨道数(occ+ExtraStates)

      Block       3 contains       1 states:       4 -       4
Info: Ground-state restart information will be read from 'restart/gs'.

           Info: Reading states: gs. 2021/03/23 at 22:18:22

ETA: .
Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

前面的并行信息也许可以参考

************************** Parallelization ***************************
Info: Octopus will run in *parallel*

      Number of processes           :      10
      Number of threads per process :       1

Info: Number of nodes in ParStates   group:     2 (       4)
Info: Number of nodes in ParDomains  group:     5 (  131995)
Info: Octopus will waste at least  0.00% of computer time.
**********************************************************************

************************** Parallelization ***************************
Info: Octopus will run in *parallel*

      Number of processes           :       4
      Number of threads per process :       1

Info: Number of nodes in ParStates   group:     4 (       4)
Info: Octopus will waste at least  0.00% of computer time.
**********************************************************************

************************** Parallelization ***************************
Info: Octopus will run in *parallel*

      Number of processes           :       5
      Number of threads per process :       1

Info: Number of nodes in ParDomains  group:     5 (  131995)
Info: Octopus will waste at least  0.00% of computer time.
**********************************************************************

************************** Parallelization ***************************
Info: Octopus will run in *parallel*

      Number of processes           :       6
      Number of threads per process :       1

Info: Number of nodes in ParStates   group:     3 (       4)
Info: Number of nodes in ParDomains  group:     2 (  131995)

** Warning:
**   Info: Octopus will waste at least 33.33% of computer time.
**   Usually a number of processors which is a multiple of small primes is best.

radius太小

算PDOS时出现,需要比自洽更大的radius

**************************** FATAL ERROR *****************************
*** Fatal Error (description follows)
*--------------------------------------------------------------------
* From node =    0
*--------------------------------------------------------------------
* The radius of an orbital set is bigger than the radius of the simulatio box.
* Increase the value of Radius or decrease the value of OrbitalsThreshold_LDAU.
* The value of the radius is  5.90000 Bohr.
**********************************************************************

*--------------------------------------------------------------------
* An orbital set has points outside of the simulatio box.
* Increase the value of Radius or decrease the value of OrbitalsThreshold_LDAU.
* The value of the radius is 11.00000 Bohr.

这里的The value of the radius is 5.90000 Bohr.的radius是轨道的半径,不是模拟体系的半径, 增加模拟尺寸,如模拟半径radius = 3.2*angstrom即可,下面可以证明设置AO开头的Atomic Orbitals参数无效.

报错程序src/basis_set/atomic_orbital_inc.F90

! ---------------------------------------------------------
!> This routine returns the atomic orbital basis -- provided
!! by the pseudopotential structure in geo.
! ---------------------------------------------------------
subroutine X(get_atomic_orbital) (geo, mesh, sm, iatom, ii, ll, jj, os, orbind, radius, d_dim)
!...
  if(sm%np == -1) then

    if(mesh%sb%box_shape == MINIMUM .and. radius > mesh%sb%rsize) then
      message(1) = "The radius of an orbital set is bigger than the radius of the simulatio box."
      message(2) = "Increase the value of Radius or decrease the value of OrbitalsThreshold_LDAU."
      write(message(3),'(a,f8.5,a,i5,a)') 'The value of the radius is ', radius, ' Bohr.'
      call messages_fatal(3)
    end if

    if(mesh%sb%box_shape == SPHERE .or. mesh%sb%box_shape == CYLINDER) then
      if(radius > mesh%sb%rsize) then
       message(1) = "The radius of an orbital set is bigger than the radius of the simulatio box."
       message(2) = "Increase the value of Radius or decrease the value of OrbitalsThreshold_LDAU."
       write(message(3),'(a,f8.5,a,i5,a)') 'The value of the radius is ', radius, ' Bohr.'
       call messages_fatal(3)
      end if
      if(mesh%sb%box_shape == CYLINDER .and. radius > mesh%sb%xsize) then
       message(1) = "The radius of an orbital set is bigger than the length of the cylinder box."
       message(2) = "Increase the value of XLength or decrease the value of OrbitalsThreshold_LDAU."
       write(message(3),'(a,f8.5,a,i5,a)') 'The value of the radius is ', radius, ' Bohr.'
       call messages_fatal(3)
      end if
    end if

    if(mesh%sb%box_shape == SPHERE ) then
      if(sqrt(sum(geo%atom(iatom)%x(1:mesh%sb%dim)**2)) + radius > mesh%sb%rsize) then
       message(1) = "An orbital set has points outside of the simulatio box."
       message(2) = "Increase the value of Radius or decrease the value of OrbitalsThreshold_LDAU."
       write(message(3),'(a,f8.5,a,i5,a)') 'The value of the radius is ', radius, ' Bohr.'
       call messages_fatal(3)
      end if
    end if

    if(mesh%sb%box_shape == CYLINDER ) then
      if(sqrt(sum(geo%atom(iatom)%x(2:mesh%sb%dim)**2)) + radius > mesh%sb%rsize) then
       message(1) = "An orbital set has points outside of the simulatio box."
       message(2) = "Increase the value of Radius or decrease the value of OrbitalsThreshold_LDAU."
       write(message(3),'(a,f8.5,a,i5,a)') 'The value of the radius is ', radius, ' Bohr.'
       call messages_fatal(3)
      end if
      if(abs(geo%atom(iatom)%x(1)) + radius > mesh%sb%xsize) then
       message(1) = "An orbital set has points outside of the simulatio box."
       message(2) = "Increase the value of Xlength or decrease the value of OrbitalsThreshold_LDAU."
       write(message(3),'(a,f8.5,a,i5,a)') 'The value of the radius is ', radius, ' Bohr.'
       call messages_fatal(3)
      end if
    end if

溯源./basis_set/orbitalset_utils_inc.F90

  subroutine X(orbitalset_utils_getorbitals)(os, geo, mesh)
      ! We obtain the orbital
      call X(get_atomic_orbital)(geo, mesh, os%sphere, os%iatom, os%ii, os%ll, os%jj, &
                                              os, iorb, os%radius, os%ndim)

溯源./system/dos.F90

              os%radius = atomic_orbital_get_radius(geo, mesh, ia, iorb, 1, &
                                OPTION__AOTRUNCATION__AO_FULL, CNST(0.01))

最终决定半径的是./basis_set/atomic_orbital.F90

  ! ---------------------------------------------------------
  FLOAT function atomic_orbital_get_radius(geo, mesh, ia, iorb, ispin, truncation, threshold) result(radius)
    type(geometry_t), target, intent(in)   :: geo
    type(mesh_t),             intent(in)   :: mesh
    integer,                  intent(in)   :: ia, iorb, ispin
    integer(8),               intent(in)   :: truncation
    FLOAT,                    intent(in)   :: threshold

    type(species_t), pointer :: spec
    integer :: ii, ll, mm

    PUSH_SUB(atomic_orbital_get_radius)

    spec => geo%atom(ia)%species
    call species_iwf_ilm(spec, iorb, ispin, ii, ll, mm)

    if(truncation == OPTION__AOTRUNCATION__AO_FULL) then
      radius = species_get_iwf_radius(spec, ii, ispin, threshold)
    else
      radius = species_get_iwf_radius(spec, ii, ispin)

      if(truncation == OPTION__AOTRUNCATION__AO_BOX) then
        ! if the orbital is larger than the size of the box, we restrict it to this size,
        ! otherwise the orbital will overlap more than one time with the simulation box.
        ! This would induces phase problem if the complete mesh is used instead of the sphere
        radius = min(radius, minval(mesh%sb%lsize(1:mesh%sb%dim)-mesh%spacing(1:mesh%sb%dim)*CNST(1.01)))
      else
        !If asked, we truncate the orbital to the radius on the projector spheres
        !of the NL part of the pseudopotential.
        !This is a way to garanty no overlap between orbitals of different atoms.
        if(species_is_ps(spec)) &
          radius = min(radius,species_get_ps_radius(spec))
        end if

      end if
      ! make sure that if the spacing is too large, the orbitals fit in a few points at least
      radius = max(radius, CNST(2.0)*maxval(mesh%spacing(1:mesh%sb%dim)))

    POP_SUB(atomic_orbital_get_radius)
  end function atomic_orbital_get_radius

可以发现./system/dos.F90强行限制了AOTruncation为ao_full,按照手册

ao_full: The full size of the orbitals used. The radius is controled by variable AOThreshold.

然而我们发现dos.F90也设置了threshold=CNST(0.01),所以无法通过参数调整,唯一的途径就只有调模拟尺寸的大小

缺少赝势

并行运行时，没有显示

Input: [PseudopotentialSet = standard]
Reading Coordinates from Coordinates block
Abort(269076484) on node 0 (rank 0 in comm 0): Fatal error in PMPI_Recv_init: Invalid tag, error stack:
PMPI_Recv_init(147): MPI_Recv_init(buf=0x7c0dd80, count=1, MPI_INTEGER, src=1, tag=1620299, MPI_COMM_WORLD, request=0x7c0dd40) failed
PMPI_Recv_init(99).: Invalid tag, value is 1620299

单核运行发现是缺少赝势

************************** Calculation Mode **************************
Input: [CalculationMode = gs]
**********************************************************************

Input: [PseudopotentialSet = standard]
Reading Coordinates from Coordinates block

**************************** FATAL ERROR *****************************
*** Fatal Error (description follows)
*--------------------------------------------------------------------
* From node =    0
*--------------------------------------------------------------------
* Species Ag not found in default pseudopotential set.
* ( /home/apps/octopus10.4/share/octopus/pseudopotentials/PSF )
**********************************************************************

Abort(999) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 999) - process 0

添加赝势到/home/apps/octopus10.4/share/octopus/pseudopotentials/PSF就可以了

            Calculation started on 2021/12/16 at 10:46:32


************************** Calculation Mode **************************
Input: [CalculationMode = gs]
**********************************************************************

Input: [PseudopotentialSet = standard]
Reading Coordinates from Coordinates block

****************************** Species *******************************
  Species 'Ag'
    type             : pseudopotential
    file             : '/home/apps/octopus10.4/share/octopus/pseudopotentials/PSF/Ag.psf'
    file format      : PSF
    valence charge   : 11.0
    atomic number    :  47
    form on file     : semilocal
    orbital origin   : calculated
    lmax             : 3
    llocal           : 0
    projectors per l : 1

本文首发于我的博客@cndaqiang.
本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议，转载请注明出处！