moldft: virtual orbitals are broken

[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 https://github.com/m-a-d-n-e-s-s/madness/issues/127#issuecomment-123873991 .

[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 https://github.com/m-a-d-n-e-s-s/madness/issues/127#issuecomment-123873991 .

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 https://github.com/m-a-d-n-e-s-s/madness/issues/127#issuecomment-125239139 .

[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 https://github.com/m-a-d-n-e-s-s/madness/issues/127#issuecomment-125239139 .

� 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 https://github.com/m-a-d-n-e-s-s/madness/issues/127#issuecomment-125253890 .

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 https://github.com/m-a-d-n-e-s-s/madness/issues/127#issuecomment-128148482 .

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 https://github.com/m-a-d-n-e-s-s/madness/issues/127#issuecomment-128148482.

[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 https://github.com/m-a-d-n-e-s-s/madness/issues/127#issuecomment-128464997 .

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.