Source code for pyscf.scf.diis

#!/usr/bin/env python
# Copyright 2014-2018 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#

"""
DIIS
"""

from functools import reduce
import numpy
import scipy.linalg
import scipy.optimize
from pyscf import lib
from pyscf.lib import logger

DEBUG = False

# J. Mol. Struct. 114, 31-34
# PCCP, 4, 11
# GEDIIS, JCTC, 2, 835
# C2DIIS, IJQC, 45, 31

# error vector = SDF-FDS
# error vector = F_ai ~ (S-SDS)*S^{-1}FDS = FDS - SDFDS ~ FDS-SDF in converge
class CDIIS(lib.diis.DIIS):
    def __init__(self, mf=None, filename=None):
        lib.diis.DIIS.__init__(self, mf, filename)
        self.rollback = False
        self.space = 8

    def update(self, s, d, f, *args, **kwargs):
        errvec = get_err_vec(s, d, f)
        logger.debug1(self, 'diis-norm(errvec)=%g', numpy.linalg.norm(errvec))
        xnew = lib.diis.DIIS.update(self, f, xerr=errvec)
        if self.rollback > 0 and len(self._bookkeep) == self.space:
            self._bookkeep = self._bookkeep[-self.rollback:]
        return xnew

    def get_num_vec(self):
        if self.rollback:
            return self._head
        else:
            return len(self._bookkeep)

SCFDIIS = SCF_DIIS = DIIS = CDIIS

[docs]def get_err_vec(s, d, f): '''error vector = SDF - FDS''' if isinstance(f, numpy.ndarray) and f.ndim == 2: sdf = reduce(numpy.dot, (s,d,f)) errvec = sdf.T.conj() - sdf elif isinstance(f, numpy.ndarray) and f.ndim == 3 and s.ndim == 3: errvec = [] for i in range(f.shape[0]): sdf = reduce(numpy.dot, (s[i], d[i], f[i])) errvec.append((sdf.T.conj() - sdf)) errvec = numpy.vstack(errvec) elif f.ndim == s.ndim+1 and f.shape[0] == 2: # for UHF nao = s.shape[-1] s = lib.asarray((s,s)).reshape(-1,nao,nao) return get_err_vec(s, d.reshape(s.shape), f.reshape(s.shape)) else: raise RuntimeError('Unknown SCF DIIS type') return errvec
[docs]class EDIIS(lib.diis.DIIS): '''SCF-EDIIS Ref: JCP 116, 8255 ''' def update(self, s, d, f, mf, h1e, vhf): if self._head >= self.space: self._head = 0 if not self._buffer: shape = (self.space,) + f.shape self._buffer['dm' ] = numpy.zeros(shape, dtype=f.dtype) self._buffer['fock'] = numpy.zeros(shape, dtype=f.dtype) self._buffer['etot'] = numpy.zeros(self.space) self._buffer['dm' ][self._head] = d self._buffer['fock'][self._head] = f self._buffer['etot'][self._head] = mf.energy_elec(d, h1e, vhf)[0] self._head += 1 ds = self._buffer['dm' ] fs = self._buffer['fock'] es = self._buffer['etot'] etot, c = ediis_minimize(es, ds, fs) logger.debug1(self, 'E %s diis-c %s', etot, c) fock = numpy.einsum('i,i...pq->...pq', c, fs) return fock
def ediis_minimize(es, ds, fs): nx = es.size nao = ds.shape[-1] ds = ds.reshape(nx,-1,nao,nao) fs = fs.reshape(nx,-1,nao,nao) df = numpy.einsum('inpq,jnqp->ij', ds, fs).real diag = df.diagonal() df = diag[:,None] + diag - df - df.T def costf(x): c = x**2 / (x**2).sum() return numpy.einsum('i,i', c, es) - numpy.einsum('i,ij,j', c, df, c) def grad(x): x2sum = (x**2).sum() c = x**2 / x2sum fc = es - 2*numpy.einsum('i,ik->k', c, df) cx = numpy.diag(x*x2sum) - numpy.einsum('k,n->kn', x**2, x) cx *= 2/x2sum**2 return numpy.einsum('k,kn->n', fc, cx) if DEBUG: x0 = numpy.random.random(nx) dfx0 = numpy.zeros_like(x0) for i in range(nx): x1 = x0.copy() x1[i] += 1e-4 dfx0[i] = (costf(x1) - costf(x0))*1e4 print((dfx0 - grad(x0)) / dfx0) res = scipy.optimize.minimize(costf, numpy.ones(nx), method='BFGS', jac=grad, tol=1e-9) return res.fun, (res.x**2)/(res.x**2).sum()
[docs]class ADIIS(lib.diis.DIIS): ''' Ref: JCP, 132, 054109 ''' def update(self, s, d, f, mf, h1e, vhf): if self._head >= self.space: self._head = 0 if not self._buffer: shape = (self.space,) + f.shape self._buffer['dm' ] = numpy.zeros(shape, dtype=f.dtype) self._buffer['fock'] = numpy.zeros(shape, dtype=f.dtype) self._buffer['dm' ][self._head] = d self._buffer['fock'][self._head] = f ds = self._buffer['dm' ] fs = self._buffer['fock'] fun, c = adiis_minimize(ds, fs, self._head) if self.verbose >= logger.DEBUG1: etot = mf.energy_elec(d, h1e, vhf)[0] + fun logger.debug1(self, 'E %s diis-c %s ', etot, c) fock = numpy.einsum('i,i...pq->...pq', c, fs) self._head += 1 return fock
def adiis_minimize(ds, fs, idnewest): nx = ds.shape[0] nao = ds.shape[-1] ds = ds.reshape(nx,-1,nao,nao) fs = fs.reshape(nx,-1,nao,nao) df = numpy.einsum('inpq,jnqp->ij', ds, fs).real d_fn = df[:,idnewest] dn_f = df[idnewest] dn_fn = df[idnewest,idnewest] dd_fn = d_fn - dn_fn df = df - d_fn[:,None] - dn_f + dn_fn def costf(x): c = x**2 / (x**2).sum() return (numpy.einsum('i,i', c, dd_fn) * 2 + numpy.einsum('i,ij,j', c, df, c)) def grad(x): x2sum = (x**2).sum() c = x**2 / x2sum fc = 2*dd_fn fc+= numpy.einsum('j,kj->k', c, df) fc+= numpy.einsum('i,ik->k', c, df) cx = numpy.diag(x*x2sum) - numpy.einsum('k,n->kn', x**2, x) cx *= 2/x2sum**2 return numpy.einsum('k,kn->n', fc, cx) if DEBUG: x0 = numpy.random.random(nx) dfx0 = numpy.zeros_like(x0) for i in range(nx): x1 = x0.copy() x1[i] += 1e-4 dfx0[i] = (costf(x1) - costf(x0))*1e4 print((dfx0 - grad(x0)) / dfx0) res = scipy.optimize.minimize(costf, numpy.ones(nx), method='BFGS', jac=grad, tol=1e-9) return res.fun, (res.x**2)/(res.x**2).sum()