Jacek Kobus

2 Dimensional Finite Difference

Hartree-Fock Program

User's Guide

ver. 2.0 (June 2010)

Contents

• Description of input data
  • Mandatory labels
  • Optional labels
  • Deprecated labels
  
  
• Examples of input data
• Program's data files
• How to run the program?
• Some hints

Description of input data

An input data file consists of separate lines each containing

-

a label

-

a label followed by a string of characters, integer(s) and/or real number(s)

-

a string of characters, integer(s) and/or real number(s)

Real numbers can be written in a fix point or scientific notation. Note that

-

labels and strings can be in upper or lower case,

-

the compulsory labels must follow the order given below; the optional ones can be inserted anywhere between the title and stop labels,

-

optional parameters are enclosed in square brackets,

-

$r$ denotes a real number, $i$ - an integer, $c$ - a string of characters,

-

an exclamation mark or a hash placed anywhere in an input line starts a comment and what follows ``!'' or ``#'' is ignored.

Mandatory labels

The following labels must be specified in the given order:

TITLE

Format:

title $c$
$c$ is any string of up to 74 characters describing the current case. This string is added as a header to a text file with extension dat that contains basic data identifying a given case, i.e. atomic numbers of nuclei, grid size and the number of electrons and orbital and exchange functions.

NUCLEI

Format:

nuclei $Z_A \;\;Z_B\;\; R \;\;[\; c\;]$
Set the nuclei charges and the bond length.

$Z_A$ :

nuclear charge of centre A (real)

$Z_B$ :

nuclear charge of centre B (real)

$R$ :

bond length (real)

$c$ :

$angstrom$ - the internuclear separation can be given in angstrom units if this string is included (the conversion factor 0.529177249 is used)

If $\vert Z_A-Z_B\vert<10^{-6}$ then the molecule is considered to be a homonuclear one (this threshold can be changed by redefining HOMOLEVL variable in blk_data).

CONFIG

Format:

config $i$

$i$ :

the total charge of a system

The following cards define molecular orbitals and their occupation. Note that the last orbital description card must contain the $end$ label.

The possible formats are:

Format:

$i\;\;\;c$

$i$ :

number of fully occupied orbitals of a given irreducible representation (irrep) of the $C_{\infty\,v}$ group; two electrons make $\sigma$  orbitals fully occupied and four electrons are needed for filling the orbitals of other symmetries

$c$ :

symbol of the $C_{\infty\,v}$ irrep to which the orbitals belong (sigma, pi, delta or phi)

Format:

$i\;\;\;c_1 \;\;c_2$

$i$ :

number of fully occupied orbitals of a given irrep of the $D_{\infty\,h}$ group

$c_1$ :

symbol of the $C_{\infty\,v}$ irrep to which the orbitals belong (sigma, pi, delta or phi)

$c_2$ :

symbol for the inversion symmetry of the $D_{\infty\,h}$ irreps (u or g)

Use this format for a homonuclear molecule unless break card is included (see below).

Format:

$i\;\;\;c_1\;\;c_2\;\;[c_3\;\;[\;c_4\;\;[\;c_5\;]\;]\;]$

$i$ :

number of orbitals of a given irrep of the $C_{\infty\,v}$ group

$c_1$ :

symbol for the $C_{\infty\,v}$ irreps to which the orbitals belong (sigma, pi, delta, phi)

$c_2$ -$c_5$ :

$+,\;-$ or . (a dot); $+/-$ denotes spin up/down electron and . denotes an unoccupied spin orbital

Format:

$i\;\;\;c_1\;\;c_2\;\;c_3\;\;[\;c_4\;\;[\;c_5\;\;[\;c_6\;]\;]\;]$

$i$ :

number of orbitals of a given irrep of the $D_{\infty\,h}$ group

$c_1$ :

symbol for the $C_{\infty\,v}$ irrep to which the orbitals belong (sigma, pi, delta, phi)

$c_2$ :

symbol for the inversion symmetry of the $D_{\infty\,h}$ irrep (u or g)

$c_3$ -$c_6$ :

$+,\;-$ or . (a dot); $+/-$ denotes spin up/down electron and . denotes an unoccupied spin orbital

GRID

Two possible formats are (the second one is deprecated and is retained for backward compatibility):

Format:

grid $n_{\nu}$ $R_{\infty}$
An integer and a real define a single two-dimensional grid.

$n_{\nu}$ :

the number of grid points in $\nu$ variable

$R_{\infty}$ :

the practical infinity

$n_{\mu}$ is calculated so as to make the step size in $\mu$ variable equal to the stepsize in $\nu$ variable. $n_{\nu}$ and $n_{\mu}$ have to meet special conditions. If the conditions are not fulfilled the nearest (but smaller) appropriate values are used.

Format:

grid $n_{\nu}$ $n_{\mu}$ $R_{\infty}$
Two integers and one real define a single two-dimensional grid.

$n_{\nu}$ :

the number of grid points in $\nu$ variable

$n_{\mu}$ :

the number of grid points in $\mu$ variable

$R_{\infty}$ :

the practical infinity

$n_{\nu}$ and $n_{\mu}$ have to meet special conditions. If the conditions are not fulfilled the nearest (but smaller) appropriate values are used.

ORBPOT

Format:

orbpot $c$
where $c$ a character string determining the initial source of orbitals and potentials. Its allowed values are:

• $hydrogen$ - molecular orbitals are formed as a linear combination of hydrogenic functions on centres $A$ and $B$ as defined via the lcao label. In the case of HF or HFS calculations Coulomb (exchange) potentials are approximated as a linear combination of Thomas-Fermi ($1/r$ ) potentials at the two centres. If the OED method is chosen the potential function is approximated as a linear combination of $Z_A/r_1$ and $Z_B/r_2$ terms and the exchange potentials are set to zero.

• $gauss$ - GAUSSIAN94/GAUSSIAN98 output is used to retrieve exponents and expansion coefficients of (uncontracted) molecular orbitals (it is assumed that the output is contained in gaussian.out and gaussian.pun files) and Coulomb and exchange potentials are initialized as in the hydrogen case; see routine PREPGAUSSIAN for more details.

• $gauss-c$ - GAUSSIAN94/GAUSSIAN98 customized output is used to retrieve exponents and expansion coefficients of molecular orbitals (it is assumed that the output is contained in gaussianc.out file) and Coulomb and exchange potentials are initialized as in the hydrogen case; see routine prepGaussCust for details.

• $old$ - initial orbitals, Coulomb and exchange potentials are retrieved from disk files (2dhf_input.orb, 2dhf_input.coul and 2dhf_input.exch, respectively) created in a previous run. Data defining the case are retrieved from a 2dhf_input.dat textfile.

• $qrhf$ - radial Hartree-Fock orbitals for the centre A and B obtained from the qrhf program¹ are retrieved from disk files 1dhf_inputA.orb and 1dhf_inputB.orb, respectively, and Coulomb and exchange potentials are initialized as in the hydrogen case (see routine initHF for details).

• $noexch$ - orbitals and Coulomb potentials are retrieved from disk files and exchange potentials are initialized as in $hydrogen$ case; this is useful when going from HFS to HF calculations.

• $nodat$ - initial orbitals and potentials are retrieved from disk files but the content of a 2dhf_input.dat file is retrieved from a 2dhf_input.orb (binary) file. Use this value when reading binary data generated by the ealier than 2.0 versions of the program.

STOP

Format:

stop
This label indicates the end of input data.

Optional labels

The following additional labels can be specified in any order:

BREAK

Format:

break
When this label is present homonuclear molecules are calculated in $C_{\infty\,v}$ symmetry and the $D_{\infty\,h}$ symmetry labels (u or g) are superfluous.

CONV

Format:

conv $[\;i_1\;[\;i_2\;[\;i_3\;]\;]\;]$
Sometimes the convergence threshold (for energy and/or normalization) is set too low and cannot be satisfied on a given grid and as a result the SCF/SOR process continues in vain. The iterations are stopped if orbital energies or orbital norms display no improvment over a given number of $i_2$ and $i_3$ most recent iterations, respectively (the default values are set to 20). The monitoring begins after $i_1$  initial iterations (default 600).

DEBUG

Format:

debug $[\; i_1 \; [\;i_2 \ldots \;] \;i_{40}\;]$

Up to 40 different debug flags can be set at a time². If the integer $i_k$ is encountered the debug flag $i_k$ is set, i.e. idbg $(i_k)=1\;\; (1 \leq i_k<999)$ . These are used to generate additional printouts by adding where needed the lines of the form

 
       if (idbg(ik).eq.1) then
           print *, ``debuging something ...''
           ...
      endif

These are used to generate debug information.

DFT

Format:

dft $[\;c_1\;]\;\; [\;c_2\;]$

$c_1:$

specifies the type of DFT exchange potential to be used in Fock equations

• $c_1=lda$ - the local density approximation with the potential

  $$V_X(\alpha)=-\frac{3}{2} \alpha \left(\frac{3}{\pi}\right)^{1/3} 2^{-2/3} \sum_{\sigma} \rho_{\sigma}^{1/3}$$

  
  
  where $\alpha$ is by default set to 2/3 (the Slater exchange potential). To change this value use the xalpha label.

• $c_1=b88$ - the Becke exchange potential

$c_2:$

selects the type of correlation potential to be used in Fock equations

• $c_2=lyp$ - the correlation potential of Lee, Yang and Parr

• $c_2=vwn$ - the correlation potential of VWN

When the bare label is present and the method selected is HF then the exchange contributions (LDA, B88, PW86 and PW91) and the correlation contributions (LYP and VWN) to the total energy are calculated upon completion of the SCF iterations.

EXCHIO

Format:

exchio $c_1$ $c_2$
This is a replacement for $i_2$ parameter of INITIAL card (especially when ORBPOT card is used to define the source of orbitals and potentials). These parameters specify how exchange potentials are to be read/written and manipulated (stored in core memory). The program always keeps all orbitals and Coulomb potentials in the memory. If computer resources are adequate all exchange potentials can also be all kept there. However, during the relaxation of a particular orbital only a fraction of exchange potentials is needed. Thus all exchange potentials can be kept on disk as separate files (named fort.31, fort.32, $\ldots$ during a run) and only relevant ones are being retrieved when necessary.

The possible values of $c_1$ and $c_2$ are in-one, in-many, out-one and out-many. Possible combinations are

• exchio in-many out-many - read exchange potentials as separate files and write them back as separate files (equivalent to setting $i_2=0$ on INITIAL card)

• exchio in-one out-many - read all exchange potentials in one file but write them out as separate files (same as $i_2=1$ )
• exchio in-many out-one - read all exchange potentials separately but write them out as a single file (same as $i_2=2$ )
• exchio in-one out-one - read and write exchange potentials in the form of a single file (same as $i_2=3$ ); this is the default

FEFIELD

Format:

fefield $r$

$r$ :

a strength of an external static electric field directed along the internuclear axis (in atomic units)

FERMI

Format:

fermi $r_A$ $r_B$
When this label is present, the Fermi nuclear charge distribution is used. Optional parameters $r_A$ and $r_B$ define the atomic masses (in amu) of nuclei A and B. If omitted the corresponding values are taken from the table of atomic masses compiled by Wapstra and Audi (see blk_data).

FIXORB

Format:

fixorb $[\; i_1 \; [\;i_2 \ldots \;] \;i_{40}\;]$
This label is used to specify orbitals to be kept frozen during SCF/SOR process. $i_1$ , $i_2$ , $\ldots$ are the numbers of these orbitals as they appear on the program's listing, i.e. their order is reversed to that used when defining the electronic configuration (see the config card). Up to 40 different orbitals can be set at a time. Use the bare label to keep all orbitals frozen but allow for the relaxation of potentials.

FIXCOUL

Format:

fixcoul
If this label is present then all Coulomb potentials are kept frozen during the SCF/SOR process.

FIXEXCH

Format:

fixexch
If this label is present then all exchange potentials are kept frozen during the SCF/SOR process.

GAUSS

Format:

gauss $r_A$ $r_B$
When this label is present, the Gauss nuclear charge distribution is used. Optional parameters $r_A$ and $r_B$ define the atomic masses (in amu) of nuclei A and B. If omitted the corresponding values are taken from the table of atomic masses compiled by Wapstra and Audi (see blk_data).

HOMO

Format:

homo
This label is used to impose explicitly $D_{\infty\,h}$ symmetry upon orbitals of homonuclear molecules in order to improve SCF/SOR convergence.

INOUT

Format:

inout $c_1$ $c_2$
The x2dhf program can be compiled to support calculation using three different combinations of integer/real data types: i32 (4-byte integers, 8-byte reals), i64 (8-byte integers, 8-byte reals) and r128 (8-byte integers, 16-byte reals); see src/Makefile for details. Strings $c_1$ and $c_2$ determine the combination appropriate for the format of input and output data, respectively, and each string can be i32, i64 or r128.

In order to facilitate exchange of binary data generated on machines of different architectures or using different compilers additional formats are available, namely i32f, i64f or r128f which allow to export/import data in formatted instead of unformatted form.

INTERP

Format:

interp
Use this label to change the grid between separate runs of the program. The restriction is that only a number of grid points in one of the variables or $R_{\infty}$ can be changed at a time.

LCAO

Format:

lcao $[\;i\;]$
If the source of orbitals is declared as hydrogen then this card must be present. In such a case the initialization of each of the orbitals has to be defined in terms of the linear combination of atom centred hydrogen-like functions. For each orbital include a card of the following format (make sure that the order of orbitals should match the order specified under the config label):

Format:

$c_A\;\;n_A \;\;l_A \;\;\zeta_A \;\;\;\;c_B\;\;n_B\;\;l_B\;\; \zeta_B\;\;\;\;i_1\;\;\;\;[\;i_2\;]$
where

$c_A$ - relative mixing coefficient for a hydrogenic orbital on the $Z_A$ centre (real),

$n_A$ - its principle quantum number (integer)

$l_A$ - its orbital quantum number (integer)

$\zeta_A$ - the effective nuclear charge if $i=1$ (default) or a screening parameter if $i=2$ (real)

$c_B$ - relative mixing coefficient for a hydrogenic orbital on the $Z_B$ centre (real),

$n_B$ - its principle quantum number (integer)

$l_B$ - its orbital quantum number (integer)

$\zeta_B$ - the effective nuclear charge if $i=1$ (default) or a screening parameter if $i=2$ (real)

$i_1$ - set to 1 to freeze the orbital during scf; otherwise set to 0 (integer)

$i_2$ - a number of successive over-relaxations for a given orbital (integer); if omitted is set to 10

The mixing coefficients are normalized so that $\vert c_A\vert+\vert c_B\vert=1$

MCSOR

Format:

mcsor $[\;i_1\;[\;i_2\;]\;]$
Selects the MCSOR method for solving the Poisson equations for orbitals and potentials (default) and changes the value of the MCSOR relaxation sweeps during a single SCF cycle for orbitals ($i_1$ ) and potentials ($i_2$ ); by default $i_1=i_2=10$ .

METHOD

Format:

method $c$
Select the type of calculation.

$c$ :

HF - the Hartree-Fock method

$c$ :

DFT - the Hartree-Fock method with the $X\alpha$ exchange potential ( $\alpha=2/3$ ); see the dft label to choose another exchange or correlation potential

$c$ :

HFS - the Hartree-Fock-Slater method (Hartree-Fock with the $X\alpha$ exchange potential) with an optimum value of the $\alpha$ parameter (see blk-data.inc for details)

$c$ :

OED - One Electron Diatomic ground and excited states can be calculated for the Coulomb potential in the prolate spheroidal coordinates (default). It is also possible to specify the Coulomb and Krammers-Henneberger potentials in cylindrical coordinates (see the poth3 and potkh labels, respectively). When more than one orbital is specified calculations are carried out as if in case of a multielectron system.³

$c$ :

SCMC - the Hartree-Fock method with $X\alpha$ exchange where the $\alpha$ parameter is calculated according to the Self-Consistent Multiplicative Constant method⁴

MULTIPOL

Format:

multipol $r$ [ $i$ ]

$r$ :

if $r>0$ multipole moment expansion coefficients are recalculated when the maximum error in orbital energy is changed by $r$ (the default value is 1.15; see setDefaults). When the maximum error is less than $r$ the coefficients are recalculated at least every $i_2$ SCF iterations (see SCF label). To suppress recalculation of the coefficients set $r$ to a negative real number. This is useful when generating potentials from a set of fixed orbitals, e.g from GAUSSIAN orbitals.

$i$ :

number of terms in the multipole expansion used to calculate boundary value for potentials ( $2 \le i\le 8$ and the default is 4)

OMEGA

Format:

omega $\omega_{orb}$ [ $\omega_{pot}$ ]
One or two real numbers setting over-relaxation parameters for relaxation of orbitals and potentials. The negative value of a given parameter indicates that its value should be set to a near optimum value obtained from a semiempirical formula (see initCBlocks and setomega for details). If this card is omitted $\omega_{pot}$ is estimated from the semiempirical formula and $\omega_{orb}=0.98 \, \omega_{pot}$ (see setDefaults to change the default value of the scaling factor).

For backward compatibility the following format is also supported:

Format:

omega
$\omega_{orb}$
$\omega_{pot}$

OMEGAOPT

Format:

omegaopt [ i [ $r_{orb}$ [ $r_{pot}$ ] ] ]
Optional integer parameter can be set to 1 (default) or 2. In the former case rather conservative (safe) omega values are chosen (this is equivalent to using the omega card with the negative values of the omega parameters). In the latter case somewhat better omega values are chosen but faster convergence is to be expected when good initial estimates of orbitals and potentials are available or when calculations with fixed orbitals or potentials are performed. The near-optimal omega values obtained from a semi-empirical formula are scaled down to produce final values used by the program. The default values of the scaling factors for the orbital and potential overrelaxation parameters (0.986 and 0.997, respectively) can be changed by setting $r_{orb}$ and $r_{pot}$ to their desired values.

ORDER

Format:

order $[\;i\;]$
An integer defining the ordering of mesh points: 1 - natural column-wise ordering, 2 - 'middle' type of sweep (default), 3 - natural row-wise, 4 - reversed natural column-wise (see mesh routine for details)

POTGSZ

Format:

potgsz
When the OED method is chosen then this label selects a model potential due to Green, Sellin and Zachor. For a given atom this potential produces HF-like orbitals but it was found useful in finding decent starting orbitals for any molecular system.

POTGSZG

Format:

potgszg
When the OED method is chosen then this label selects a model potential due to Green, Sellin and Zachor and the Gauss nuclear charge distribution. For a given atom this potential produces HF-like orbitals but it was found useful in finding decent starting orbitals for any molecular system.

POTH3

Format:

poth3 $m\;\;a\;\;V_0$
When the OED method is chosen then this label selects a two-dimensional model potential of the form $-V_0/\sqrt(a^2+r^2)$ (see kh.c for details). The following parameters can be set

$m:$

magnatic quantum number of a state (integer)

$a:$

width of the model potential (real)

$V_0:$

depth of the model potential (real)

In order to get the hydrogen Coulomb potential set $a=0.0$ and $V_0=1.0$ . Set $a>0$ and $V>0$ to choose its smoothed variant.

POTKH

Format:

potkh $m\;\;\varepsilon\;\;\omega\;\;\;[\;a\;[\;V_0\;[\;N\;]\;]\;]$
When the OED method is chosen then this label selects the Krammers-Henneberger potential (see routine kh.c for details). The following parameters can be set

$m:$

magnatic quantum number of a state being calculated

$\varepsilon:$

laser intensity

$\omega:$

laser cycle frequency

$a:$

original (before averaging over one laser cycle) soft-core potential width (a positive real number, by default $a=1.0$ )

$V_0:$

original soft-core potential depth (by default $V_0=1.0$ )

$N:$

number of intervals in the Simpson quadrature (an integer, $N=1000$ by default)

PRINT

Format:

print $[\; i_1 \; [\;i_2 \ldots \;] \;i_{40}\;]$
Up to 40 different printing flags can be set at a time. If the integer $i_k$ is encountered the printing flag $i_k$ is set, i.e. iprint $(i_k)=1\;\; (1 \leq i_k<999)$ . These are used to generate additional printouts by adding where needed the lines of the form

 
       if (iprint(ik).eq.1) then
           print *, ``printing something ...''
           ...
      endif

Set

• $i_1=10$ to print a total radial density relative to the centre A along the internuclear axis ( $-R_{\infty}\le z\le -R/2$ )

• $i_2=11$ to print a total radial density relative to the centre B along the internuclear axis ( $R/2\le z\le R_{\infty}$ )

See inputData routine for a list of used flags.

PRTEVERY

Format:

prtevery $i_1$ $i_2$
Subroutine pmtx can be used to output two-dimensional arrays in a tabular row-wise form with every $i_1$ -th row and $i_2$ -th column being printed (by default every 10th row and column is selected)

SCF

Format:

scf $[\;i_1\;[\;i_2\;[\;i_3\;[\;i_4\;[i_5]\;]\;]\;]\;]$

$i_1$ :

maximum number of scf iterations (default 1000); to skip the scf step set $i_1$ to a negative integer,

$i_2$ :

every $i_2$ scf iterations orbitals and potentials are saved on disk (default 20). If $i_2=0$ functions are saved on disk upon completion of the scf process. If $i_2<0$ functions are never written to disk,

$i_3$ :

if the maximum error in orbital energy is less than $10^{-i_3}$ than the scf process is terminated (the default value is 10),

$i_4$ :

if the maximum error in orbital norm is less than $10^{-i_4}$ than scf process is terminated (the default is 10),

$i_5$ :

the level of output during scf process

• $i_5=1$ - the orbital energy, the difference between its current and previous value, the normalization error and the (absolute) value of the largest overlap integral between the current orbital and all the lower lying ones of the same symmetry (the value is zero for the lowest orbitals of each symmetry) is printed for every orbital in every scf iteration
• $i_5=2$ - the orbital energy, the difference between its current and previous value and the normalization error is printed for the worst converged orbital in energy (first line) and norm (second line) in every scf iteration (default)
• $i_5=3$ - the orbital energy, the difference between its current and previous value and the normalization error is printed for the worst converged orbital in energy (first line) and norm (second line) every $i_2$ iterations. Printing of ``... multipole moment expansion coefficients (re)calculated ...'' communique is suppressed

Total energy is printed every $i_2$ iterations.

SOR

Format:

sor $[\;i_1\;[\;i_2\;]\;]$
Selects the SOR method for solving the Poisson equations for orbitals and potentials (default) and changes the value of the SOR relaxation sweeps during a single SCF cycle for orbitals ($i_1$ ) and potentials ($i_2$ ); by default $i_1=i_2=10$ .

XALPHA

Format:

xalpha $\alpha$
This label allows to change the $\alpha$ parameter of the LDA potential (see the HFS method and dft label); 2/3 is its default value.

Deprecated labels

The following labels have been replaced by others but are supported for backward compatibility with the previous version of the input data:

INITIAL (deprecated, use ORBPOT and LCAO instead)

Format:

initial $i_1\;[\;i_2\;[\;i_3\;]\;]$

$i_1$ :

determine the initial source of orbitals and potentials:

• $i_1=1$ - molecular orbitals are formed as a linear combination of hydrogenic functions on centres $A$ and $B$ (see the description of $i_3$ for further details); in the case of HF or HFS calculations Coulomb (exchange) potentials are approximated as a linear combination of Thomas-Fermi ($1/r$ ) potentials at the two centres; if method OED is chosen the potential function is approximated as a linear combination of $Z_A/r_1$ and $Z_B/r_2$ terms and the exchange potentials are set to zero

• $i_1=2$ - GAUSSIAN94/GAUSSIAN98 output is used to retrieve exponents and expansion coefficients of (uncontracted) molecular orbitals (it is assumed that the output is contained in gaussian.out and gaussian.pun files) and Coulomb and exchange potentials are initialized as in $i_1=1$ case; see routine prepGauss for more details

• $i_1=3$ - GAUSSIAN94/GAUSSIAN98 customized output is used to retrieve exponents and expansion coefficients of molecular orbitals (it is assumed that the output is contained in gaussianc.out file) and Coulomb and exchange potentials are initialized as in $i_1=1$ case; see routine prepGaussCust for details

• $i_1=5$ - initial orbitals, Coulomb and exchange potentials are retrieved from disk files (2dhf_input.orb, 2dhf_input.coul and 2dhf_input.exch, respectively) created in a previous run. Data defining the case are retrieved from a 2dhf_input.dat textfile.

• $i_1=6$ - orbitals and Coulomb potentials are retrieved from disk files and exchange potentials are initialized as in $i_1=1$ case (convenient when going from HFS to HF calculations)

• $i_1=11$ - radial Hartree-Fock orbitals for the centre A and B are retrieved from disk files 1dhf_inputA.orb and 1dhf_inputB.orb, respectively, and Coulomb and exchange potentials are initialized as in the hydrogen case (see routine initHF for details).

• $i_1=55$ - initial orbitals and potentials are retrieved from disk files but the content of a 2dhf_input.dat file is retrieved from a 2dhf_input.orb (binary) file. Use this value when reading binary data generated by older versions of the program.

$i_2$ :

specifies how exchange potentials are to be read/written and manipulated (stored in memory). The program always keeps all orbitals and Coulomb potentials in memory. If computer resources are adequate all exchange potentials can also be all kept in core memory. However, during a relaxation of a particular orbital only a fraction of exchange potentials is in fact needed. Thus all exchange potentials can be kept on disk as separate files (named fort.31, fort.32, ... during a run) and only relevant ones are being retrieved when necessary.⁵

• $i_2=0$ - read exchange potentials as separate files and write them back as separate files
• $i_2=1$ - read all exchange potentials in a file but write them out as separate files
• $i_2=2$ - read all exchange potentials separately but write them out as a single file
• $i_2=3$ - read and write exchange potentials in the form of a single file (default)

$i_3$ :

if $i_1=1$ then this parameter must be set to 1 or 2 (if omitted it is set to 1). In such a case the initialization of each of the orbitals has to be defined in terms of the linear combination of atom centered hydrogen-like functions For each orbital include a card of the following format (make sure that the order of orbitals should match the order specified under the config label):

Format:

$c_A\;\;n_A \;\;l_A \;\;\zeta_A \;\;\;\;c_B\;\;n_B\;\;l_B\;\; \zeta_B\;\;\;\;i_1\;\;\;\;[\;i_2\;]$
where

$c_A$ - relative mixing coefficient for a hydrogenic orbital on the $Z_A$ centre (real),

$n_A$ - its principle quantum number (integer)

$l_A$ - its orbital quantum number (integer)

$\zeta_A$ - the effective nuclear charge if $i_3=1$ or a screening parameter if $i_3=2$ (real)

$c_B$ - relative mixing coefficient for a hydrogenic orbital on the $Z_B$ centre (real),

$n_B$ - its principle quantum number (integer)

$l_B$ - its orbital quantum number (integer)

$\zeta_B$ - the effective nuclear charge if $i_3=1$ or a screening parameter if $i_3=2$ (real)

$i_1$ - set to 1 to freeze the orbital during scf; otherwise set to 0 (integer)

$i_2$ - a number of successive over-relaxations for a given orbital (integer); if omitted is set to 10

For other values of $i_1$ than 1 the orbital cards can be omitted but then the $i_3$ parameter must be set to 0.

FIX (deprecated, use FIXORB, FIXCOUL or FIXEXCH instead)

Format:

fix $[\;i_1\;[\;i_2\;[\;i_3\;]\;]\;]$
If $i_1$ , $i_2$ or $i_3$ are set to 1 then orbitals, Coulomb potentials or exchange potentials, respectively, are kept frozen during the scf/sor process (the respective default values are 0, 0 and 2). If $i_3=2$ then exchange potentials are relaxed only once during an scf cycle. $i_2$  and $i_3$  cannot be set to 1 if hydrogenic orbitals are used to initiate the orbitals.

Examples of input data

1. $^2$ S ground state of the Th$^{+89}$ one-electron system (see examples/oed/th+89/th+89_1s.lst).

   ../examples/oed/th+89/th+89_1s.data

2. First excited $^2$ S state of the Th$^{+89}$ one-electron system (see examples/oed/th+89/th+89_2s.lst).

   ../examples/oed/th+89/th+89_2s.data

3. Hartree-Fock ground state of the beryllium atom calculated in two consecutive steps (see examples/be/be.lst and examples/be/be-1.lst).

   ../examples/be/be.data

   ../examples/be/be-1.data

4. Hartree-Fock ground state energy of the hydrogen molecule (see examples/h2/h2.lst).

   ../examples/h2/h2.data

5. Hartree-Fock ground state of the BF molecule (see examples/bf/bf_init2.lst).

   ../examples/bf/bf_init2.data

6. HF calculations for the lowest $^3P$ state of the carbon atom (see examples/c/c.lst).

   ../examples/c/c.data

7. HF calculations for the lowest $^2P$ state of the $C^+$ ion (see examples/c+/c+.lst).

   ../examples/c+/c+.data

8. HF calculations for the lowest state of the $C_2$ molecule (see examples/c2/c2a.lst and examples/c2/c2b.lst).

   ../examples/c2/c2a.data

   ../examples/c2/c2b.data

9. HF calculations for the lowest state of the $N_2$ molecule (see examples/n2/n2.lst).

   ../examples/n2/n2.data

10. HF calculations for the lowest state of the $F_2$ molecule (see examples/f2/f2.lst).

    ../examples/f2/f2.data

11. A series of HF calculations for the the $FH$ in external static electric field.
    1. no external field (see examples/fh/fh-0.lst)

       ../examples/fh/fh-0.data

    2. field strength -0.0001 a.u. (see examples/fh/fh-m1.lst)

       ../examples/fh/fh-m1.data

    3. field strength -0.0002 a.u. (see examples/fh/fh-m2.lst)

       ../examples/fh/fh-m2.data

    4. field strength +0.0001 a.u. (see examples/fh/fh-p1.lst)

       ../examples/fh/fh-p1.data

    5. field strength 0.0002 a.u. (see examples/fh/fh-p2.lst)

       ../examples/fh/fh-p2.data

    6. DFT calculations with LDA and LYP functionals (see examples/dft/be-1.lst and examples/dft/be-2.lst).

       ../examples/dft/be-1.data

       ../examples/dft/be-2.data

    7. Two lowest states of the 2D hydrogenic harmonic potential (see examples/oed/h3/h3-1.lst and examples/oed/h3/h3-2.lst).

       ../examples/oed/h3/h3-1.data

       ../examples/oed/h3/h3-2.data

    8. Lowest state of the Krammers-Henneberger potential (see examples/oed/kh/kh.lst).

       ../examples/oed/kh/kh.data

    9. SCMC ground state calculations for the beryllium atom (see examples/scmc/be-scmc.lst).

       ../examples/scmc/be-scmc.data

Program's data files

There are several standard names used by the program to keep track of its input and output disk files. Normally the program writes out the data in the course of computations and upon its completion into the following disk files:

2dhf_output.dat (a text (ASCII) file) containing the title of a case, the time and date of its commencement, the number of mesh points, the internuclear distance, the charges of nuclei and the number of orbitals, electrons and exchange potentials (see writeDisk for details),

2dhf_output.orb (a binary file) containing molecular orbitals (in the order specified by the input data following config label) followed by their normalization factors, orbital energies, Lagrange multipliers and multipole moment expansion coefficients (see write&sstarf#star;),

2dhf_output.coul (a binary file) containing corresponding Coulomb potentials and

2dhf_output.exch (a binary file) containing all exchange potentials if the exchio [in-one$\vert$ in-many] out-one card is present

fort.31, fort.32, ... (binary files) each containing the exchange potential required for a particular pair of orbitals if the exchio in-many [out-one$\vert$ out-many] or exchio [in-one$\vert$ in-many] out-many card is present

These files can be used to restart a given case or run another with slightly modified parameters. If orbpot old card is present orbitals are retrieved from 2dhf_input.orb file, Coulomb potentials from 2dhf_input.coul and exchange potentials from 2dhf_input.exch file (or fort.31, fort 32, ...files, if exchio in-many [out-one$\vert$ out-many]). Note that there is only one set of fort files.

How to run the program?

In order to simplify the usage of the program, the xhf script (see tests/xhf) is provided to facilitate handling of the disk files. The command xhf can be envoked with one, two or three parameters. There are two basic modes of its usage:

./xhf file1 file2
runs x2dhf reading input data from file1.data file and writing text data describing the case into file2.dat file and binary data with orbitals and potentials into file2.orb, file2.coul and file2.exch files.

./xhf file1 fil2 file3
runs x2dhf reading input data from file1.data and initial orbitals and potentials from file2.dat, file2.orb, file2.coul and file2.exch files and writing resulting data into file3.dat, file3.orb, file3.coul and file3.exch files.

If, for example, be.data file contains input data for the beryllium atom (see Example 3) then

./xhf be be-1

starts and performs the first 300 scf iterations. Type

./xhf be-1 be-1 be-2

to continue calculations. In order to converge the SCF process even better increase the convergence parameters (see the scf label) and use the following command

./xhf be-1 be-2 be-1

In addition, the xhf script can be used to perform the following tasks:

./xhf stop
creats stop_x2dhf file in a current directory to stop a running program (see Section 5)

./xhf mkgauss filename
creates symbolic links gaussian.out and gaussian.pun to files filename.out and filename.pun, respectively (see e.g. Example 5).

./xhf rmgauss
removes gaussian.out and gaussian.pun files from a current directory

./xhf clean
removes *.[dat|orb|coul|exch] files

Some hints

1. The program should be easy to use provided you can start a calculation for a specific system. You should not encounter any serious problems when the system contains atoms from the first two rows of the Periodic Table. Then even the rough hydrogenic estimates of the orbitals should prove adequate and after the initial couple of dozen of iterations a smooth convergence should set in.

   If, however, a system contains more than 15-20 electrons the initial estimates of the orbitals have to be good enough to avoid divergences. Then, you have to choose the parameters of the hydrogenic orbitals carefully or perform the finite basis set calculations using the Gaussian94 to provide the initialization data for orbitals (see Example 5).

   One can also use HFS method to produce initial estimates of orbitals and Coulomb potentials. For example, to start calculations for the neon atom one can use the following input data:

   ../examples/ne/ne-hfs.data

   This input produces good enough HFS initial estimates so that the calculations can be continued at the HF level with the corresponding input:

   ../examples/ne/ne-hf.data

   One can also use HF method with some model potential, e.g. the model potential of Green, Sellin, Zachor (label POTGSZ).

2. At the very beginning set the maximum number of scf iterations to something between 20 and 50 and/or impose crude convergence criteria for the orbital energy and normalization.

3. In case of convergence problems try to perform calculations on a sparser grid. For example, the $[61\times 79;30]$ grid is sufficient to check the quality of the initial data for the Ne$_2$ system.

4. Choose smaller values of the relaxation parameters ( $1.7\leq \omega\leq 1.85$ ) to avoid divergences in the first few dozens of SCF iterations (rarely the values as small as 1.2 may be needed). Subsequently the parameters should be increased to their (near) optimum value (see the omega label and Example 8).

   It is possible to set $\omega_{pot}$ to its near-optimal value by calculating it from a semiempirical formula; see the omega label. As a rule of thumb the optimal value of the orbital relaxation parameter is somewhat smaller and, by default, is obtained by scaling the $\omega_{pot}$ value by 0.98 (see setDefaults).

5. How to stop the program gracefully during a lengthy calculation without killing the process and interrupting disk read/write operations? All you have to do is to create a (zero length) file named stop_x2dhf in a working directory by typing ./xhf stop (you can also use the Unix touch command to this end). The program stops whenever this file is detected upon the completion of a current orbital/potential relaxation.

About this document ...

This document was generated using the LaTeX2HTML translator Version 2008 (1.71)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -split 0 users-guide.tex

The translation was initiated by on 2010-05-19

Footnotes

... program¹

J.Kobus and Ch. Froese Fischer, Quasi-Relativistic Hartree-Fock program for Atoms, to be published.

... time²

The maximum number of flags can be changed by adjusting the parameter maxflags (see inputData).

... system.³

In this type of calculations convergence rates differ greatly between orbitals. Therefore, if for a given orbital the orbital energy threshold is reached it is being frozen.

... method⁴

V.V.Karasiev and E.V.Ludenia, Self-consistent multiplicative constant method for the exchange energy in density functional theory, Phys. Rev. A 65 (2002) 062510.

... necessary.⁵

A note of warning for the users of the g77 compiler. You might encounter an I/O error when trying to run cases requiring more than 70 exchange potentials. By default g77 accepts file unit numbers in the range 0-99. If you need more files to be opened you have to edit f/runtime/libI77/fio.h in the g77 source tree, changing the line: #define MXUNIT 100. Change the line so that the value of MXUNIT is defined to be at least one greater than the maximum unit number needed.