cndaqiang Web Linux DFT

octopus 快速上手

2019-04-11
cndaqiang
RSS

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

参考

octopus-Tutorial

Gcc OPENMPI 编译octopus-8 遇到的问题和解决方案

centos6.5 Gcc Openmpi 编译octopus-7.1

简单运行

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之类的

控制参数

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
%

性质计算

输入

#基态计算
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
    

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强行限制了AOTruncationao_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 协议 ,转载请注明出处!



评论


广告

目录

广告
访客数据