14. DFT for Periodic Systems (GAPSS)

NOTE: A user is working with a prototype version of the GAPSS code. Encountered problems should be reported to J.E. Jaffe, john.jaffe@pnl.gov, or M. Gutowski, maciej.gutowski@pnl.gov.

The NWChem density functional theory module for periodic systems (GAPSS) uses the Gaussian basis set approach to compute electron densities and Kohn-Sham orbitals in the:

• local density approximation (LDA),
• non-local density approximation (NLDA),
• local spin-density approximation (LSD),
• non-local spin-density approximation (NLSD).

GAPSS input is provided using the compound GAPSS directive

  GAPSS
    ...
  END

The actual DFT calculation will be performed when the input module encounters the TASK directive (Section 5.10).

  TASK GAPSS

There are sub-directives which allow for customized application; those currently provided as options for the GAPSS module are:

  CORRELATION [(vwn||lyp||perdew86||pbe96) default vwn]
  EXCHANGE [(slater||becke88||pbe96) default slater]
  VXCACC [(low||medium||high||extrahigh||superhigh) default high]
  WEIGHT [(lin||becke||delley) default lin]
  MONKHORST <integer monkhorst default 3>
  ISCREENLAT <integer iscreenlat default 7>
  IC1C3CONV <integer ic1c3conv default 7>
  IC2C4CONV <integer ic2c4conv default 6>
  INTACC [(standard||high) default standard]
  SCFTYPE [(conventional||direct) default conventional]
  IPOL <integer ipol default 1>
  ITRSCF <integer itrscf default 30>
  SCFCONV <integer scfconv default 6>
  MIXRATE <integer mixrate default 30>
  NCYCMIX <integer ncycmix default 0>
  NEWMIXRATE <integer newmixrate default 0>
  LEVELSHIFT <real levelshift default 0.0>
  NCLSHIFT <integer nclshift default 0>
  GAPRESTART [(intsok||noints) default intsok] \
             [(covecs||atomic) default covecs]

The following sections describe these keywords and optional sub-directives that can be specified for a GAPSS calculation in NWChem.

14.1 Specification of Basis Sets for the GAPSS Module

The GAPSS module requires a basis set for the Kohn-Sham crystal orbitals. This basis set must be in the default basis set named "ao basis", or it must be assigned to this default name using the SET directive (see Section 5.7). The valence and polarization shells of standard molecular AO basis sets usually require some modification for solids often by deleting or increasing exponents that are $\leq$ 0.2 au$^{-2}$ initially. Only s, p, and 6d type functions are allowed at present.

The formal scaling of the computation is reduced by choosing to use auxiliary Gaussian basis sets to fit the charge density (CD). In addition to the basis set for the Kohn-Sham orbitals, the charge density fitting basis set must also be specified in the input directives for the GAPSS module. This basis set is used for the evaluation of the Coulomb potential in the Dunlap scheme^14.1, slightly modified for periodic systems^14.2. The charge density fitting basis set must have the name "cd basis". This can be the actual name of a basis set, or a basis set can be assigned this name using the SET directive, as described in Section 5.7. The molecular CD basis sets can often be used unchanged for periodic systems.

For the GAPSS module, the input options for defining the basis sets in a given calculation can be summarized as follows;

• "ao basis" - Kohn-Sham crystal orbitals; required for all calculations
• "cd basis" - charge density fitting basis set; required for all calculations

14.2 CORRELATION and EXCHANGE -- Exchange-Correlation Potentials

  CORRELATION [(vwn||lyp||perdew86||pbe96) default vwn]
  EXCHANGE [(slater||becke88||pbe96) default slater]

The user has the option of specifying the exchange-correlation treatment in the GAPSS Module. The default exchange-correlation functional is defined as the local density approximation (LDA) for closed shell systems and its counterpart the local spin-density (LSD) approximation for open shell systems. Within this approximation the exchange functional is the Slater $\rho^{1/3}$ functional (from J.C. Slater, Quantum Theory of Molecules and Solids, Vol. 4: The Self-Consistent Field for Molecules and Solids (McGraw-Hill, New York, 1974)), and the correlation functional is the Vosko-Wilk-Nusair (VWN) functional (functional V) (S.J. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 58, 1200 (1980)). The parameters used in this formula are obtained by fitting to the Ceperley and Alder^14.3Quantum Monte-Carlo solution of the homogeneous electron gas.

These defaults can be invoked explicitly by specifying the following keywords within the GAPSS module input directive,

  CORRELATION vwn
  EXCHANGE slater

Several alternative exchange and correlation functionals are available to the user. The following sections describe these options.

14.2.1 Optional Functionals

There are two exchange functionals in addition to the default Slater exchange functional. These are the Becke gradient-corrected functional (see A.D. Becke, J. Chem. Phys. 88, 3098 (1988)), and the Perdew, Burke, Ernzerhof generalized gradient approximation to the exchange-correlation functional (see J.P. Perdew, K. Burke, M. Ernzerhof,  Phys. Rev. Lett. 77, 3865 (1996)).

The Becke gradient-corrected functional is invoked by specifying the input line

   EXCHANGE becke88

The Perdew, Burke, Ernzerhof exchange functional is invoked by specifying the input line

   EXCHANGE pbe96

There are three correlation functionals in addition to the default vwn correlation functional. These are the LYP gradient-corrected functional (C. Lee, W. Yang and R. G. Parr, Phys. Rev. B 37, 785 (1988)), Perdew86 gradient-corrected functional (J. P. Perdew, Phys. Rev. B 33, 8822 (1986)), and the Perdew, Burke, Ernzerhof generalized gradient approximation to the exchange-correlation functional (J.P. Perdew, K. Burke, M. Ernzerhof,  Phys. Rev. Lett. 77, 3865 (1996)).

The LYP gradient-corrected functional is invoked by specifying the input line

   CORRELATION lyp

The Perdew86 gradient-corrected functional is invoked by specifying the input line

   CORRELATION perdew86

The Perdew, Burke, Ernzerhof correlation functional is invoked by specifying the input line

   CORRELATION pbe96

14.3 Numerical Integration

  VXCACC [(low||medium||high||extrahigh||superhigh) default high]

A numerical integration is necessary for the evaluation of the exchange-correlation contribution to the total energy and the Fock matrices. The user can specify the level of accuracy with the keywords; low, medium, high, extrahigh, superhigh. The default is high which corresponds to accuracy of ca. 1e-4 for the XC potential and energy. Each level of higher accuracy improves the accuracy of integration by a factor $\approx$ 10 and increases the cost by a factor of $\approx$ 3.

  WEIGHT [(lin||becke||delley) default lin]

The three-dimensional integration is reduced to a sum of one-center, atomic-like integrations using either the Lin and Hess scheme (Z. Lin, J.E. Jaffe, A.C. Hess, J. Phys. Chem. A 103, 2117 (1999))

  WEIGHT lin

or the Becke scheme (A.D. Becke, J. Chem. Phys. 88, 2547 (1988))

  WEIGHT becke

or the Delley scheme (B. Delley, J. Chem. Phys. 92, 508 (1990))

  WEIGHT delley

14.4 Summations and Integrations for a Periodic System

  MONKHORST <integer monkhorst default 3>
  ISCREENLAT <integer iscreenlat default 7>
  IC1C3CONV <integer ic1c3conv default 7>
  IC2C4CONV <integer ic2c4conv default 6>
  INTACC [(standard||high) default standard]

The wave vectors in the first Brillouin zone are sampled following the Monkhorst-Pack scheme^14.4 MONKHORST = 5 is recommended for systems with 1-2 atoms per primitive unit cell, smaller values may be used for larger cells. Metals at present require MONKHORST $\geq$ 10. The convergence of the C$^1$-C$^3$ and C$^2$-C$^4$ series is controlled by the IC1C3CONV and IC2C4CONV directives, respectively, see J.E. Jaffe, A.C. Hess, J. Chem. Phys. 105, 10983 (1996), and the accuracy of H(k) and S(k) with respect to the summation over unit cells is controlled by the ISCREENLAT directive. ISCREENLAT = 7 produces errors $\approx$ 10$^{-6}$ E$_h$, with each higher level reducing this error by ca. 1 order of magnitude at little extra cost. IC1C3CONV = ISCREENLAT and IC2C4CONV = ISCREENLAT - 1 are recommended. Finally, the accuracy of two-electron three-center integrals is controlled by the directive INTACC.

14.5 SCF Iterative Procedure

  SCFTYPE [(conventional||direct) default conventional]
  IPOL <integer ipol default 1>
  ITRSCF <integer itrscf default 30>
  SCFCONV <integer scfconv default 6>
  MIXRATE <integer mixrate default 30>
  NCYCMIX <integer ncycmix default 0>
  NEWMIXRATE <integer newmixrate default 0>
  LEVELSHIFT <real levelshift default 0.0>
  NCLSHIFT <integer nclshift default 0>
  GAPRESTART [(intsok||noints) default intsok] \
             [(covecs||atomic) default covecs]

The GAPSS code can work in the conventional mode, with integrals calculated once and written to disk

  SCFTYPE conventional

or in the direct mode, with integrals evaluated at every SCF cycle.

  SCFTYPE direct

The direct directive is mandatory in parallel runs.

Both closed-shell systems

  IPOL 1

and open-shell systems

  IPOL 2

may be studied with the GAPSS code.

The default optimization in the GAPSS module is to iterate on the Kohn-Sham equations for a specified number of iterations (default 30). The directive that controls the number of iterations is ITRSCF, and has the following general form,

   ITRSCF <integer iterations default 30>

The optimization procedure will stop when the specified number of iterations is reached or convergence is met. Convergence is satisfied when the total energy at iteration N, iteration N-1, and iteration N-2 differ by a value less than some value (the default is 1e-6). This value can be modified using the directive

  SCFCONV <integer scfconv default 6>

The final $\Delta E$ is often almost one order of magnitude below the cutoff due to the requirement of two successive $\Delta E$'s being below it. A user may control convergence of the self-consistent field procedure by mixing density matrices from the previous and current SCF cycle

  MIXRATE <integer mixrate default 30>
  NCYCMIX <integer ncycmix default 0>
  NEWMIXRATE <integer newmixrate default 0>

The directive mixrate stands for the precentage of the new density, ncycmix tells for how many SCF cycles this mixing should be applied, and newmixrate tells what is the percentage of the new density matrix after ncycmix cycles. If NCYCMIX = 0 then the mixing rate is always MIXRATE (not NEWMIXRATE). NCYCMIX = 4 and NEWMIXRATE = 55 are recommended for most problems.

A user may also apply ``level shifting''^14.5 to ensure convergence of the SCF procedure in the cases when the HOMO-LUMO separation (band gap) is small

  LEVELSHIFT <real levelshift default 0.0>
  NCLSHIFT <integer nclshift default 0>

and levelshift defines the amount of shift applied to the diagonal elements of the unoccupied block of the Fock matrix whereas nclshift tells for how many cycles the level shifting will be applied.

Finally, the GAPSS calculation may be restarted in the case of a conventional run.

  GAPRESTART [(intsok||noints) default intsok] \
             [(covecs||atomic) default covecs]

If the integrals are still available, the intsok directive should be applied. Otherwise use noints. When restarting, one may use a density formed from crystal orbitals as an initial density (directive covecs) or a superposition of atomic densities (directive atomic). Runs may be restarted (with covecs and intsok) when the SCF process was interrupted or when one needs to change the SCF conditions or the MONKHORST settings. Runs can be restarted with covecs and noints when IC1C3CONV or IC2C4CONV are changed, since this will change the values of some integrals but not their number or orders. Restarting with atomic and intsok may succeed when the SCF process has failed to converge previously. Restart is not possible at present when ISCREENLAT, AO or CD basis sets, or geometry are changed as these result in different sets of integrals being computed and stored.