Trace of the RHF density matrix does not equal the number of electrons

[br2575]

We are trying to compute natural orbitals by direct diagonalization. As a test case, we are looking at water (see below) and are finding that the trace of the RHF density matrix (6.524993937115451) does not equal the number of electrons (10). Is this a bug or are we doing something wrong? Alternatively, is there a keyword to get these natural orbitals from a mean-field method?

from scipy import linalg 
from pyscf import scf, gto

mol         = gto.Mole()
mol.verbose = 4
mol.atom    = ''' 
O          0.00000        0.00000        0.11779
H          0.00000        0.75545       -0.47116
H          0.00000       -0.75545       -0.47116
'''
mol.basis   = 'ccpvdz'
mol.charge  = 0
mol.spin    = 0
mol.build()
nbasis = mol.nao_nr()

m = scf.RHF(mol)
m.conv_tol = 1.0E-10
m.kernel()

dm = m.make_rdm1()
print "dm trace", np.trace(dm)

evals, evecs = linalg.eigh(dm)

for f in range(nbasis): 
    print "orbital number=", f, "occupation number=", evals[f]
    print evecs[:,f]

[gkc1000]

The DM is being computed in the AO matrix, so the correct norm is trace(DS) where S is the AO overlap metric.

Similarly, when diagonalizing to get the natural orbitals, you have to include the AO overlap matrix, i.e. D C = S C e.

On Tue, Jun 11, 2019 at 12:46 PM br2575 notifications@github.com wrote:

We are trying to compute natural orbitals by direct diagonalization. As a test case, we are looking at water (see below) and are finding that the trace of the RHF density matrix (6.524993937115451) does not equal the number of electrons (10). Is this a bug or are we doing something wrong? Alternatively, is there a keyword to get these natural orbitals from a mean-field method?

from scipy import linalg from pyscf import scf, gto

mol = gto.Mole() mol.verbose = 4 mol.atom = ''' O 0.00000 0.00000 0.11779 H 0.00000 0.75545 -0.47116 H 0.00000 -0.75545 -0.47116 ''' mol.basis = 'ccpvdz' mol.charge = 0 mol.spin = 0 mol.build() nbasis = mol.nao_nr()

m = scf.RHF(mol) m.conv_tol = 1.0E-10 m.kernel()

dm = m.make_rdm1() print "dm trace", np.trace(dm)

evals, evecs = linalg.eigh(dm)

for f in range(nbasis): print "orbital number=", f, "occupation number=", evals[f] print evecs[:,f]

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/pyscf/pyscf/issues/359?email_source=notifications&email_token=AAN5NRHZEFWUMWO7MZNWD5LPZ76K3A5CNFSM4HXCZAU2YY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4GY4ZHTQ, or mute the thread https://github.com/notifications/unsubscribe-auth/AAN5NRHUTYEDN22SVPPOUATPZ76K3ANCNFSM4HXCZAUQ .

[br2575]

Thank you for clarifying that point.

We found the eig function for solving the generalized eigenvalue problem at https://sunqm.github.io/pyscf/scf.html and tried to use it as below. Did we follow your instructions to find the natural orbitals in the AO basis properly?

To give more context, we are trying to generate the natural orbitals to more easily identify which orbitals have significant d character so that we may add them to the initial active space for a CASSCF calculation. Orca has an automated way to select the second d shell; does pyscf have something similar?

thanks,

import numpy as np
from scipy import linalg 
from pyscf import scf, gto

mol         = gto.Mole()
mol.verbose = 4
mol.atom    = ''' 
O          0.00000        0.00000        0.11779
H          0.00000        0.75545       -0.47116
H          0.00000       -0.75545       -0.47116
'''
mol.basis   = 'ccpvdz'
mol.charge  = 0
mol.spin    = 0
mol.build()
nbasis = mol.nao_nr()

m = scf.RHF(mol)
m.conv_tol = 1.0E-10
m.kernel()

s1e = m.get_ovlp(mol)
dm = m.make_rdm1()
evals_eig , evecs_eig = m.eig( dm, s1e  )

for f in range(nbasis):
    print "orbital number=", f
    print "evec in AO basis", evecs_eig[:,f]

[gkc1000]

Yes, this will find the NO, unfortunately, you are finding the NO of restricted HF theory. These will just be the HF occupied or virtual orbitals.

There are a couple of simple ways to find a second d-shell orbital. For example, by localizing the DFT occupied and virtuals and looking at them, you can often identify the 4d orbitals. However, the easiest way would be to use the automated valence active space construction, see e.g.

https://github.com/pyscf/pyscf/blob/master/examples/mcscf/43-avas.py

described in more detail in the paper by Sayfutyarova et al.

You can insert "M 4d" to add the 4d shell to the space of the atom called M. It should just return an active space containing the 4d as well as whatever other valence orbitals you specify. I assume you want to do this for a TM complex, not the water molecule in your example!

On Wed, Jun 12, 2019 at 1:37 PM br2575 notifications@github.com wrote:

Thank you for clarifying that point.

We found the eig function for solving the generalized eigenvalue problem at https://sunqm.github.io/pyscf/scf.html and tried to use it as below. Did we follow your instructions to find the natural orbitals in the AO basis properly?

To give more context, we are trying to generate the natural orbitals to more easily identify which orbitals have significant d character so that we may add them to the initial active space for a CASSCF calculation. Orca has an automated way to select the second d shell; does pyscf have something similar?

thanks,

import numpy as np from scipy import linalg from pyscf import scf, gto

mol = gto.Mole() mol.verbose = 4 mol.atom = ''' O 0.00000 0.00000 0.11779 H 0.00000 0.75545 -0.47116 H 0.00000 -0.75545 -0.47116 ''' mol.basis = 'ccpvdz' mol.charge = 0 mol.spin = 0 mol.build() nbasis = mol.nao_nr()

m = scf.RHF(mol) m.conv_tol = 1.0E-10 m.kernel()

s1e = m.get_ovlp(mol) dm = m.make_rdm1() evals_eig , evecs_eig = m.eig( dm, s1e )

for f in range(nbasis): print "orbital number=", f print "evec in AO basis", evecs_eig[:,f]

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/pyscf/pyscf/issues/359?email_source=notifications&email_token=AAN5NRDSSSQSDFRYCCBKB43P2FM75A5CNFSM4HXCZAU2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODXRXKHY#issuecomment-501445919, or mute the thread https://github.com/notifications/unsubscribe-auth/AAN5NRB46E25FTLDQZO34PTP2FM75ANCNFSM4HXCZAUQ .

[br2575]

Hi Garnet, thanks. And yes we intend to study TM complexes!

It appears the MINAO basis doesn't have 4d AOs. In the paper by Sayfutyarova et al the ANO-RCC and other basis sets are suggested as an alternative, if 4d AOs are desired.

Could you please explain how we can use such alternatives in PySCF? Is it as simple as replacing MINAO = getattr(__config__, 'mcscf_avas_minao', 'minao') with something like ANORCC = getattr(__config__, 'mcscf_avas_ano-rcc', 'ano-rcc'), and then using: norb, ne_act, orbs = avas.avas(mf, ao_labels, canonicalize=False, minao=ANORCC)?

[sunqm]

ANORCC in this case needs to be truncated. Otherwise the entire ANO set will be used as the reference basis to generate active space.

There are many ways to truncate the basis. The simplest one is https://github.com/pyscf/pyscf/blob/master/examples/gto/04-input_basis.py#L98 . Actually, all the methods to input basis as shown in that example can be used to define the avas reference basis

[jshee]

Thanks Qiming. So if we want to combine the minao basis with the 4d functions from the ano basis, can't we use something like the following? basis = ['minao', 'ano@4d']

When printing with gto.spheric_labels(mol), it does add the 4d shell, but also 5-7d:
['0 Ti 1s ', '0 Ti 2s ', '0 Ti 3s ', '0 Ti 4s ', '0 Ti 2px ', '0 Ti 2py ', '0 Ti 2pz ', '0 Ti 3px ', '0 Ti 3py ', '0 Ti 3pz ', '0 Ti 3dxy ', '0 Ti 3dyz ', '0 Ti 3dz^2 ', '0 Ti 3dxz ', '0 Ti 3dx2-y2', '0 Ti 4dxy ', '0 Ti 4dyz ', '0 Ti 4dz^2 ', '0 Ti 4dxz ', '0 Ti 4dx2-y2', '0 Ti 5dxy ', '0 Ti 5dyz ', '0 Ti 5dz^2 ', '0 Ti 5dxz ', '0 Ti 5dx2-y2', '0 Ti 6dxy ', '0 Ti 6dyz ', '0 Ti 6dz^2 ', '0 Ti 6dxz ', '0 Ti 6dx2-y2', '0 Ti 7dxy ', '0 Ti 7dyz ', '0 Ti 7dz^2 ', '0 Ti 7dxz ', '0 Ti 7dx2-y2']

Am I misinterpreting this, or is the @4d not working as expected?

thanks, James

[sunqm]

In the current version, "@4d" will take 4 sets of d orbitals from ANO basis, not just the 4d shell.

I think it may be useful to have a way to select a specific basis from a basis set. I opened a new issue #376 about this feature

[bogdanoff]

Hi James, You can try to specify minao option for avas without changing basis and other things norb, ne_act, orbs = avas.avas(mf, ao_labels, canonicalize=False, minao='ano')