DHF convergence woes

[Peterson, Kirk]

Dear Dirac experts,

I generally hate emails like this, but I’m at my wits end. I’ve been trying to get a ECP-based, 2-component HF to converge all day today but to no avail. I’ve tried every trick I can think of but in every case the iterations will seem like they are converging but then in 1 iteration the energy will jump up by hundreds of milliHartrees or even a few Hartrees, then fight its way down but then eventually bounce back up. The system in question is UF6- . The single open-shell electron occupies the lowest 5f spinor on the central U. I’ve been using the converged closed-shell UF6 orbitals as the starting vectors. I was able to successfully get the spin free case to converge by first doing an average of configuration SCF of 1 electron in the 14 f-type spinors and then using that solution for the case of 1 electron in just 2 spinors (this solution nicely matches what I get from an ROHF calculation in Molpro). Using any of these orbitals (UF6 or spin-free UF6-) as initial vectors for the 2-component case doesn’t yield convergence. Level shifts, changing DIIS parameters, etc. doesn’t seem to effect the end result.

In principle, what would cause such huge bounces in the energy? The orbitals look reasonable and there do not seem to be any real near-degeneracy problems. Any tricks I can try that maybe aren’t so obvious to me? I’ve attached a representative output file.

thanks for any advice,

best regards,

-Kirk

[Trond SAUE]

Hi Kirk,
I can have a look at this one, but can you send me the basis set file
ECP60MDF ?
All the best,
Trond

--

Trond SAUE
Laboratoire de Chimie et Physique Quantiques
Universit� de Toulouse 3 (Paul Sabatier),
118 route de Narbonne, 31062 Toulouse (FRANCE)
t�l: 00 33 (0)561556031 FAX: (0)561556065

Email: trond...@irsamc.ups-tlse.fr
Web: http://dirac.ups-tlse.fr/saue
DIRAC: http://dirac.chem.vu.nl

[Peterson, Kirk]

Hi Trond,

certainly.  attached is the basis set file and the ECP file.

best regards,

-Kirk

On Nov 3, 2013, at 8:01 AM, Trond SAUE <trond...@irsamc.ups-tlse.fr> wrote:

> Hi Kirk,
> I can have a look at this one, but can you send me the basis set file ECP60MDF ?
>    All the best,
>       Trond
>
> --
>
>                          Trond SAUE
>         Laboratoire de Chimie et Physique Quantiques

>            Université de Toulouse 3 (Paul Sabatier),

>          118 route de Narbonne, 31062 Toulouse (FRANCE)

>            tél: 00 33 (0)561556031 FAX: (0)561556065

>
>              Email: trond...@irsamc.ups-tlse.fr
>               Web: http://dirac.ups-tlse.fr/saue
>                 DIRAC: http://dirac.chem.vu.nl
>

> --
> You received this message because you are subscribed to the Google Groups "dirac-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dirac-users...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.

[Trond SAUE]

Hi Kirk,
I tried to run you UF6 case with your inputs and got

=== RECP Error ===
Maximum symmetry of RECP integral is D2h

Do you have a special version ?

All the best,
Trond

--

Trond SAUE
Laboratoire de Chimie et Physique Quantiques

UMR 5626 CNRS --- Universit� Toulouse III-Paul Sabatier

118 route de Narbonne, 31062 Toulouse (FRANCE)

t�l: +33/(0)561556031 FAX: +33/(0)561556065

Email: trond...@irsamc.ups-tlse.fr
Web: http://dirac.ups-tlse.fr/saue

DIRAC: http://www.diracprogram.org
ESQC: http://www.esqc.org

[Peterson, Kirk]

Hi Trond,

there was a fix for this quite a while back (Dec. 2012) that involves commenting out one line, line 353, of /src/ecp/recp_lnk.F that performs this test.

best wishes and thanks again for looking into this, I really appreciate it.

-Kirk

On Nov 4, 2013, at 9:37 AM, Trond SAUE <trond...@irsamc.ups-tlse.fr> wrote:

> Hi Kirk,
> I tried to run you UF6 case with your inputs and got
>
> === RECP Error ===
> Maximum symmetry of RECP integral is D2h
>
> Do you have a special version ?
>
> All the best,
> Trond
>
> --
>
> Trond SAUE
> Laboratoire de Chimie et Physique Quantiques

> UMR 5626 CNRS --- Université Toulouse III-Paul Sabatier

> 118 route de Narbonne, 31062 Toulouse (FRANCE)

> tél: +33/(0)561556031 FAX: +33/(0)561556065

>
> Email: trond...@irsamc.ups-tlse.fr
> Web: http://dirac.ups-tlse.fr/saue
> DIRAC: http://www.diracprogram.org
> ESQC: http://www.esqc.org
>

[Trond SAUE]

HI Kirk,
thanks for the hint. I will try this...

All the best,
Trond

--

Trond SAUE
Laboratoire de Chimie et Physique Quantiques

Universit� de Toulouse 3 (Paul Sabatier),

118 route de Narbonne, 31062 Toulouse (FRANCE)

t�l: 00 33 (0)561556031 FAX: (0)561556065

Email: trond...@irsamc.ups-tlse.fr
Web: http://dirac.ups-tlse.fr/saue

DIRAC: http://dirac.chem.vu.nl

[Peterson, Kirk]

Thanks!

On Nov 4, 2013, at 2:30 PM, Trond SAUE <trond...@irsamc.ups-tlse.fr> wrote:

> HI Kirk,
> thanks for the hint. I will try this...
> All the best,
> Trond
>
> --
>
> Trond SAUE
> Laboratoire de Chimie et Physique Quantiques

> Université de Toulouse 3 (Paul Sabatier),

> 118 route de Narbonne, 31062 Toulouse (FRANCE)

> tél: 00 33 (0)561556031 FAX: (0)561556065

>
> Email: trond...@irsamc.ups-tlse.fr
> Web: http://dirac.ups-tlse.fr/saue
> DIRAC: http://dirac.chem.vu.nl
>

[Trond SAUE]

Hi Kirk,
I have worked a bit on UF6. If I run the neutral molecule it converges
like a charm, as expected.
Using projection analysis I find that the fluorines have charges -0.57
and uranium +3.45.
More imporant, uranium has electronic configuration
5f^1.23_6d^1.32_7s^0.06 in the molecule, so some 5f electrons are
already in use. Looking at the lower virtuals in ungerade symmetry the
first seven orbitals are indeed basically f orbitals.
The ungerade HOMO has energy -0.64206314 and the ungerade LUMO comes in
at -0.11488246.
However, if I try to calculate the anionic species
.CLOSED SHELL
44 42
.OPEN SHELL
1
1/0,14
starting from the coefficients of the neutral molecule, the SCF does not
converge. Yet, if I project the coefficients of the neutral species onto
the anionic (unconverged) one I find that the coefficients match rather
well, in particular that the open shell orbitals are still f. However,
now I observe that the highest inactive ungerade MO has energy
--0.40353915 and the lowest active MO comes in at -0.39634956. In other
words, they are almost degenerate and you will get artificial mixing and
convergence problems. I solved this by introducing a level shift of the
open shell
.OLEVEL
0.2
Now the SCF converges in 38 iterations. However, the convergence gets a
bit sluggish towards the end, which is typically an indication of
numerical noise, so I reduced the screening threshold on two-electron
integrals
**INTEGRALS
*TWOINT
.SCREEN
1.0D-16
(default is 1E-12), but this did not change the convergence pattern.
Then I decided to play with different values of the level shift

[saue@lpqsv3 kirk]$ ls *shift*
UF6_UF6.anion_shift_0.15 UF6_UF6.anion_shift_0.195
UF6_UF6.anion_shift_0.20 UF6_UF6.anion_shift_0.203
UF6_UF6.anion_shift_0.22
UF6_UF6.anion_shift_0.16 UF6_UF6.anion_shift_0.196
UF6_UF6.anion_shift_0.200 UF6_UF6.anion_shift_0.204
UF6_UF6.anion_shift_0.23
UF6_UF6.anion_shift_0.17 UF6_UF6.anion_shift_0.197
UF6_UF6.anion_shift_0.200.0 UF6_UF6.anion_shift_0.205
UF6_UF6.anion_shift_0.24
UF6_UF6.anion_shift_0.18 UF6_UF6.anion_shift_0.198
UF6_UF6.anion_shift_0.201 UF6_UF6.anion_shift_0.21
UF6_UF6.anion_shift_0.25
UF6_UF6.anion_shift_0.19 UF6_UF6.anion_shift_0.199
UF6_UF6.anion_shift_0.202 UF6_UF6.anion_shift_0.210
UF6_UF6.anion_shift_0.30

To my great surprise the SCF converges only for VERY specific values of
the level shift

[saue@lpqsv3 kirk]$ grep 'Convergence after' UF6_UF6.anion_shift_0.*
UF6_UF6.anion_shift_0.195:* Convergence after 50 iterations.
UF6_UF6.anion_shift_0.196:* Convergence after 36 iterations.
UF6_UF6.anion_shift_0.199:* Convergence after 47 iterations.
UF6_UF6.anion_shift_0.200:* Convergence after 38 iterations.

so I was apparently exceedingly lucky with my first choice ! I have no
clue why this is so.

Just for completeness, the electron configuration of uranium in the
anion is 5f^3.40_6d^1.17_7s^0.17.
Then, another surprise, the gross populations of each atom are

* Total gross contributions:
AFUXXX 30.4118 E - 30.4118 P - 0.0000
AFF1XX 9.2884 E - 9.2884 P - 0.0000
AFF2XX 9.2884 E - 9.2884 P - 0.0000
AFF3XX 9.4625 E - 9.4625 P - 0.0000
AFF4XX 9.4625 E - 9.4625 P - 0.0000
AFF5XX 9.5432 E - 9.5432 P - 0.0000
AFF6XX 9.5432 E - 9.5432 P - 0.0000

In other words, uranium has charge +1.58, but the fluorine atoms get
slightly different charges ranging from -0.29 to -0.54, although the
geometry is perfectly octahedric:

7

U z 0.0000000000 0.0000000000 0.0000000000
F z 0.0000000000 0.0000000000 2.0775205092
F z 0.0000000000 0.0000000000 -2.0775205092
F z 2.0775205092 0.0000000000 0.0000000000
F z -2.0775205092 0.0000000000 0.0000000000
F z 0.0000000000 2.0775205092 0.0000000000
F z 0.0000000000 -2.0775205092 0.0000000000

A possible explanation is that although DIRAC detects the full symmetry,
it treats the problem in lower symmetry in which the six fluorine atoms
are not symmetry equivalent
Full group is: O(h)
Represented as: D2h

All the best,
Trond




--

Trond SAUE
Laboratoire de Chimie et Physique Quantiques

Universit� de Toulouse 3 (Paul Sabatier),

118 route de Narbonne, 31062 Toulouse (FRANCE)

t�l: 00 33 (0)561556031 FAX: (0)561556065

Email:trond...@irsamc.ups-tlse.fr
Web:http://dirac.ups-tlse.fr/saue

DIRAC:http://dirac.chem.vu.nl

[Kenneth Dyall]

This sounds a bit familiar to me. The definition of the orbital eigenvalue is arbitrary for the diagonal blocks (the orbital energy doesn't depend on it). I've had several cases with my own code where I've had to give different level shifts to the open-shell orbitals and the virtuals, to separate them. Otherwise the open-shell spinor goes below a closed-shell spinor and they swap, with the Aufbau filling of the shells, hence the jumps. A Gershgorin disk analysis might be useful for cases like this.

Trond's comments about the symmetry made me think of another issue: have you done an analysis of the spinor occupations in the double group? Are you hitting symmetry-breaking due to the Jahn--Teller effect at the O_h symmetry?  If the fermion irrep for the open f shell is quadruply degenerate, then there will be a Jahn-Teller effect. You might consider looking at the virtuals of the neutral to see what the degeneracy is.

Ken.

            Université de Toulouse 3 (Paul Sabatier),

          118 route de Narbonne, 31062 Toulouse (FRANCE)

            tél: 00 33 (0)561556031 FAX: (0)561556065

              Email:trond...@irsamc.ups-tlse.fr
               Web:http://dirac.ups-tlse.fr/saue

                 DIRAC:http://dirac.chem.vu.nl

--
You received this message because you are subscribed to the Google Groups "dirac-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to dirac-users+unsubscribe@googlegroups.com.

[Peterson, Kirk]

Ken and Trond,

For the neutral, the 14 low-lying virtual spinors in E1u (essentially 5f on U) are:

-0.11488245423  ( 2)       -0.10109790673  ( 4)       -0.08640061786  ( 2)       -0.06890302668  ( 4)       -0.05865205543  ( 2)

So there is a quadruply-degenerate spinor, but it is not the lowest energy one (but quite close at this geometry) and shouldn't correspond to the ground state of the anion (which in the absence of SO the extra electron occupies a f_xyz orbital). But even so, for the anion I request an average-of-configuration DHF involving 1 electron in all 14 of these spinors, so would the JT effect really affect the HF convergence (especially at a frozen O_h geometry)?  I'm still surprised that while I could get the spinfree case to converge, using those converged coefficients as a guess for the 2-component case didn't help matters any. Though perhaps there is a clue in the spin-free results.  I proceeded to try a open-shell CCSD calculation (single reference) for the latter case but the CCSD equations would never come close to converging. Strangely I don't have any problems at that level with Molpro with the same basis and geometry but I'm not sure how to interpret that - the use of symmetry is quite different between these two calculations.

best regards,

-Kirk

To unsubscribe from this group and stop receiving emails from it, send an email to dirac-users...@googlegroups.com.

[Peterson, Kirk]

Dear Trond,

I did want to say I really appreciate all the work you put into helping me understand what is going on in this system. Would you mine sharing your input files for the projection analysis you carried out? This is certainly something I'd like to learn how to do. For UF6 I always assumed the U really wasn't in a +6 oxidation state, but it's interesting the charge is as small as ~3.5 though.

As an update, I've been running UF6 and UF6- in all-electron DHF calculations (using Ken's VDZ on U). Surprisingly UF6 was tough to converge - it took an initial level shift of 5.0 to obtain convergence. I'm running UF6- from these coefficients now and it looks like it might be smoothly converging with just a virtual orbital level shift (LSHIFT) of 0.1. Is there any reason why the ECP case would converge so differently from the all-electron case? I have another situation from a colleague (UO2+2 with 4 waters), where the all-electron converges easily (it's closed-shell) but the ECP calculation bounces like crazy (even spinfree).

best regards,

-Kirk

> Université de Toulouse 3 (Paul Sabatier),

> 118 route de Narbonne, 31062 Toulouse (FRANCE)

> tél: 00 33 (0)561556031 FAX: (0)561556065

>
> Email:trond...@irsamc.ups-tlse.fr
> Web:http://dirac.ups-tlse.fr/saue
> DIRAC:http://dirac.chem.vu.nl
>

[Peterson, Kirk]

Hi Trond,

one more question: how does applying an open-shell level shift (OSHIFT) affect the virtual orbital energies (which I would shift with LSHIFT)? In the below, if you shift the open-shell up by 0.2 doesn't that then put them above the lowest virtuals? Apologies for me being a bit dense on this.

best,

-Kirk

On Nov 6, 2013, at 1:21 PM, Trond SAUE <trond...@irsamc.ups-tlse.fr> wrote:

> Université de Toulouse 3 (Paul Sabatier),

> 118 route de Narbonne, 31062 Toulouse (FRANCE)

> tél: 00 33 (0)561556031 FAX: (0)561556065

>
> Email:trond...@irsamc.ups-tlse.fr
> Web:http://dirac.ups-tlse.fr/saue
> DIRAC:http://dirac.chem.vu.nl
>

[Jeff Hammond]

This seems like the sort of molecule that really needs to be treated
with the full symmetry, not just the highest non-degenerate subgroup.
Is there a reason why non-degenerate groups are implemented? I know
that it is considered nontrivial for correlated methods* but I though
it was relatively easy for SCF. I don't know about too many other
codes, but at least NWChem can do Td and Oh in the SCF modules. Of
course, NWChem doesn't do four-component methods, so perhaps there is
something I still need to learn to understand the implementation
details of symmetry in DHF-SCF solvers.

* I don't fully agree with this but since I've not yet implemented
non-degenerate, I'll try not to be too harsh a critic of the codes
that don't have it.

Jeff

> Université de Toulouse 3 (Paul Sabatier),

> 118 route de Narbonne, 31062 Toulouse (FRANCE)

> tél: 00 33 (0)561556031 FAX: (0)561556065

>
> Email:trond...@irsamc.ups-tlse.fr
> Web:http://dirac.ups-tlse.fr/saue
> DIRAC:http://dirac.chem.vu.nl
>
>

> --
> You received this message because you are subscribed to the Google Groups
> "dirac-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dirac-users...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.

--
Jeff Hammond
jeff.s...@gmail.com

[Young Choon Park]

Dear Kirk,

I read this topic yesterday (I couldn't frequently check this mailing list in these days).

I also tried your input and I found same results you had.

I found scf converged with ".FOCC" option
but the AOC calculation using DFCOEF obtained from ".FOCC" also failed.

In my experience (or feeling), some combinations of ECP and ECP basis fail or give erroneous results in 2c-molecluar calculation.

Evenif the ECP gave good numbers in spin-orbit calculation based on 1c-calculation (eg. SOCI in molpro or columbus),

2c-calculation couldn't give correct numbers (in energies or spin-orbit splittings).

And, compared to Christiensen ECP, I think I had experienced this matters more when I used Stuttgart ECP.

I havn't analyzed in detail at that time, I just thought some fittings or basis functions are not flexible for that target system.

(And, I remember SO parts in old Stuttgart ECP were fitted separately after AREP fitted, and this potential can be a problem in 2c-calculation)   

I haven't had time to check detail about your system yet.

If I can find any idea, I'll also post here.

with best regards,

Young Choon

2013/11/8 Jeff Hammond <jeff.s...@gmail.com>

[Peterson, Kirk]

Dear Young Choon,

thanks for taking a look at this.  I've since done some all-electron DHF and while I thought things were better, in the end they suffered from similar

problems as the ECP calculations.  Your point about the poor performance of some ECPs is well taken.  In particular I know that most of the older

Stuttgart-style PPs are not appropriate for 2-c calculations since the SO effects of the outer-core electrons are parameterized into the PP, i.e., they are

designed for valence-only SO calculations, and 2-c calculations will over-count SO effects with these PPs.  The newer ones based on MCDHF calculations 

definitely do not have this limitation.  I also always use uncontracted basis sets in 2-c SO calculations, since the original contracted sets generally aren't 

flexible enough since they are based on 1-c calculations.

Trond's used of OLEVEL seems to help with UF6-, although I'm still a little confused by how or why this works here.

best regards,

-Kirk

[Trond SAUE]

Hi Ken,

On 11/07/2013 05:18 AM, Kenneth Dyall wrote:

This sounds a bit familiar to me. The definition of the orbital eigenvalue is arbitrary for the diagonal blocks (the orbital energy doesn't depend on it). I've had several cases with my own code where I've had to give different level shifts to the open-shell orbitals and the virtuals, to separate them. Otherwise the open-shell spinor goes below a closed-shell spinor and they swap, with the Aufbau filling of the shells, hence the jumps. A Gershgorin disk analysis might be useful for cases like this.

If I understand you correctly your proposal is to calculate Gershgorin disks
http://en.wikipedia.org/wiki/Gershgorin_circle_theorem
in order to see if there is possible overlap between eigenvalues of open- and closed-shell orbitals ? Are there library routines for doing this analysis ?

Trond's comments about the symmetry made me think of another issue: have you done an analysis of the spinor occupations in the double group? Are you hitting symmetry-breaking due to the Jahn--Teller effect at the O_h symmetry?  If the fermion irrep for the open f shell is quadruply degenerate, then there will be a Jahn-Teller effect. You might consider looking at the virtuals of the neutral to see what the degeneracy is.

In octahedral symmetry, spinfree f orbitals split into three irreps: A2u (xyz), T2u and T1u whereas the alpha/beta spin functions span E_1/2g. Taking direct products we find
   A_2u  x E_1/2u = E_5/2u
   T_2u  x E_1/2u = E_5/2u + F_3/2u
   T_1u  x E_1/2u = E_1/2u + F_3/2u
so we see that the f orbitals split into five levels. This is indeed what I observe when looking at the virtuals of the neutral complex
    -0.11488246188  ( 2)      
    -0.10109792799  ( 4)     
    -0.08640061965  ( 2)      
    -0.06890302457  ( 4)      
    -0.05865205810  ( 2)
but also in the anion
  * Open shell #1, f = 0.0714
    -0.41640997755  ( 2)       -0.39406826730  ( 4)       -0.38459787949  ( 2)       -0.35595552793  ( 4)       -0.34783685745  ( 2)
I agree with Kirk that I can not quite see how a possible Jahn-Teller effect should affect the SCF convergence...I think the problem is that DIRAC detects Oh symmetry, but runs the problem in D2h and so possibly induces symmetry breaking.

All the best,
    Trond

Ken.

To unsubscribe from this group and stop receiving emails from it, send an email to dirac-users...@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

--
You received this message because you are subscribed to the Google Groups "dirac-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to dirac-users...@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

[Trond SAUE]

On 11/07/2013 09:35 PM, Jeff Hammond wrote:
> This seems like the sort of molecule that really needs to be treated
> with the full symmetry, not just the highest non-degenerate subgroup.
> Is there a reason why non-degenerate groups are implemented? I know
> that it is considered nontrivial for correlated methods* but I though
> it was relatively easy for SCF. I don't know about too many other
> codes, but at least NWChem can do Td and Oh in the SCF modules. Of
> course, NWChem doesn't do four-component methods, so perhaps there is
> something I still need to learn to understand the implementation
> details of symmetry in DHF-SCF solvers.
>
> * I don't fully agree with this but since I've not yet implemented
> non-degenerate, I'll try not to be too harsh a critic of the codes
> that don't have it.

Hi Jeff,
there are several issues here. Many quantum chemistry codes limit
symmetry to D2h and subgroups, what I call binary groups because the
effect on any symmetry operation on the coordinates is to possibly
change sign, but never mix coordinates. Symmetry can therefore be
handled efficiently and elegantly for these groups by logical bit
operations.
Things become much more complicated once you go to higher groups. There
are two basic strategies for the implementation of symmetry at the SCF
level: 1) You generate integrals in symmetry-adapted basis or 2) you
initially work in the original AO basis, but when constructing the Fock
matrix you first generate a skeleton matrix which is a matrix built from
a minimal set of integrals, the full matrix is afterwards generated by
applying symmetry operations on the skeleton matrix. There is a very
nice book chapter by Peter Taylor written for the European Summerschool
of Quantum Chemistry about this. Higher symmetry is easier to do with
skeleton matrices than at the integral level.
A complication in a relativistic framework is that spin-orbit
coupling mixes spin and spatial degrees of freedom, so you can not treat
spin and spatial symmetry separately. This is one reason why DIRAC does
not feature an open-shell Hartree-Fock code (only
average-of-configuration --- AOC) since the so-called coupling
coefficients appearing in the open-shell HF energy expression will now
depend on the point group symmetry. Note also that the convergence
behaviour of open-shell HF and AOC may be different.

All the best,
Trond

--

Trond SAUE
Laboratoire de Chimie et Physique Quantiques

Universit� de Toulouse 3 (Paul Sabatier),

118 route de Narbonne, 31062 Toulouse (FRANCE)

t�l: 00 33 (0)561556031 FAX: (0)561556065

[Jeff Hammond]

Thanks for the info. I cherish the ESQC lecture notes (in part
because the US Postal Service lost my first set and Bjorn Roos kindly
replenished them at no charge).

The NWChem SCF code is AO-based since that appears to be more
efficient for large systems when one is using screening to its full
advantage. I suppose DIRAC (and Dalton) use the "old way" (i.e.
symmetrized integrals).

I don't suppose anyone is working on the SCF solvers along such a
vector so as to realize the full set of point group support...

Best,

Jeff

--
Jeff Hammond
jeff.s...@gmail.com

[ueks...@gmail.com]

> I don't suppose anyone is working on the SCF solvers along such a
> vector so as to realize the full set of point group support...

I have a pretty fast code for rotating 3d polynomials that's looking for
a good application :-)

Regards,
Ulf

[Trond SAUE]

On 11/11/2013 09:37 PM, Jeff Hammond wrote:
> The NWChem SCF code is AO-based since that appears to be more
> efficient for large systems when one is using screening to its full
> advantage. I suppose DIRAC (and Dalton) use the "old way" (i.e.
> symmetrized integrals).

Hi Jeff,
in fact we use both, see
http://diracprogram.org/doc/release-12/manual/integrals.html#sofock

[Paul Bagus]

Dear Trond,

Actually there is another way to perform a transformation of the Fock matrix
into a higher symmetry, for example, going from D2h to Oh symmetry than
using a "skeleton" matrix. One can diagonalize one-electron matrices and
from the degeneracy and orthogonality patterns of the eigenvectors of these
one-electron matrices determine a set of orthogonal vectors in the low
symmetry basis set that belong to irreducible representations of the high
symmetry. One then does similarity transformations of the Fock matrix with
these sets of orbitals and obtains high symmetry blocked Fock matrices. Ulf
Wahlgren and I wrote a paper on this some time ago and it is attached. The
only catch is that one has to avoid accidental degeneracies and
othogonalities which are confused with different symmetries. Thus the
transformation to high symmetry is not completely black box and one has to
manually check if the symmetry reduction is correct.

I thought that Hans-Jorgen had done something similar. It could be you mean
something similar with the skeleton matrix idea.

I believe that the same kind of similarity transformation from low symmetry
Fock matrices to high symmetry blocked Fock matrices should work for Dirac
as well as for non-relativistic codes.

Best Regards, Paul

Université de Toulouse 3 (Paul Sabatier),

118 route de Narbonne, 31062 Toulouse (FRANCE)

tél: 00 33 (0)561556031 FAX: (0)561556065

Email: trond...@irsamc.ups-tlse.fr
Web: http://dirac.ups-tlse.fr/saue
DIRAC: http://dirac.chem.vu.nl

[Paul Bagus]

Dear Kirk,

Trond's experience with the need to fine tune the level shift started me
wondering about the possibility that the different orbitals require
different shifts to avoid the near degeneracies that lead to divergence. In
principle, it is possible to use a different level shift for each orbital
and this is what I do for the occupied orbitals, closed and open shell, in
my non-relativistic code.

One interesting test would be not to define the open shell as 1/0,14 but to
restrict the open shell orbitals to the lowest 5f's rather than to the full
set. The 5f's are spin-orbit split into 5/2 and 7/2 by about 0.7 eV. There
is then a further ligand field splitting into two sets of Oh* groups of G(7)
and G(8) in one group and G(6), G(7), and G(8) in the other. Would it be
possible for you to try on open occupation of 1/0,6 restricting the open
shell orbitals to the dominantly 5/2 group. It would be interesting if this
would eliminate or reduce the strong dependence of the convergence on the
precise values of the shifts.

Best Regards, Paul

-----Original Message-----
From: dirac...@googlegroups.com [mailto:dirac...@googlegroups.com] On
Behalf Of Trond SAUE
Sent: Wednesday, November 06, 2013 3:22 PM
To: dirac...@googlegroups.com
Subject: Re: DHF convergence woes

Université de Toulouse 3 (Paul Sabatier),

118 route de Narbonne, 31062 Toulouse (FRANCE)

tél: 00 33 (0)561556031 FAX: (0)561556065

Email:trond...@irsamc.ups-tlse.fr
Web:http://dirac.ups-tlse.fr/saue
DIRAC:http://dirac.chem.vu.nl

[Peterson, Kirk]

Dear Paul,

actually in my many attempts to get this to work, the idea below is one that I also tried, i.e., defining
the open shell as 1/0,6 and starting from converged DHF orbitals of the neutral. No luck after 75
iterations with the energy bouncing by as much as 200 mEh all the way to the end. I did not try any
level shifts with this attempt however. To my untrained eye, however, it seemed to me there is a good
separation between the highest closed-shell, 1st open shell, and 1st virtual orbital in the starting guess
from the neutral:

* Fermion symmetry E1u
* Closed shell, f = 1.0000
-26.36326247638 ( 4) -26.36324656361 ( 2) -10.64426609479 ( 2) -8.61707240295 ( 4) -1.76937338733 ( 2)
-1.65899732917 ( 4) -1.53424601624 ( 2) -1.26115997319 ( 4) -0.71306208952 ( 2) -0.71120008531 ( 4)
-0.69944091809 ( 2) -0.69738374670 ( 4) -0.67715236564 ( 2) -0.64206313529 ( 4)
* Virtual eigenvalues, f = 0.0000
-0.11488245423 ( 2) -0.10109790673 ( 4) -0.08640061786 ( 2) -0.06890302668 ( 4) -0.05865205543 ( 2) <--- last of f's
0.04729311910 ( 2) 0.04760535111 ( 4) 0.17882583239 ( 2) 0.18150231739 ( 4) 0.27339788162 ( 4)

So I admit to still being a bit confused by the need for so much level shifting.

In the end out of frustration, I've since moved to a FS-CCSD calculation where I could use the UF6 solution as the
00 sector.

best regards,

-Kirk

[Paul Bagus]

Dear Kirk,

I take the anthropomorphic view that there is some reason for the
convergence problems and that once we solve the problems, we will learn the
reason.

I have three suggestions for calculations with an open shell setting of
1/0,6:

1. Turn on the select occupied orbitals by overlap, .OVLSEL. This will track
the 5f and avoid flipping of the 5f and 7p. In my version of the code, I
print out the orbital that are selected at each iteration. I don't remember
if I added this or if it was in the base code.

2. Shift the virtual space up by a large amount, say ~ 2-3 hartrees and run
for 100 iterations or so. The logic of this is similar to that of the last
suggestion. This will force the virtual space up and out of the way and
dramatically slow mixing with the virtual space. The convergence will be
very slow but if it is monotonic, you can restart with a smaller shift. This
worked for me for core excited states of cluster models of UO6 clusters.

3. If 1 and 2 do not converge, try to combine them; both select by overlap
and shift the virtual space up.

Please let me know if these tricks work.

Good luck, Paul

[Trond SAUE]

On 11/12/2013 01:33 PM, Paul Bagus wrote:
> Shift the virtual space up by a large amount, say ~ 2-3 hartrees and run
> for 100 iterations or so.

Hi,
as a matter of principle any SCF calculation, if executed properly,
should converge i less thant 100 iterations.

All the best,
Trond

--

Trond SAUE
Laboratoire de Chimie et Physique Quantiques

UMR 5626 CNRS --- Universit� Toulouse III-Paul Sabatier

118 route de Narbonne, 31062 Toulouse (FRANCE)

t�l: +33/(0)561556031 FAX: +33/(0)561556065

Email: trond...@irsamc.ups-tlse.fr
Web: http://dirac.ups-tlse.fr/saue

DIRAC: http://www.diracprogram.org
ESQC: http://www.esqc.org

[Paul Bagus]

Dear Trond,

I agree. However, when one introduces a large virtual shift, which
artificially slows the mixing of the occupied and virtual spaces, the
convergence will be slow and it will become slower as the shift is made
larger.

I don't suggest using large shifts as an efficient way to obtain a converged
solution. If the device of a large shift leads to a stable and monotonic
path to convergence, then one can introduce an automatic algorithm to lower
the shift as long as the calculation continues to converge.

Regards, Paul

-----Original Message-----
From: dirac...@googlegroups.com [mailto:dirac...@googlegroups.com] On
Behalf Of Trond SAUE
Sent: Tuesday, November 12, 2013 6:35 AM
To: dirac...@googlegroups.com
Subject: Re: DHF convergence woes

On 11/12/2013 01:33 PM, Paul Bagus wrote:
> Shift the virtual space up by a large amount, say ~ 2-3 hartrees and
> run for 100 iterations or so.
Hi,
as a matter of principle any SCF calculation, if executed properly, should

converge in less than 100 iterations.

All the best,
Trond

--

Trond SAUE
Laboratoire de Chimie et Physique Quantiques

UMR 5626 CNRS --- Université Toulouse III-Paul Sabatier

118 route de Narbonne, 31062 Toulouse (FRANCE)

tél: +33/(0)561556031 FAX: +33/(0)561556065

Email: trond...@irsamc.ups-tlse.fr
Web: http://dirac.ups-tlse.fr/saue
DIRAC: http://www.diracprogram.org
ESQC: http://www.esqc.org

[Trond SAUE]

On 11/12/2013 01:50 PM, Paul Bagus wrote:
> Dear Trond,
>
> I agree. However, when one introduces a large virtual shift, which
> artificially slows the mixing of the occupied and virtual spaces, the
> convergence will be slow and it will become slower as the shift is made
> larger.
>
> I don't suggest using large shifts as an efficient way to obtain a converged
> solution. If the device of a large shift leads to a stable and monotonic
> path to convergence, then one can introduce an automatic algorithm to lower
> the shift as long as the calculation continues to converge.
>
> Regards, Paul

Hi Paul,
I feel that SCF convergence strategies based on a very large number of
iterations is not the way to go and a waste of computer time...

All the best,
Trond

--

Trond SAUE
Laboratoire de Chimie et Physique Quantiques

UMR 5626 CNRS --- Universit� Toulouse III-Paul Sabatier

118 route de Narbonne, 31062 Toulouse (FRANCE)

t�l: +33/(0)561556031 FAX: +33/(0)561556065

[Knut F{gri]

For what it may be worth, open-shell RHF systems are notoriously
difficult, and from experience require damping, shifting and sometimes
very many iterations. I agree that many iterations are not a desirable
feature of a calculation, but I have met cases where there was no
alternative.

As for the uranium calculations, I have also had a project om uranium that
refused to converge. The problem "resolved itself" when I became Dean of
the faculty of maths and sciences (and stopped doing active research), but
that is not a solution I would recommend in general, and it certainly did
not contribute much to uranium research.

Knut

[Peterson, Kirk]

Dear Paul,

I definitely agree with your first comment below. I just needed to get my project going, hence
the punt on my half. I will try your suggestions. Your strategy #2 below is exactly what I used to get the
neutral UF6 to converge in the first place. Otherwise it would never converge. I will play with your
option #1 today with the anion. I think I had tried this once, but perhaps not with neutral orbital starting guesses.

In a broader context, is there an intrinsic reason why a AOC DHF will be more difficult to converge
than a ROHF calculation, even in the spinfree case? This same molecule converges in 10 iterations
in Molpro with the same ECP/basis set combination (using no symmetry). My spinfree calculation in Dirac
using UF6 starting orbitals still took about 75 iterations (using a 1,0/14 open shell - trying 1,0/2 wouldn't converge
for me). I know it is in bad taste to compare codes, but I'd like to understand this. Molpro uses an atomic density guess,
maybe similar to the atomic start option in Dirac, but without having to do any additional calculations. Instead it stores
so-callled "guess basis sets" which are small, minimally contracted sets (using general contractions).

best regards,

-Kirk

[Paul Bagus]

Dear Kirk,

You ask a good question about why the SCF iterations converge for some
programs and not for others. I can give you some general considerations that
could point to an answer. My comments will indicate that I do not think that
it is a generic problem AOC against ROHF but rather how the Fock operators
are constructed in the different programs.

A reason why one does not obtain convergence is that there are either: (1)
incorrect rotations between the closed and open shell orbitals or (2)
incorrect rotations between the occupied and virtual spaces. It is the
second case that may be the problem for your case. Now different
formulations of the Fock operators generate different diagonal orbital
energies at each iteration and at convergence; i.e., different diagonal Fock
matrix elements for the Fock operator in the MO basis. Now, one can shift
the diagonal matrix elements as much as one wants and, if one converges, one
will converge to the same solution for any shift. This is why we are able to
use shifts. Thus, different treatments of the Fock operator may change the
diagonal matrix elements; see Roothaan's early papers on open shell SCF for
his treatment of the off-diagonal Lagrange multipliers to see how they
change diagonal matrix elements, especially for a one Fock operator
formalism. Sometimes these changes turn out to make calculations converge
faster and sometimes they will make the calculations diverge. If you are
interested, I can show you an example where one formulation leads to
convergence and another leads a converged solution to diverge after enough
iterations.

This does not answer your specific question but it suggests that looking at
the closed, open, and virtual orbital energies for the different programs
may give you the answer.

Best Regards, Paul

[Trond SAUE]

Hi Kirk,
DIRAC is not a commercial code and distributed as a service to the community. We have absolutely no qualms about being compared to other codes, rather I believe we can learn from it.
With the ECP setup you provided neutral UF6 smoothly in 26 iterations. If convergence is problematic with some other basis/setup, please send me inputs so I can test.
Some code, and I think MOLPRO is one of them, have a small level shift included from the start. The atomic start of DIRAC is efficient, but a bit clumsy to set up, something we may fix in the future.
For better convergence it should be coupled to the augmented Roothan-Hall method of the Aarhus group, see
The augmented Roothaan–Hall method for optimizing Hartree–Fock and Kohn–Sham density matrices,
S. Høst, J. Olsen, B. Jansík, L. Thøgersen, P. Jørgensen, and T. Helgaker,
J. Chem. Phys. 129, 124106 (2008)
also on my (very long) TODO list....

I agree with Paul that convergence problems are not due to ROHF vs AOC.

    All the best,
       Trond

-- 

			  Trond SAUE
         Laboratoire de Chimie et Physique Quantiques
    UMR 5626 CNRS --- Université Toulouse III-Paul Sabatier
          118 route de Narbonne, 31062 Toulouse (FRANCE)
            tél: +33/(0)561556031 FAX: +33/(0)561556065

              Email: trond...@irsamc.ups-tlse.fr
               Web: http://dirac.ups-tlse.fr/saue

              DIRAC: http://www.diracprogram.org
                  ESQC: http://www.esqc.org

[Luis Alvarez-Thon]

Dear all

I've been following this thread about convergence. I have not tried to do my own calculations on the problematic systems, but I think it could we interesting to apply an idea I have (perhaps it is not new!). The idea is to change the value of light a bit (C=130, C=140, etc) so that you can have a "level shift" due to spin-orbit. After convergence you can start from the coefficients with the normal C value. I have tried this method for open shell atoms and it works most of the time. Perhaps DIRAC can do this automatically.

I would like you to comment about this or perhaps try it on the UF6 and see what happens.

Best regards

Luis

[Peterson, Kirk]

Hi Trond,

you are exactly right, I must have been remembering some other situation.  UF6 indeed converged for me within 27 iterations without any level shifts or special guess vectors.  You are also correct that for ROHF calculations in Molpro there is a -0.3 Eh level shift applied to the occupied orbitals, which I take as being equivalent to a +0.3 open-shell shift in Dirac. Perhaps I will take the time to set-up an atomic start in Dirac and see if the behavior really is similar in that case.

best,

-Kirk

[Jeff Hammond]

What initial guess does DIRAC use? NWChem and GAMESS-UK use the
atomic density guess, which is a much better guess than the Huckel or
Hcore guess used by other codes.

MPQC has an ever better method, which is to solve the SCF with a tiny
basis like 3-21G that is usually trivial to converge, and then
projecting into the real basis. It takes about 10 lines of input file
options to do this in NWChem and I use it any time I'm having
convergence issues or using a large diffuse basis. When I run
something like aug-cc-pVQZ for anything with more than 4 atoms, I
telescope up from cc-pVDZ, projecting all the way.

In short, the only time I cannot converge an SCF with NWChem is
because the basis is so large that I get killed by linear dependencies
or the limits of double precision. I really only have serious
problems for situations like hexacene with aug-cc-pVQZ, which is
really nasty in a few respects.

Anyways, I think hand-wringing over level shifts is not the right
place to start. That should not be necessary unless you're SCF is
fundamentally unstable and oscillating between two solutions.
Otherwise, starting from a better guess should be the focus.

Best,

Jeff

--
Jeff Hammond
jeff.s...@gmail.com

[Radovan Bast]

hi Jeff,

corrected bare nucleus is the default. an option is to use

atomic density guess which is not black box but really easy to use.

and it is perhaps an advantage that it's not black box

because it gives more flexibility because you have more control.

it is very robust in most cases.

best greetings,

  radovan

[Trond SAUE]

HI Jeff,
an important aspect spelled out in the Aarhus paper I quoted in my
previous mail is that most SCF calculations use only gradient
information, which means that you may not even be in a local minimum,
but some saddle point, hence the interest of the augmented Roothan-Hall
approach.

All the best,
Trond

--

Trond SAUE
Laboratoire de Chimie et Physique Quantiques

Universit� de Toulouse 3 (Paul Sabatier),

118 route de Narbonne, 31062 Toulouse (FRANCE)

t�l: 00 33 (0)561556031 FAX: (0)561556065

Email: trond...@irsamc.ups-tlse.fr
Web: http://dirac.ups-tlse.fr/saue

DIRAC: http://dirac.chem.vu.nl

[Peterson, Kirk]

Dear Paul,

below are the results of the tests you recommended this morning. As you'll note, I did have success but it took both open-shell
and virtual orbital shifts. With those options, however, I could easily converge both 1/0,6 and 1/0,14 open-shell cases.

best regards,

-Kirk


1) Using just .OVLSEL produced a diverging solution quite quickly (UF6 start vectors)

2) Using a big virtual orbital level shift (.LSHIFT of 3.0) seem to start things in the right direction but by around iteration 30 it
started to diverge by bouncing up by about 2 Eh and then never coming back down within 75 iterations. (UF6 start vectors)

3) Combining 1 and 2 produced results similar to #1

4) Upon looking at the orbital energies from Molpro for both UF6 and UF6-, I realized that the closed-shell orbital energies in the anion
were significantly higher than in the neutral. While in the former they were lower than the virtuals by more than 0.2 Eh, in the anion they were essentially
isoenergetic with the open-shell orbital. This of course matches what Trond had also observed. So that obviously argues for a open-shell shift, not a virtual orbital
shift. Unfortunately use of just .OLEVEL = 0.3 started out good but then eventually diverged. As Trond noted, convergence is very sensitive
to the choice of shift in this case.

5) Then came success. I combined a .OLEVEL value of 0.3 (which is the same that Molpro uses in their ROHF program) with a .LSHIFT of 0.3.
This converged both the 1/0,6 and 1/0,14 cases in 22 and 35 iterations, respectively (both starting from the same UF6 converged vectors).

6) Feeling perhaps overconfident, I then tried the same input as #5 but used the default bare nucleus Hamiltonian guess. To my surprise the 1/0,14 case
converged in just 24 iterations.

I don't know if this is the end of this particular story, thanks to all for a very enlightening thread. I've got some additional reading to do now :^)

best regards,

-Kirk

On Nov 12, 2013, at 4:33 AM, Paul Bagus <pba...@austin.rr.com> wrote:

[Kenneth Dyall]

Looks like I returned to this thread a bit too late. The solution that Kirk found to be successful was exactly what I found earlier for other open-shell cases with my own code (DREAMS). You have to separate both the occupied orbitals and the virtual orbitals from the open-shell orbitals, which you can do by shifting the open-shell orbitals up and the virtuals up as well. Sometimes I've used 0.5 for the oshift and 1.0 for the vshift, which creates a gap of 0.5 between closed and open, and a gap of 0.5 between open and virtual orbitals. Sorry I didn't make this clear in my earlier post, as it could have saved a lot of trouble!

Regarding the Gershgorin disks. They are trivial to calculate from the Fock matrix (as you can see from any of the articles e.g. on Wikipedia). The idea would be that you calculate the disks for all of the diagonal elements. Then you would apply a shift to the open shell orbitals to ensure there was no overlap with the closed shell orbitals, and a shift to the virtuals to ensure there was no overlap with the open shell (or the closed shell) orbitals. This then ensures that the closed, open, and virtual orbitals can't swap. It also (incidentally) would guard against pathologies in which numerical noise created very large off-diagonal matrix elements such that the negative-energy states had overlapping disks with the positive-energy states, for example.

All the best,
Ken.

[Kenneth Dyall]

Just to follow up on my last comment about Gershgorin disks. I am fairly confident that implementing such as scheme would do away with many convergence problems. You could evaluate the disks for each iteration (it's N^2, so not expensive) and adjust the shifts to separate the groups of orbitals that are in the same space (i.e. are independent of rotations within the space) so that the overlap is removed by a very small amount. This would keep the convergence speed up while removing the pathologies. You might even use a negative shift to accelerate convergence, if my idea is correct.

Ken.

[Kenneth Dyall]

Dear Paul and Trond,

Determining the irreps of a higher group from those of a lower group can be aided by considering atoms at the symmetry center or on the axis if there is no center, as the classification of these into the high-symmetry irreps are well known. Combined with the diagonalization and blocking that you mention, this can help to resolve the classifications.

Ken.

[Kenneth Dyall]

Regarding the possibility of a Jahn-Teller effect, as Paul points out the f_5/2 splits into the doubly-degenerate Gamma_7_-  irrep and the quadruply-degenerate Gamma_8^- irrep.  If the latter is lower, then you might see a Jahn-Teller effect.

Ken.