moldft: virtual orbitals are broken

[lratcliff]

No description provided.

[robertjharrison]

With canonical orbitals, the problem is due to use of symmetric orthogonalization (via Newton iteration) which does not respect the priority of occupied over virtual orbitals.

In addition, localized orbitals need to modify the Lagrange multipliers in BSH update (but am not clear how).

I will make a fix for canonical orbitals. Long term I think we need a separate solver for virtuals.

[lratcliff]

I've been having a quick look at this for canonical and I think it's more than just the orthogonalization causing problems - I did a bit of a hack to make the orthogonalization asymmetric (not very efficient, but it more or less works), and that seemed to help things for some systems but not others, i.e. it gets closer to the right solution than before but still never converges. I think we might have a similar problem going on in get_fock_transformation/diag_fock_matrix where there's also some mixing between occupied and virtual, but I've yet to figure out exactly what we should be doing differently.

[robertjharrison]

The mixing in the diagonalization should be good for convergence but to
debug that just zero out the fock matrix elements that mix occupied and
virtual orbitals.
On Jul 22, 2015 5:36 PM, "lratcliff" notifications@github.com wrote:

I've been having a quick look at this for canonical and I think it's more
than just the orthogonalization causing problems - I did a bit of a hack to
make the orthogonalization asymmetric (not very efficient, but it more or
less works), and that seemed to help things for some systems but not
others, i.e. it gets closer to the right solution than before but still
never converges. I think we might have a similar problem going on in
get_fock_transformation/diag_fock_matrix where there's also some mixing
between occupied and virtual, but I've yet to figure out exactly what we
should be doing differently.

—
Reply to this email directly or view it on GitHub
#127 (comment)
.

[robertjharrison]

What system (molecule and functional) are you testing on?

On Wed, Jul 22, 2015 at 8:51 PM, Robert Harrison rjharrison@gmail.com
wrote:

The mixing in the diagonalization should be good for convergence but to
debug that just zero out the fock matrix elements that mix occupied and
virtual orbitals.
On Jul 22, 2015 5:36 PM, "lratcliff" notifications@github.com wrote:

I've been having a quick look at this for canonical and I think it's more
than just the orthogonalization causing problems - I did a bit of a hack to
make the orthogonalization asymmetric (not very efficient, but it more or
less works), and that seemed to help things for some systems but not
others, i.e. it gets closer to the right solution than before but still
never converges. I think we might have a similar problem going on in
get_fock_transformation/diag_fock_matrix where there's also some mixing
between occupied and virtual, but I've yet to figure out exactly what we
should be doing differently.

—
Reply to this email directly or view it on GitHub
#127 (comment)
.

Robert J. Harrison
tel: 865-272-9262

[lratcliff]

I tried that, and it fixed things in some cases, but in a couple of cases it made convergence much slower and in another case nothing changed. I've been testing an assortment of small molecules using LDA.

[robertjharrison]

Using nwchem in a large basis or another code what is the virtual orbital
eigenvalue? Is it negative? Even if it is negative do you need a larger
box to represent it?
On Jul 27, 2015 11:05 AM, "lratcliff" notifications@github.com wrote:

I tried that, and it fixed things in some cases, but in a couple of cases
it made convergence much slower and in another case nothing changed. I've
been testing an assortment of small molecules using LDA.

—
Reply to this email directly or view it on GitHub
#127 (comment)
.

[fbischoff]

Laura,

in our experience LDA virtual orbitals should be unbound (especially for small molecules). We did TDA calculations on helium and we used box sizes of about 1000 bohr and it was still not converged. If you use HF it's better, but the higher virtuals still may be unbound. What exactly do you want to do with the virtuals?

Florian

On Jul 27, 2015, at 5:42 PM, robertjharrison notifications@github.com wrote:

Using nwchem in a large basis or another code what is the virtual orbital
eigenvalue? Is it negative? Even if it is negative do you need a larger
box to represent it?
On Jul 27, 2015 11:05 AM, "lratcliff" notifications@github.com wrote:

I tried that, and it fixed things in some cases, but in a couple of cases
it made convergence much slower and in another case nothing changed. I've
been testing an assortment of small molecules using LDA.

�
Reply to this email directly or view it on GitHub
#127 (comment)
.

�
Reply to this email directly or view it on GitHub.

[lratcliff]

I checked each system with BigDFT first to make sure I'm only looking for the negative eigenvalues (i.e. just the first few states). I haven't modified the box sizes so that's probably the next thing to check, but I think moldft usually uses larger box sizes than BigDFT. I'm just looking to calculate dipole matrix elements between core and virtual states as a first approximation of core spectra.

[fbischoff]

Unfortunately I'm not familiar with core spectra.. Do you need actual virtual orbitals for that or do you only need a complete set of orthogonal functions? In this case you could just generate a set of gaussian functions and project out the occupied space. As for virtuals we also had some problems with generating suitable guess functions. In the end we would compute a large Fock matrix and diagonalize that to get a guess.

Florian

On Jul 27, 2015, at 5:57 PM, lratcliff notifications@github.com wrote:

I checked each system with BigDFT first to make sure I'm only looking for the negative eigenvalues (i.e. just the first few states). I haven't modified the box sizes so that's probably the next thing to check, but I think moldft usually uses larger box sizes than BigDFT. I'm just looking to calculate dipole matrix elements between core and virtual states as a first approximation of core spectra.

�
Reply to this email directly or view it on GitHub.

[lratcliff]

There might be a more sophisticated approach that doesn't need them explicitly, but the method I'm using needs the actual orbitals and their energies. For the moment my aim is just to use it as an application of the mixed all electron/pseudopotential implementation, so I wasn't looking to implement anything more advanced.
I did try using a larger basis set for the initial guess, but that sometimes breaks things even for a normal moldft calculation without virtuals, so either I made a mistake importing the basis sets, or this might be a general robustness problem with the solver. I also tried increasing the box size and that made no difference.

[kottmanj]

I changed a few lines in the TDA response code that it is now able to
solve for (bounded) virtual orbitals.
Tests on cations (Li+, B+, Na+) converged to the right values.

For which systems are you trying to get the virtuals ?

I have added an example input for the B+ cation in src/examples/
The response code is executed over src/examples/tdhf (libraries at
chem/TDA.h)

For now this does only CIS but since for the virtuals only the XC
potential is needed it should be enough to exchange the K operator with
the vxc potential of moldft.

Best wishes from Berlin,
Jakob

[robertjharrison]

I think that even if BigDFT uses a smaller box size than MADNESS the way
that BigDFT iterates to convergence will cause it to not really experience
convergence problems, whereas MADNESS will constantly be trying to find a
solution with the correct decay only to truncate the solution when it does
not fit in the box.

What is the BigDFT eigenvalue for the virtual?

What molecule is this? Can you please post the MADNESS input and output?

On Mon, Jul 27, 2015 at 11:57 AM, lratcliff notifications@github.com
wrote:

I checked each system with BigDFT first to make sure I'm only looking for
the negative eigenvalues (i.e. just the first few states). I haven't
modified the box sizes so that's probably the next thing to check, but I
think moldft usually uses larger box sizes than BigDFT. I'm just looking to
calculate dipole matrix elements between core and virtual states as a first
approximation of core spectra.

—
Reply to this email directly or view it on GitHub
#127 (comment)
.

Robert J. Harrison
tel: 865-272-9262

[lratcliff]

Thanks for the input Jakob, I'll have a look at your code to see if there's anything obvious we're missing. Mostly though it already seems to be working already for very small molecules, so I'm not sure how much we can conclude anything from tests on cations. It seems to be that the larger the system, the worse the problem, at least from my small sampling of a dozen or so molecules.

One example where things are broken even with the asymmetric orthogonalization is SO2. I'm just trying to calculate the LUMO, for which BigDFT gives an energy of -0.18. Without the virtual state, moldft converges in 10 iterations, whereas with the virtual state it's still not converged after 80 iterations. I tried nearly doubling the cell size, and it made no difference. Here's the input file I'm using:

dft
xc lda
protocol 1e-4 1e-6
maxiter 40
aobasis 6-31g
canon
nvalpha 1
nvbeta 1
end

geometry
units angstrom
S .000000 .000000 .370268
O .000000 1.277617 -.370268
O .000000 -1.277617 -.370268
end

[robertjharrison]

Can you send the output? Have the occupied states converged but not the
virtual? Or are the occupied screwed up?

On Wed, Aug 5, 2015 at 5:06 PM, lratcliff notifications@github.com wrote:

Thanks for the input Jakob, I'll have a look at your code to see if
there's anything obvious we're missing. Mostly though it already seems to
be working already for very small molecules, so I'm not sure how much we
can conclude anything from tests on cations. It seems to be that the larger
the system, the worse the problem, at least from my small sampling of a
dozen or so molecules.

One example where things are broken even with the asymmetric
orthogonalization is SO2. I'm just trying to calculate the LUMO, for which
BigDFT gives an energy of -0.18. Without the virtual state, moldft
converges in 10 iterations, whereas with the virtual state it's still not
converged after 80 iterations. I tried nearly doubling the cell size, and
it made no difference. Here's the input file I'm using:

dft
xc lda
protocol 1e-4 1e-6
maxiter 40
aobasis 6-31g
canon
nvalpha 1
nvbeta 1
end

geometry
units angstrom
S .000000 .000000 .370268
O .000000 1.277617 -.370268
O .000000 -1.277617 -.370268
end

—
Reply to this email directly or view it on GitHub
#127 (comment)
.

Robert J. Harrison
tel: 865-272-9262

[kottmanj]

I have included the lda functional for the calculation of virtuals in
the tdhf code.
The LUMO of the SO2 molecule converged in 25 iterations (dconv 1.e-4,
econv 1.e-6, thresh 1.e-6, input in the appendix).
This was also possible with a small guess for the LUMO so the guess size
should not be the problem.

I guess it is more stable to first converge all occupied states and then
do the virtuals (this is what the modified tdhf code basically does),
this should also prevent 'mixing' of occupied and virtual states.

Here is my input for the SO2 calculations

plot
plane x2 x3
zoom 1.0
npoints 50
origin 0.5 0.0 0.0
end

dft
xc lda
k 8
econv 1.e-6
dconv 1.e-4
canon
nuclear_corrfac none
end

excite
guess_thresh 1.e-4
solve_thresh 1.e-5
solve_seq_thresh 1.e-6
end

TDA
dft
compute_virtuals
guess big_fock_3
guess_excitations 1
iterating_excitations 1
excitations 1
iter_max 30
econv 1.e-3
dconv 2.e-2
guess_econv 1.e-2
guess_dconv 1.e-1
hard_dconv 1.e-4
hard_econv 1.e-6
no_otf
end

geometry
units angstrom
S .000000 .000000 .370268
O .000000 1.277617 -.370268
O .000000 -1.277617 -.370268
end

Am 05/08/15 um 23:06 schrieb lratcliff:

Thanks for the input Jakob, I'll have a look at your code to see if
there's anything obvious we're missing. Mostly though it already seems
to be working already for very small molecules, so I'm not sure how
much we can conclude anything from tests on cations. It seems to be
that the larger the system, the worse the problem, at least from my
small sampling of a dozen or so molecules.

One example where things are broken even with the asymmetric
orthogonalization is SO2. I'm just trying to calculate the LUMO, for
which BigDFT gives an energy of -0.18. Without the virtual state,
moldft converges in 10 iterations, whereas with the virtual state it's
still not converged after 80 iterations. I tried nearly doubling the
cell size, and it made no difference. Here's the input file I'm using:

dft
xc lda
protocol 1e-4 1e-6
maxiter 40
aobasis 6-31g
canon
nvalpha 1
nvbeta 1
end

geometry
units angstrom
S .000000 .000000 .370268
O .000000 1.277617 -.370268
O .000000 -1.277617 -.370268
end

—
Reply to this email directly or view it on GitHub
#127 (comment).

[lratcliff]

Sorry, I should have attached the output yesterday, I'm on vacation now without access to my desktop at ANL. But both the occupied and virtuals were messed up, essentially the energy was just fluctuating all over the place, and the residuals weren't going down at all.

That sounds really promising Jakob, thanks for doing the test. I had thought that doing occupied then virtuals would be a good route to go down, and this seems to confirm it. As you say, this would avoid problems with unwanted mixing. When I get back from vacation I'll try using your input for a few other molecules and if everything works then I suggest we adopt that approach in moldft.

[robertjharrison]

I agree that for the purpose of computing the virtuals, the right approach
is to converge the occupied to get the potential and then to converge the
virtuals.

However, adding a few virtuals should make convergence of the occupied
faster if we are doing diagnonalization (no effect if doing localized) and
should also make convergence more robust since it facilitates larger
changes in occupation number etc. Thus, I would still like to nail down
the problem.

It would also be good to do fractional occupation which is sort of the same
problem.

On Thu, Aug 6, 2015 at 2:20 PM, lratcliff notifications@github.com wrote:

Sorry, I should have attached the output yesterday, I'm on vacation now
without access to my desktop at ANL. But both the occupied and virtuals
were messed up, essentially the energy was just fluctuating all over the
place, and the residuals weren't going down at all.

That sounds really promising Jakob, thanks for doing the test. I had
thought that doing occupied then virtuals would be a good route to go down,
and this seems to confirm it. As you say, this would avoid problems with
unwanted mixing. When I get back from vacation I'll try using your input
for a few other molecules and if everything works then I suggest we adopt
that approach in moldft.

—
Reply to this email directly or view it on GitHub
#127 (comment)
.

Robert J. Harrison
tel: 865-272-9262

[lratcliff]

I've tested my latest fix (commit 2015578) on some larger molecules (coronene, phthalocyanine, various fullerenes), and for canonical orbitals everything seems to be working now.

Sometimes some of the higher energy (but still bound) virtual orbitals are missed by MADNESS, so I borrowed a trick that I previously used for virtual states in ONETEP: start with more virtual orbitals than needed and optimize for a few iterations before restarting with the number you actually want. This reorders some of the orbitals from the initial guess and the missing states are recovered.

For localized orbitals the original problem persists, i.e. as soon as you add any virtual states convergence is broken.