"SCF energy in HF" and "SCF energy in RELCC" don't match each other.

[Tomohiro Ichibha]

Dear all,

Q(a). 

I have run HF -> RELCC and found that the SCF energies 

written in HF and RELCC parts are different: 

    SCF energy in HF        : -3383.6651884103194

    SCF energy in RELCC: -3383.772933531707167

Is it fine? If so what is the reason of the difference?

(BTW, what is the energy unit? Hartree?)

I have attached the output file as chg0.out.  

Q(b). 

I have run HF -> RELCC for +1 charged state for the same system

and found that the energy difference becomes larger:

    SCF energy in HF        : -3383.5576026759754

    SCF energy in RELCC: -3383.774443897787933

The "SCF energy in RELCC" is actually very close to that in Q(a),

so RELCC did not recognize .CHARGE in *CHARGE?

I have attached the output file as chg1.out.

Best Regards,

Tomohiro

[Lorenzo Lodi]

I have also noticed your first point, and I also would like an answer to it.

I.e., why the SCF energies printed right after the SCF calculation are not the same as after the correlated, eg CCSD(T), step.
In fact, I think I can partially answer this question. The key is in the following bit in your outputm during the "Set-up for Coupled Cluster calculations" step:

QUOTE

 Checking the orbital energies, the program computes the diagonal elements of the

 reconstructed Fock matrix. Differences with the reference orbital energies

 are given if above a treshold or if iprnt > 1

 Spinor   Abelian Rep.         Energy   Recalc. Energy

  O    1    1    1g       -6.4132446820   -6.6496300111

........

UNQUOTE

The Fock matrix is at this step "reconstructed" (whatever that means) and, orbital energies change and therefore presumably also the total SCF energy changes.

The manual mentions an option .USEOE and it says 

QUOTE

Ignore recomputed Fock matrix and use orbital energies supplied by the SCF program. This option is sometimes useful for degenerate open shell cases in which case the perturbation theory for the unrestricted formalism is not invariant for rotations among degenerate orbitals. It should only change the outcome of the [T], (T) and -T energy corrections.

UNQUOTE

I think this option would do the trick and keep the SCF energy the same in the CCSD(T) step. However, I ignore the technical details and I don't know if this is good or bad.

Also, I'm not sure how to use this option in the input file (I tried but failed).

I'm trying to compute net relativistic contributions as differences of relativitistic and non-relativistic energies, and the issue of which SCF energy I should use becomes relevant.

Any advice?

Lorenzo Lodi

On Thursday, 16 January 2020 16:21:32 UTC+1, Tomohiro Ichibha wrote:

Dear all,

Q(a). 

I have run HF -> RELCC and found that the SCF energies 

written in HF and RELCC parts are different: 

    SCF energy in HF        : -3383.6651884103194

    SCF energy in RELCC: -3383.772933531707167

Is it fine? If so what is the reason of the difference?

(BTW, what is the energy unit? Hartree?)

I have attached the output file as chg0.out.  

[...]

[Visscher, L.]

Dear Lorenzo and Tomohiro,

It is tricky to define the uncorrelated reference energy in CC calculations. 

What is printed is actually the expectation value of the determinant that is specified via the input keyword .NELEC or (like in the Yttrium example of Tomohiro) the keyword .NELEC_OPEN. For closed shell molecules this will be identical to the HF energy and this is why its is called SCF energy. In such cases it just serves to check that no errors were made in defining the integral transformation, and it should be identical to the energy printed by the Hartree-Fock module. An exception is when approximations are made in the integral transformation such as in the XC-MMF Hamiltonian in which two-electron relativistic corrections on valence orbitals are neglected (picture change error). In that case it is better to use the correct orbital energies and SCF energy as taken from the HF code and for this reason the recomputed Fock matrix is ignored (we know it will be slightly different but also know the correct value and can use that instead). This does remove a check, but as long as you work with closed shell miolecules and do not specify .NELEC explicitly (the program will calculate this for you) you should be fine.

For open shell systems things are different and no single advice can be given. You need to understand how the energy average is defined in the open shell SCF procedure and what is the best single determinant approximation to the electronic state that you are interested in. In the January email that was quoted here the open shell concerned a total average of the 5s,4d)3 configuration of the Yttrium atom. This average takes all possible occupation of these orbitals, so averages over doublet as well as many quartet states (using non-relativistc nomenclature just for convenience, they are miixing due to spin-orbit coupling of course). This energy will be quite a bit higher than the expectation value of a properly chosen determinant like the 5s)2 4d_{3/2}^1 that was selected by Tomohiro and this is exactly what you see. Tricking the program in using the HF computed energy is not a good idea here because this represents something else.

So take home messages:

1) SCF energy printed by HF program represents an average of configuration energy.

2) SCF energy printed by RELCCSD represents an expectation value of a single determinant

1) and 2) are only strictly equal for a closed shell molecule. 

hope this helps, best regards,

Luuk

--
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.
To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/80a0aed7-fb61-4dd1-8510-f5b97a069629o%40googlegroups.com.

[Lorenzo Lodi]

Dear Luuk,

Thanks a lot for your explanation, it clarifies things a lot. 

In my case I'm doing calculations on closed-shell molecules and I noticed differences between the two SCF energies (SCF printed by the HF program and as reported by RELCCSD) when then Gaunt interaction is included. I attach an output for the hydrogen fluoride diatomic molecule.

The SCF value after the SCF run is:

Total energy                             :     -99.95486257546007

while the value reported at the of CCSD(T) is

SCF energy :                                -99.955591019913243 

I guess this is due to the following fragment from your previous message:

An exception is when approximations are made in the integral transformation such as in the XC-MMF Hamiltonian in which two-electron relativistic corrections on valence orbitals are neglected (picture change error).

Am I correct?

Dirac also says, even before the SCF run, 

 * WARNING: Gaunt only available for Hartree-Fock,  and will be ignored in 2-electron integral transformation.

 * INFO: Gaunt only implemented for AO Fock builder ---> switched to AOFOCK

I want to compute the net effect of the Gaunt term (if I understand correctly this makes sense only at the SCF level). 

Which SCF value should I take as the correct Gaunt-correct one? It makes a big difference!

Thanks,

Lorenzo

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

[Visscher, L.]

Dear Lorenzo,

With your input you get the most accurate energies by using the SCF energy as printed by the HF program and manually adding the correlation energy printed by RELCCSD. The difference between the two indeed comes from the neglect of the Gaunt interaction in the index transform, you essentially calculating a DC only energy at the CC level and a DCG energy at the HF level. Gaunt should be quite well described at the mean field level of theory,  which is one reason why we have so far not implemented this in the index transformation (another reason is that this implementation is quite challenging to get efficient, as you are essentially transforming 3-dimensional  current density matrices rather than the 1-dimensional scalar density matrix).

If you really want to make sure you can try the old MOLFDIR code in which I did implement the Gaunt transformation. The program is not so efficient as DIRAC, but for a molecule like HF calculations are easily feasible. This code is available at 

https://github.com/MOLFDIR/MOLFDIR

We do plan to introduce the Gaunt transformation functionaliity in DIRAC as well of course, but this will probably be done in the new massively parallel CC code that we are working on.

best regards,

Luuk

To unsubscribe from this group and stop receiving emails from it, send an email to dirac-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/8cdb7575-3966-479d-a672-c7744ab06f6eo%40googlegroups.com.
<output_HF.txt>

[Lorenzo Lodi]

Dear Luuk,

thanks, that helps a lot! Thanks for the link to MOLFDIR; I'll hope to give it a try soon.

I hope you get round to implementing the Gaunt correction transformation in Dirac; in fact, it'd be great also to see a full Breit implementation. 

I believe that the program Bagel by Shiozaki fully implements it, but I've been unable to run that software on my system so far.

I've noticed something else which I don't understand when the .Gaunt option is included in my hydrogen fluoride tests.

Let me briefly explain. 

When I do calculations with .GAUNT using the fully uncontracted aug-cc-pCVTZ basis set (all electrons correlated), then the energies of correlated methods (MP2, CCSD, CCSD(T) ) is almost identical (within 1e-6 Eh) to the energies obtained without the .Gaunt keyword. And in light of what we've said, this is okay.

On the other hand, when using the (uncontracted) aug-cc-pVTZ basis set and keeping the 1s core electrons frozen, the energies of MP2, CCSD, CCSD(T) when .Gaunt is included are very different from the corresponding energies obtained without the .Gaunt keyword (differences of about 1.1e-2 Eh).

Do you have any explanation as to what is going on?

Here are some numbers if you care to look. I can include output files if necessary. The values are for the same geometry of the output I attached in my previous message, ie r=1.2 a0.

I called SCF-A the energies printed right after the SCF program and SCF-B the energies printed out after running the CCSD(T) method. 

uncontracted aug-cc-pCVTZ basis set

ALL ELECTRONS CORRELATED (ALL VIRTUAL ORBITALS INCLUDED)

Gaunt No Gaunt

SCF-A:    -99.957444  -99.969591

SCF-B:    -99.969589  -99.969591

CCSD(T): -100.309989 -100.309988 <=== same energies

aug-cc-pVTZ basis set

FROZEN CORE (ALL VIRTUAL ORBITALS INCLUDED)

Gaunt No Gaunt

SCF-A:    -99.954863  -99.966987

SCF-B:    -99.955591  -99.966987

CCSD(T): -100.230268 -100.241693 <=== different energies

Best regards,

Lorenzo

Luuk

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/8cdb7575-3966-479d-a672-c7744ab06f6eo%40googlegroups.com.
<output_HF.txt>

[Lorenzo Lodi]

Incidentally, the correction to the shape of HF molecule potential energy curve due to .GAUNT comes out very large,  in fact about as large as the overall non-Gaunt relativistic correction (ie, the difference between relativistic and non-relativistic energies). I wonder if this can be right as, especially for such a light molecule, I'd expect the Gaunt term to be small.

[Visscher, L.]

Dear Lorenzo,

This indeed looks suspicious. How did you get the Gaunt difference exactly ? Most fullproof is probably to take the difference between the SCF energies as printed by the HF module.

best regards,

Luuk

On Aug 20, 2020, at 10:52 AM, Lorenzo Lodi <lodi...@gmail.com> wrote:

Incidentally, the correction to the shape of HF molecule potential energy curve due to .GAUNT comes out very large,  in fact about as large as the overall non-Gaunt relativistic correction (ie, the difference between relativistic and non-relativistic energies). I wonder if this can be right as, especially for such a light molecule, I'd expect the Gaunt term to be small.

--
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.
To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/b993a840-1398-4081-830c-dd84281ad392n%40googlegroups.com.

[Kenneth Dyall]

The 1-e relativistic corrections scale with an extra factor of Z compared to the 2-e corrections, so proportionately the 2-e corrections are larger in the light elements than in the heavy elements. However, the Gaunt term doesn't have any direct contribution to a closed shell molecule, only exchange, so I wouldn't expect the Gaunt term to be large. It does contribute a scalar term, but mostly it contributes to the spin-other-orbit and spin-spin terms. 

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/47E9BB98-CE65-4288-994A-BD8D17956FE1%40vu.nl.

[Lorenzo Lodi]

Thanks for the answers. 

Here are some number (output of the SCF method, computed with the uncontracted aug-cc-pCVDZ basis set).

First column:  bond length in bohrs;

Second column:  energy with .Gaunt and .Dossss

Third column:  energy with .Lvcorr and .Dossss

Fourth column: energy with .LEVY-LEBLOND 

r/a0  GAUNT          LVCORR        LEVY-LEBLOND

1.2 -99.926621 -99.938723 -99.847181

1.3 -100.012855 -100.024909 -99.933459

1.4 -100.065988 -100.078004 -99.986612

1.5 -100.096406 -100.108392 -100.017034

1.6 -100.111211 -100.123173 -100.031833

1.7 -100.115311 -100.127254 -100.035920

1.8 -100.112129 -100.124057 -100.032721

1.9 -100.104060 -100.115975 -100.024633

2.0 -100.092778 -100.104682 -100.013330

2.1 -100.079446 -100.091342 -99.999979

2.2 -100.064872 -100.076762 -99.985387

2.3 -100.049619 -100.061504 -99.970116

2.4 -100.034078 -100.045959 -99.954559

2.5 -100.018525 -100.030402 -99.938991

2.6 -100.003153 -100.015026 -99.923606

2.7 -99.988097 -99.999968 -99.908538

2.8 -99.973451 -99.985320 -99.893881

2.9 -99.959277 -99.971145 -99.879698

3.0 -99.945613 -99.957480 -99.866026

What I called "non-Gaunt relativistic correction" is the difference LEVY-LEBLOND - LVCORR.

What I called "Gaunt correction" is the difference LVCORR - GAUNT

The former correction has an average value of -0.0914 Eh, and the latter of 0.0119 Eh, so in absolute terms the Gaunt contribution is about 8 times smaller, but I'm concerned with the position dependence.

Shifting the correction curves so that they are 0 at r=1.7 a0 I get the following values (in mEh):

r/a0       GAUNT    NON-GAUNT

1.20 0.159 -0.208

1.30 0.111 -0.116

1.40 0.073 -0.058

1.50 0.043 -0.024

1.60 0.019 -0.006

1.70 0.000 0.000

1.80 -0.015 -0.002

1.90 -0.028 -0.008

2.00 -0.038 -0.018

2.10 -0.046 -0.030

2.20 -0.053 -0.042

2.30 -0.058 -0.054

2.40 -0.062 -0.065

2.50 -0.066 -0.077

2.60 -0.069 -0.087

2.70 -0.071 -0.096

2.80 -0.073 -0.105

2.90 -0.075 -0.113

3.00 -0.076 -0.120

And the two correction curves have then quite comparable magnitude.

Using larger basis sets doesn't change the conclusion that the Gaunt term changes the shape of the curve quite significantly. I wasn't expecting it.

Lorenzo

[Peterson, Kirk]

Dear Lorenzo,

 

I took the liberty of fitting your data and doing a standard Dunham analysis (using just the 9 points between 1.4 and 2.2 bohr):

 

                       Re (Å)      we (cm-1)     wexe (cm-1)

NR-HF         0.89939       4462.42         90.63

DC-DHF       0.89938      4460.39          90.71

DC+Gaunt   0.89950      4459.57          90.91

 

So overall, I would say that as expected the Gaunt contribution is a bit larger than I would have expected, but still somewhat smaller than the usual scalar contribution for the harmonic frequency (we).

 

best regards,  -Kirk

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/4a61cccf-ddb9-492e-ba19-834695995f94n%40googlegroups.com.

[Visscher, L.]

Dear Kirk, Lorenzo,

This fit nicely condenses the information, This does not look so strange, although maybe a bit larger than one would expect on forehand. In my experience Gaunt adds a bit of extra repulsion between the atoms, typically making bond lengths sligtly longer and frequencies slighty smaller. This is in line with what comes out of the fit. 

One thing to check is activating both .DOSSSS and .LVCORR. This should not be done as it woiuld double count the 2-center SS repulsion (.LVCORR does this approximately, .DOSSSS wiill do this explicitly). I am not sure whether we check on this combination and deactivate the correction. Can you run wiith only .DOSSSS and see if this makes a difference ?

best regards,

Luuk

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/2DA0108C-F5F0-48E1-883A-067ACA5869E1%40wsu.edu.

[Lorenzo Lodi]

Dear Luuk, 

concerning .DOSSSS: the energies I posted under the column 'LVCORR' were obtained using the .DOSSSS keyword (but not specifying LVCORR); sorry for now having been precise. The Gaunt energies used both .GAUNT and .DOSSSS keywords; in any case, using the .LVCORR keyword gives a curve which is extremely close to the one obtained with .DOSSSS.

I think the situation is clear now and that I haven't made any silly mistake. 

For those interested, here are some more data with relative corrections to the J=0 vibrational energy levels.

All values obtained with auc-cc-pCVDZ-uncontracted and SCF.

Raw energies:

  .GAUNT/.DOSSSS .DOSSSS .LVCORR .SPINFREE .LEVY-LEBLOND
1.20 -99.926621 -99.938723 -99.938732 -99.938726 -99.847181
1.30 -100.012855 -100.024909 -100.024918 -100.024912 -99.933459
1.40 -100.065988 -100.078004 -100.078012 -100.078006 -99.986612
1.50 -100.096406 -100.108392 -100.108401 -100.108395 -100.017034
1.60 -100.111211 -100.123173 -100.123182 -100.123176 -100.031833
1.70 -100.115311 -100.127254 -100.127263 -100.127256 -100.035920
1.80 -100.112129 -100.124057 -100.124065 -100.124059 -100.032721
1.90 -100.104060 -100.115975 -100.115984 -100.115977 -100.024633
2.00 -100.092778 -100.104682 -100.104691 -100.104684 -100.013330
2.10 -100.079446 -100.091342 -100.091351 -100.091344 -99.999979
2.20 -100.064872 -100.076762 -100.076771 -100.076763 -99.985387
2.30 -100.049619 -100.061504 -100.061512 -100.061505 -99.970116
2.40 -100.034078 -100.045959 -100.045967 -100.045959 -99.954559
2.50 -100.018525 -100.030402 -100.030410 -100.030402 -99.938991
2.60 -100.003153 -100.015026 -100.015035 -100.015026 -99.923606
2.70 -99.988097 -99.999968 -99.999977 -99.999968 -99.908538
2.80 -99.973451 -99.985320 -99.985329 -99.985319 -99.893881
2.90 -99.959277 -99.971145 -99.971153 -99.971143 -99.879698
3.00 -99.945613 -99.957480 -99.957489 -99.957478 -99.866026

J=0 energy levels (nuclear masses were used)

v .GAUNT/.DOSSSS .DOSSSS .LVCORR .SPINFREE .LEVY-LEBLOND
0 2211.90 2212.35 2212.35 2212.37 2213.43
1 6495.65 6496.96 6496.97 6497.01 6500.14
2 10611.12 10613.24 10613.24 10613.32 10618.43
3 14567.13 14570.01 14570.01 14570.13 14577.13
4 18371.66 18375.25 18375.26 18375.41 18384.21
5 22031.89 22036.15 22036.15 22036.35 22046.86
6 25554.24 25559.13 25559.13 25559.38 25571.50
7 28944.52 28950.00 28950.00 28950.30 28963.94
8 32208.01 32214.05 32214.05 32214.40 32229.48
9 35349.61 35356.16 35356.16 35356.58 35373.01

Net contributions (differences between consecutive columns) to J=0 energy levels

v .GAUNT/.DOSSSS .DOSSSS .LVCORR .SPINFREE

0 -0.45 0.00 -0.01 -1.06
1 -1.31 0.00 -0.05 -3.13
2 -2.12 0.00 -0.08 -5.11
3 -2.88 0.00 -0.12 -7.00
4 -3.59 0.00 -0.16 -8.80
5 -4.26 0.00 -0.20 -10.51
6 -4.89 0.00 -0.25 -12.12
7 -5.48 0.00 -0.30 -13.64
8 -6.03 0.00 -0.35 -15.08
9 -6.55 0.00 -0.41 -16.43

So my conclusions are:

1) The .DOSSSS keywords affects the energy levels by less than 1e-2 cm-1 and has therefore negligible effect  (its effect is of the order of 2e-3 cm-1, although it cannot be seen from the data above).

2) Spin-dependent terms have a small effect on energy levels.

3) The Gaunt term brings down the energy levels, and its effect is about 0.4 times the non-Gaunt relativistic effect. In other words, Gaunt increases the magnitude of the relativistic shifts to energy levels by 40%. In agreement with Kirk's calculations, its effect on the v=0 to v=1 transition is -0.86 cm-1, and grows approximately linearly for higher v.

Thank you all for the support!

Kind regards,

Lorenzo

You received this message because you are subscribed to a topic in the Google Groups "dirac-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dirac-users/SL_McWyuIwU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dirac-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/4F10728C-85D6-49E0-9ED0-FE35E0DB07BD%40vu.nl.

[Paul Bagus]

Dear Kirk.

 

There is an uncertainty in a Dunham analysis depending on the choice of points and the degree of the polynomial fit.

 

You didn’t mention the degree of fit that you used but if you get wexe, you must have at least a cubic fit (maybe quartic(

 

Would it be easy for you to change the degree of polynomial ft to see whether and by how much the choice of polynomial fit affects the comparison you made.

 

Regards, Paul

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/2DA0108C-F5F0-48E1-883A-067ACA5869E1%40wsu.edu.

[Peterson, Kirk]

Dear Paul,

 

I actually used a 7th degree polynomial in bond displacement coordinates to 9 points, so it was close to a polynomial interpolation. Lower order polynomials of course decreased the accuracy of the fit, but it seemed the general trends were unchanged since the fits were of very similar quality for each curve.  I can certainly post a comparison if this is of interest.

 

regards,  -Kirk

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/BYAPR14MB21679150A6DDDF1C7F183D32BD5B0%40BYAPR14MB2167.namprd14.prod.outlook.com.

[Peterson, Kirk]

Dear Paul,

 

see attached.

 

regards,  -Kirk

 

From: <dirac...@googlegroups.com> on behalf of Paul Bagus <ba...@paulbagus.com>

Reply-To: "dirac...@googlegroups.com" <dirac...@googlegroups.com>
Date: Friday, August 21, 2020 at 7:01 AM
To: "dirac...@googlegroups.com" <dirac...@googlegroups.com>

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/BYAPR14MB21679150A6DDDF1C7F183D32BD5B0%40BYAPR14MB2167.namprd14.prod.outlook.com.

[Paul Bagus]

Dear Kirk,

 

A caution about polynomial fits. If you have the same number of coefficients in the polynomial as data points, as you write, you get a perfect fit. The polynomial exactly passes through all the data points. However the fit is likely to be physically meaningless. Hartree and Hartree wrote a book on numerical analysis with a beautiful example of how such a fit can go wrong. The polynomial, I think for a potential curve, had several oscillations to give a perfect fit but did not have a nice r(e) or w(e).

 

I have added printing the maximum and average error of the fit to the code that I use for a Dunham analysis. Large errors can come because the range of distances I use goes beyond the region where a polynomial is valid or there is an error in the data for the points. I try to avoid polynomials higher than quartic. Often, if the fit has errors, I reduce the range of distances in the points that I use.

 

Best Regards, Paul

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/1E893BB2-CD32-411A-ABDB-410A3CA901EC%40wsu.edu.

[Peterson, Kirk]

Dear Paul,

 

this is why I look very carefully at the coefficients of the resulting polynomial. It is straightforward for these types of fits to determine if things are well behaved or not. In general the coefficients should alternate in sign and the magnitudes should be dropping.  I also only use this analysis for a very specific range of points around the minimum - generally no more than -0.4 bohr to +0.7 bohr and I use a specific distribution of points.  But you're of course exactly right that polynomial fitting has the danger of giving one oscillations and odd behavior. Fortunately most ground state potential curves are very well behaved, at least in the near-equilibrium region and inspection of the coefficients is very reliable to tell if things are ok.  But one can certainly see in what I posted for HF that even small changes in the fit can produce differences of more than a cm-1 in the harmonic frequency.

 

My code also prints average error and errors at each data point, so this can be helpful (as well as the range of data points in units of omega_e). The problem with using polynomials of only order quartic is that the fits are generally somewhat poor (generally much less than microHartree RMS for points covering the near-equilibrium region) and the anharmonicity depends on the quartic derivative, which is then not very stable without higher order terms.

 

best regards,  -Kirk

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/SN6PR14MB21758668AC548770726CD102BD560%40SN6PR14MB2175.namprd14.prod.outlook.com.

[Paul Bagus]

Kirk, We agree. Your comments and cautions are well taken. Regards, Paul

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/10BAF264-7FD8-413F-B67F-8CAC5653227F%40wsu.edu.

[Hans Jørgen Aagaard Jensen]

How about using (1/R)**n for the fit? Have you tried that?

 

Also, have you tried Trond’s polfit in the utils/ directory?

 

-- Hans Jørgen

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/SN6PR14MB217536C3A9FCA8C1A31AEDECBD560%40SN6PR14MB2175.namprd14.prod.outlook.com.

[Peterson, Kirk]

I think years ago I compared different coordinate systems (1/R, SPF, etc.) and I've compared my results to others that have used expansions in 1/R, but after 30 years of using it, I'm pretty comfortable with powers of (r-re).  At least I understand it really well. In some ways the final results can be just as sensitive to how you sample the potential as to how you fit it.

 

I haven't tried Trond's polfit, but a long time ago I compared my program to Dunning's version and the results were very similar.

 

-Kirk

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/A4271C5C-B305-4F9B-923F-D6D9D60DE6B9%40sdu.dk.

[Paul Bagus]

Hans, What would be the advantage of an expansion in (1/R)**n. It is a way of reproducing an LJ 6-12 potential relevant for vdW interactions but this is a restricted class of systems. Also if you want properties near r_e, I believe you can expand the 6-12 in a polynomial. Isn’t this correct? Paul

 

From: 'Hans Jørgen Aagaard Jensen' via dirac-users <dirac...@googlegroups.com>

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/A4271C5C-B305-4F9B-923F-D6D9D60DE6B9%40sdu.dk.

[Paul Bagus]

Hans, BTW. My understanding is that polfit gives you the polynomial expansion and minimum and derivates around the minimum. Does it also compute the spectroscopic constants or does one have to use the formulas for them in terms of the expansion coefficients? Paul

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/E8E968D7-BBDD-48F4-8998-1D679AEC322A%40wsu.edu.

[Hans Jørgen Aagaard Jensen]

I guess the advantage of the 1/R expansion is that it will not go to alternatively to +inf and -inf for R -> inf and increasing order, but on the other hand it has problems for R -> 0.

 

I just searched for “inverse R fit of diatomic potential energy curve” with google and selected the hits

 

1. A thesis https://www.collectionscanada.gc.ca/obj/s4/f2/dsk3/OWTU/TC-OWTU-113.pdf  (supervisor Le Roy)
2. A chapter by Le Roy: http://scienide2.uwaterloo.ca/~rleroy/Pubn/11Demaison-Chapt6.pdf
3. A recent article: Fitting an experimental potential energy curve for the 10(0+)[43Π043Π0] electronic state of NaCs

J. Chem. Phys. 151, 024307 (2019);  https://doi.org/10.1063/1.5100748  

 

It seems that the best fits are obtained with a generalization of John Ogilvie’s suggestion to use  (R-Re)/(R+Re) as “x”, which is finite both at R=0 and at R-> inf.

From no. 1, the thesis:

 

This was followed by discussion of other fits.

 

Hmm, a thought: the nuclear repulsion is exactly known and can be treated analytically. Why not fit the electronic energy curve  which is finite at R=0 ???? One then obtains

V( R ) = V_NN( R ) + fit to E_el( R ). Would that not be more stable with the y_p variable (or another variable)? Has that been attempted ?

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/E8E968D7-BBDD-48F4-8998-1D679AEC322A%40wsu.edu.

[Peterson, Kirk]

If one wants to fit a whole diatomic potential curve, I think the best that I know of is from Bob LeRoy, the "modified Lennard-Jones oscillator" (MLJ) model. It uses the same "x" as you have below but with some very flexible options for introducing the correct long range (JCP 112, 3949).  I'm not aware of anyone doing what you suggest with the exact nuclear-nuclear repulsion but then the electronic energy just goes to a constant at R=0?  That might be a bit more difficult to enforce.

To view this discussion on the web visit https://groups.google.com/d/msgid/dirac-users/296F3FEA-4C91-4361-906D-27666B755456%40sdu.dk.

[Mrinal Thapa]

Dear Visscher,

In this case, what SCF value is correct for open-shell systems is correct to proceed with for further calculations, the SCF one or the RELCC module one. I am using it to calculate other properties. But I have been using the SCF energy from the HF.

Best Regards

Mrinal