[gpaw-users] inverse_cholesky error

[Andrew Medford]

Hi all,

I am attempting to do some calculations with WC and am running into an
inverse_cholesky error. I am currently just trying to converge the energy
with respect to grid spacing, but this calculation on 2 atoms gives me the
error. I read from previous posts that changing the number of bands can
sometimes fix this, and managed to avoid the error by using less than 7
unoccupied bands. It seems that the energy is sufficiently converged with
even 2-3 unoccupied bands so I guess the problem is fixed, but I was
curious if anyone could shed some light on what is going on or if there is
a better way to avoid it. Specifically, should this approach be valid when
I go to larger unit cells or surfaces? Will 6 unoccupied bands always be
sufficient, or will I need more when I include more atoms in the cell?


AJ Medford

More details: The problem occurs with either the rmm-diis or cg
eigensolver. Changing the basis sets for the initial guess does not fix the
problem. The calculation seems to work in LCAO mode with the default number
of bands. I have attached the error file output and the script below. Any
help or advice is appreciated. Thanks,

Error File:
Traceback (most recent call last):
  File "opt_nbands_WC.py", line 80, in ?
    wc.get_potential_energy()
  File
"/nfs/slac/g/suncatfs/sw/gpawv9/install/lib/python2.4/site-packages/ase/atoms.py",
line 516, in get_potential_energy
    return self._calc.get_potential_energy(self)
  File
"/nfs/slac/g/suncatfs/sw/gpawv9/install/lib64/python2.4/site-packages/gpaw/aseinterface.py",
line 33, in get_potential_energy
    self.calculate(atoms, converge=True)
  File
"/nfs/slac/g/suncatfs/sw/gpawv9/install/lib64/python2.4/site-packages/gpaw/paw.py",
line 258, in calculate
    self.occupations):
  File
"/nfs/slac/g/suncatfs/sw/gpawv9/install/lib64/python2.4/site-packages/gpaw/scf.py",
line 46, in run
    wfs.eigensolver.iterate(hamiltonian, wfs)
  File
"/nfs/slac/g/suncatfs/sw/gpawv9/install/lib64/python2.4/site-packages/gpaw/eigensolvers/eigensolver.py",
line 72, in iterate
    wfs.orthonormalize()
  File
"/nfs/slac/g/suncatfs/sw/gpawv9/install/lib64/python2.4/site-packages/gpaw/wavefunctions/fdpw.py",
line 188, in orthonormalize
    self.overlap.orthonormalize(self, kpt)
  File
"/nfs/slac/g/suncatfs/sw/gpawv9/install/lib64/python2.4/site-packages/gpaw/overlap.py",
line 94, in orthonormalize
    self.ksl.inverse_cholesky(S_nn)
  File
"/nfs/slac/g/suncatfs/sw/gpawv9/install/lib64/python2.4/site-packages/gpaw/kohnsham_layouts.py",
line 144, in inverse_cholesky
    self._inverse_cholesky(S_NN)
  File
"/nfs/slac/g/suncatfs/sw/gpawv9/install/lib64/python2.4/site-packages/gpaw/kohnsham_layouts.py",
line 151, in _inverse_cholesky
    inverse_cholesky(S_NN)
  File
"/nfs/slac/g/suncatfs/sw/gpawv9/install/lib64/python2.4/site-packages/gpaw/utilities/lapack.py",
line 99, in inverse_cholesky
    raise RuntimeError('inverse_cholesky error: %d' % info)
RuntimeError: inverse_cholesky error: 14


Calculation Script (sorry its a little longer than necessary, but it makes
it easier to modify things):

#This script calculates the energy of a WC unit cell as a function of the
number of empty bands. The range of empty bands can be specified in the
bandrange va
riable. The script will successfully calculate the energy at nbands=-6, but
gives an inverse_cholesky error at nbands = -7.

from ase.lattice.spacegroup import crystal
from ase.visualize import view
from ase.io import *
from gpaw import GPAW,Mixer,FermiDirac
from gpaw.poisson import PoissonSolver
from gpaw.dipole_correction import DipoleCorrectionPoissonSolver
from gpaw import GPAW
import numpy as np

formula = 'WC'
a = 2.9
b = 2.85
bandrange = range(6,8) #specify the range of the number of empty bands
name = formula + '_nband_opt'

xc = 'RPBE' #Exchange-correlation functional
mode = None #Set to 'lcao' for lcao mode
dipole = False #Use dipole correction? True = yes, anything else = No
dipole_axis = 2 #Axis to apply dipole correction.
h = 0.1 #Approximate grid spacing
kpts = (4,4,4) #k-point density (x,y,z)
fermi_smearing = 0.10 #Smearing of Fermi population. 0.1 eV corresponds to
1160 K. Results are automatically extrapolated to 0 K
maxiter = 600 #Maximum iterations for electronic convergence
mixer_beta = 0.1 #Beta parameter for mixer (linear mixing coefficient).
Default = 0.1. No convergence - try 0.05
mixer_nmaxold = 5 #Number of old densities to use in mixer. Default = 5
mixer_weight = 50.0 #Weight parameter for mixer (weight used by the
metric?). Default = 50. No convergence - try 100
eigensolver = 'cg' #Eigensolver to use. Default = rmm-diis. No convergence
- try 'cg'
atomic_basis = 'sz(dzp)' #Atomic basis sets to use for initial guess.
Default is none (pseudo partial waves). No convergence = try 'sz(dzp)' or
'dzp'
convergence_energy = 0.0005 #Energy change must be less than this for the
last 3 SCF iterations. Default = 0.0005.
convergence_density = 1e-4 #Integrated absolute value of density change
must be less than this value. Default = 1e-4
convergence_eigenstates = 4e-8 #Integrated value of the square of the
residuals of KS equations should be less than this value. Default = 4e-8
convergence_bands = 'occupied' #Bands required to converge. Default =
'occupied'. Force convergence of all bands with 'all', lower n bands with
n, and all but
 last n bands with -n.

poisson_eps = 2e-10 #Convergence criterion for self-consistent Poisson
solution. Default = 2e-10. No convergence - try 1e-12.
poisson_maxiterations = 1000 #Maximum iterations for Poisson solver.
Default = 1000.
poisson_relaxation = 'J' #Relaxation method for poisson solver. Default is
'J' for Jacobian. Use 'GS' for Gauss-Seidel which will be faster for serial
calcula
tions.
poisson_stencil_accuracy = 3 #Accuracy of Laplacian stencil for KS and
Poisson equation solution. Default is 3.
mixer = Mixer(beta=mixer_beta, nmaxold = mixer_nmaxold, weight =
mixer_weight)
fermidirac = FermiDirac(fermi_smearing)
convergence = {'energy': convergence_energy,
        'density': convergence_density,
        'eigenstates': convergence_eigenstates,
        'bands': convergence_bands}
p = PoissonSolver(nn = poisson_stencil_accuracy, relax =
poisson_relaxation, eps=poisson_eps)
p.maxiter = poisson_maxiterations

if dipole == True:
    poisson = DipoleCorrectionPoissonSolver(p,dipole_axis)
else:
    poisson = p
calcargs = dict(xc=xc,
        h=h,
        kpts=kpts,
        mixer=mixer,
        nbands=0,
        maxiter=maxiter,
        occupations=fermidirac,
        convergence = convergence,
        txt = name+'_gpaw.txt')
if atomic_basis:
    calcargs['basis'] = atomic_basis
if mode:
    calcargs['mode'] = mode

traj = PickleTrajectory(name+'.traj','w')

for n in bandrange:
    calcargs['nbands'] = -1*n
    calc = GPAW(**calcargs)
    wc  = crystal('WC', [(0,0,0),(0.333,0.666,0.5)], spacegroup=187,
cellpar=[a, a, b, 90, 90, 120])
    wc.set_calculator(calc)
    wc.get_potential_energy()
    traj.write(wc)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://listserv.fysik.dtu.dk/pipermail/gpaw-users/attachments/20111129/4744173f/attachment.html