Trovemaster / Duo

Duo is a diatomic code for solving a fully coupled rovibronic Schroedinger equation
14 stars 11 forks source link

N.A. in output #2

Closed MartinBeseda closed 5 years ago

MartinBeseda commented 5 years ago

Hello, I'm a completely new user of Duo and I'm puzzled by one thing.

I wanted to compute energy levels and spectroscopic constants of N2+ in SigmaU+ symmetry.

The problem is, I'm not getting N.A. symbol instead of the desired spectroscopic symols.

Equilibrium properties
           imin =                      6
   r_imin / ang =           1.0595959596
V(r_imin) / cm-1 =         -391.89938844
 Has a single min?                    No
 True r_e  / ang =                  N.A.

Derivatives at true re
der0, cm-1       =                  N.A.
der1, cm-1/ang   =                  N.A.
der2, cm-1/ang^2 =                  N.A.
der3, cm-1/ang^3 =                  N.A.
der4, cm-1/ang^4 =                  N.A.

Harmonic we, cm-1=                  N.A.
  Rotat. B0, cm-1=                  N.A.
Anharm. const. xe=                  N.A.
 Coriol. ae, cm-1=                  N.A.
  Centr. De, cm-1=                  N.A.
        Y00, cm-1=                  N.A.

Approximate J=0 vibrational energy levels (no couplings) 
 given by E(v, J=0) = V(re) + Y00 + we*(v+0.5) - we*xe*(v+0.5)^2
                 v
                 0                  N.A.
                 1                  N.A.
                 2                  N.A.
                 3                  N.A.

            Spin =                   0.5
        |Lambda| =                   1.0
  Physical J_min =                   0.5

Approximate J=J_min vibrational-rotational energy levels (no couplings)
 given by E(v, J) = E(v, J=0) B0*J*(J+1)- ae*(v+0.5)*J*(J+1) - De*[J(J+1)]^2
                 v
                 0                  N.A.
                 1                  N.A.
                 2                  N.A.

I suppose, that it's caused by a detection of multiple minima in my data. The problem is, there should be just one real minima, as can be seen in the plot. sigmau+

Could you, please, tell me, what am I doing wrong in the input?


Complete output file

(Transcript of the input --->)
(N2+ 2SigmaU+)
masses 14 14
nstates 1
jrot -2 2
symmetry C2h

grid
    npoints 100
    range 0.85 5
    type 0
end

vibrationalbasis
    vmax 1
end

eigensolver
    nroots 20
end

potential 1
    name "2 SigmaU+"
    symmetry u +
    lambda 1
    mult 2
    spin 0.5
    type grid
    interpolationtype Cubicsplines
    units angstrom eH
    values
        0.8 0.4663431999999972
        0.85    0.2716072000000054
        0.9 0.14446329999999818
        0.95    0.06631509999999707
        1   0.023365799999993442
        1.05    0
        1.1 -0.00121289999999874
        1.15    0.008618799999993598
        1.2 0.0252145999999982
        1.25    0.04443120000000533
        1.3 0.06253599999999437
        1.35    0.0777470999999963
        1.4 0.08978539999999668
        1.45    0.099500699999993
        1.5 0.1078013999999996
        1.55    0.11524179999999262
        1.6 0.12209389999999587
        1.65    0.1284830000000028
        1.7 0.13446980000000508
        1.75    0.14007789999999432
        1.8 0.14535119999999324
        1.85    0.1503034999999926
        1.9 0.15495169999999803
        1.95    0.1593098000000026
        2   0.16338980000000447
        2.05    0.16720200000000318
        2.1 0.17075570000000084
        2.15    0.17405920000000208
        2.2 0.1771202000000045
        2.25    0.1799457999999987
        2.3 0.18254299999999546
        2.35    0.18491869999999722
        2.4 0.18708049999999332
        2.45    0.1890367000000026
        2.5 0.19079659999999876
        2.55    0.19237090000000023
        2.6 0.1937710999999922
        2.65    0.1950096000000059
        2.7 0.1960983999999968
        2.75    0.1970500000000044
        2.8 0.19787619999999606
        2.85    0.198588799999996
        2.9 0.1991991000000013
        2.95    0.19971769999999367
        3   0.20015519999999754
        3.05    0.20054960000000221
        3.1 0.20084439999999404
        3.15    0.2010855000000049
        3.2 0.20128060000000403
        3.25    0.20143659999999386
        3.3 0.2015597999999983
        3.35    0.20165559999999516
        3.4 0.20172870000000387
        3.45    0.20178300000000604
        3.5 0.20182219999999518
        3.55    0.20184899999999573
        3.6 0.20186619999999778
        3.65    0.20167109999999866
        3.7 0.20167419999999936
        3.75    0.20167290000000548
        3.8 0.20166840000000263
        3.85    0.2016617000000025
        3.9 0.20165350000000615
        3.95    0.20164450000000045
        4   0.20163519999999835
        4.05    0.20162589999999625
        4.1 0.20161690000000476
        4.15    0.20160850000000607
        4.2 0.20160079999999425
        4.25    0.20159400000000005
        4.3 0.201588000000001
        4.35    0.2015829999999994
        4.4 0.20157890000000123
        4.45    0.20157580000000053
        4.55    0.20157229999999515
        4.6 0.20157179999999641
        4.65    0.20157209999999282
        4.7 0.2015731000000045
        4.75    0.20157480000000305
        4.8 0.20157720000000268
        4.85    0.20157999999999276
        4.9 0.20158340000000408
        4.95    0.2015872999999999
        5   0.20159160000000043
        30  0.20215919999999699
end
(<--- End of the input.)

 ____                                        ____  _   _  ___
|  _ \ _ __ ___   __ _ _ __ __ _ _ __ ___   |  _ \| | | |/ _ \
| |_) | '__/ _ \ / _` | '__/ _` | '_ ` _ \  | | | | | | | | | |
|  __/| | | (_) | (_| | | | (_| | | | | | | | |_| | |_| | |_| |
|_|   |_|  \___/ \__, |_|  \__,_|_| |_| |_| |____/ \___/ \___/
                 |___/

Please refer to:
 Sergei N. Yurchenko, Lorenzo Lodi, Jonathan Tennyson and Andrey V. Stolyarov, 
 `DUO: a general program for calculating spectra of diatomic molecules',
  Computer Physics Communication, (to be submitted), 2015.
  Contacts: s.yurchenko@ucl.ac.uk; l.lodi@ucl.ac.uk; j.tennyson@ucl.ac.uk; 
            avstol@phys.chem.msu.ru

Values of physical constants used by DUO:
                    Planck constant h =   6.626069570000E-27    erg*second
                     Speed of light c =     29979245800.0000     cm/second
                       Bohr radius a0 =      0.5291772109200     angstroms
                    Hartree energy Eh =    219474.6313708000         cm^-1
           Unified atomic mass unit u =   1.660538921000E-24         grams
           Unified atomic mass unit u =   1822.8884861185961            me
                    Electron charge e =   1.602176565000E-19      Coulombs
                Boltzmann constant kB =   1.380648800000E-16    erg/Kelvin
                    Avogadro constant =   6.022141290000E+23        mol^-1
      Fine structure constant 1/alpha =   137.03599907466793
                     Electron mass me =   9.109382903261E-28         grams
                  a.u. of dipole e*a0 =       2.541746363812        debyes
                                 1 eV =      8065.5442959967         cm^-1

Reduced mass is      7.00000000000000         atomic mass units (Daltons)
Reduced mass is      12760.2194028302         atomic units (electron masses)

Using a uniformly spaced grid with npoints =             100
rmin, rmax, step (in ang)   =          0.85000000000000         5.00000000000000         0.04191919191919
rmin, rmax, step (in bohrs) =          1.60626720588030         9.44863062282531         0.07921579209035

Generate a grid representation for all Hamiltonian terms
...done!

Potential functions:
     1                       2 SigmaU+  type =                 GRID
            r(Ang)                     1
        0.85000000        5.96108901E+04
        0.89191919        3.53189271E+04
        0.93383838        1.91364265E+04
        0.97575758        8.99058349E+03
        1.01767677        2.82452296E+03
        1.05959596       -3.91899388E+02
        1.10151515       -2.25374736E+02
        1.14343434        1.51274429E+03
        1.18535354        4.36117328E+03
        1.22727273        7.82374126E+03
        1.26919192        1.13348204E+04
        1.31111111        1.45283006E+04
        1.35303030        1.72424716E+04
        1.39494949        1.94656822E+04
        1.43686869        2.13149816E+04
        1.47878788        2.29153448E+04
        1.52070707        2.43544961E+04
        1.56262626        2.56831723E+04
        1.60454545        2.69279109E+04
        1.64646465        2.81026012E+04
        1.68838384        2.92150897E+04
        1.73030303        3.02679469E+04
        1.77222222        3.12667133E+04
        1.81414141        3.22153268E+04
        1.85606061        3.31149423E+04
        1.89797980        3.39680011E+04
        1.93989899        3.47762295E+04
        1.98181818        3.55412089E+04
        2.02373737        3.62643273E+04
        2.06565657        3.69468352E+04
        2.10757576        3.75898847E+04
        2.14949495        3.81945240E+04
        2.19141414        3.87617458E+04
        2.23333333        3.92924672E+04
        2.27525253        3.97875985E+04
        2.31717172        4.02480264E+04
        2.35909091        4.06746694E+04
        2.40101010        4.10685374E+04
        2.44292929        4.14307117E+04
        2.48484848        4.17623801E+04
        2.52676768        4.20648984E+04
        2.56868687        4.23397119E+04
        2.61060606        4.25883754E+04
        2.65252525        4.28124940E+04
        2.69444444        4.30136088E+04
        2.73636364        4.31933546E+04
        2.77828283        4.33532576E+04
        2.82020202        4.34948664E+04
        2.86212121        4.36196226E+04
        2.90404040        4.37290901E+04
        2.94595960        4.38245410E+04
        2.98787879        4.39065934E+04
        3.02979798        4.39826049E+04
        3.07171717        4.40460431E+04
        3.11363636        4.40955668E+04
        3.15555556        4.41384255E+04
        3.19747475        4.41740461E+04
        3.23939394        4.42035995E+04
        3.28131313        4.42279262E+04
        3.32323232        4.42477068E+04
        3.36515152        4.42636507E+04
        3.40707071        4.42762234E+04
        3.44898990        4.42860332E+04
        3.49090909        4.42937717E+04
        3.53282828        4.42976945E+04
        3.57474747        4.43073206E+04
        3.61666667        4.42916287E+04
        3.65858586        4.42583558E+04
        3.70050505        4.42624477E+04
        3.74242424        4.42625506E+04
        3.78434343        4.42612499E+04
        3.82626263        4.42604591E+04
        3.86818182        4.42589779E+04
        3.91010101        4.42574437E+04
        3.95202020        4.42557702E+04
        3.99393939        4.42540594E+04
        4.03585859        4.42523441E+04
        4.07777778        4.42506594E+04
        4.11969697        4.42490515E+04
        4.16161616        4.42475431E+04
        4.20353535        4.42461493E+04
        4.24545455        4.42448973E+04
        4.28737374        4.42437651E+04
        4.32929293        4.42427836E+04
        4.37121212        4.42419480E+04
        4.41313131        4.42412548E+04
        4.45505051        4.42407173E+04
        4.49696970        4.42403199E+04
        4.53888889        4.42400550E+04
        4.58080808        4.42399173E+04
        4.62272727        4.42399059E+04
        4.66464646        4.42400113E+04
        4.70656566        4.42402216E+04
        4.74848485        4.42405411E+04
        4.79040404        4.42409728E+04
        4.83232323        4.42414664E+04
        4.87424242        4.42420418E+04
        4.91616162        4.42427079E+04
        4.95808081        4.42434465E+04
        5.00000000        4.42442421E+04

Equilibrium properties
           imin =                      6
   r_imin / ang =           1.0595959596
V(r_imin) / cm-1 =         -391.89938844
 Has a single min?                    No
 True r_e  / ang =                  N.A.

Derivatives at true re
der0, cm-1       =                  N.A.
der1, cm-1/ang   =                  N.A.
der2, cm-1/ang^2 =                  N.A.
der3, cm-1/ang^3 =                  N.A.
der4, cm-1/ang^4 =                  N.A.

Harmonic we, cm-1=                  N.A.
  Rotat. B0, cm-1=                  N.A.
Anharm. const. xe=                  N.A.
 Coriol. ae, cm-1=                  N.A.
  Centr. De, cm-1=                  N.A.
        Y00, cm-1=                  N.A.

Approximate J=0 vibrational energy levels (no couplings) 
 given by E(v, J=0) = V(re) + Y00 + we*(v+0.5) - we*xe*(v+0.5)^2
                 v
                 0                  N.A.
                 1                  N.A.
                 2                  N.A.
                 3                  N.A.

            Spin =                   0.5
        |Lambda| =                   1.0
  Physical J_min =                   0.5

Approximate J=J_min vibrational-rotational energy levels (no couplings)
 given by E(v, J) = E(v, J=0) B0*J*(J+1)- ae*(v+0.5)*J*(J+1) - De*[J(J+1)]^2
                 v
                 0                  N.A.
                 1                  N.A.
                 2                  N.A.
                 3                  N.A.

Construct the J=0 matrix
Solving one-dimentional Schrodinger equations using : SINC

Vibrational (contracted) energies: 
    N        Energy/cm    State v

    1          0.000000 [    1   0 ] 2 SigmaU+

j =       2.0

Define the quanta book-keeping
The total number sigma/lambda states (size of the sigma-lambda submatrix) = 4

Sigma-Lambda basis set:
     i     jrot  state   spin    sigma lambda   omega
     1      2.0    1      0.5     -0.5    1      0.5   2 SigmaU+
     2      2.0    1      0.5      0.5    1      1.5   2 SigmaU+
     3      2.0    1      0.5     -0.5   -1     -1.5   2 SigmaU+
     4      2.0    1      0.5      0.5   -1     -0.5   2 SigmaU+
...done!

Contracted basis set:
     i     jrot ilevel ivib state v     spin    sigma lambda   omega   Name
     1      2.0    1    1    1    0      0.5     -0.5    1      0.5   2 SigmaU+
     2      2.0    2    1    1    0      0.5      0.5    1      1.5   2 SigmaU+
     3      2.0    3    1    1    0      0.5     -0.5   -1     -1.5   2 SigmaU+
     4      2.0    4    1    1    0      0.5      0.5   -1     -0.5   2 SigmaU+

 Memory Report:
 Active Arrays                                        size (GiB)
 GRID                    0.000004                                                
 r-field                 0.000003                                                
 contrfunc               0.000075                                                
 Total memory   =                                     0.00008427 GiB
 Maximal memory =                                     0.00020479 GiB (        100000.0)

 (       19 arrays contributing less than 1% are not shown)

Construct the hamiltonian matrix
...done!

Non-zero matrix elements of the coupled Sigma-Lambda matrix:
Threshold for printing is   2.22E-16 cm^-1
RV == Rotational-vibrational
SO == Spin-Orbit interaction
SS == Spin-Spin interaction
JS == J.S interaction (aka S-uncoupling)
LS == L.J interaction (aka L-uncoupling)
LS == L.S interaction (spin-electronic)
RV=        677.001;  JS(  1  2)=     -4.7036 ;
RV=        672.895;
RV=        672.895;  JS(  3  4)=     -4.7036 ;
RV=        677.001;

Eigenvalues for J =      2.0

       J      N        Energy/cm  State   v  lambda spin   sigma   omega  parity
       2.0    1        669.816022   1     0   1     0.5     0.5     1.5   +    ||2 SigmaU+
       2.0    2        680.080211   1     0   1     0.5    -0.5     0.5   +    ||2 SigmaU+

       J      N        Energy/cm  State   v  lambda spin   sigma   omega  parity
       2.0    1        669.816022   1     0  -1     0.5    -0.5    -1.5   -    ||2 SigmaU+
       2.0    2        680.080211   1     0  -1     0.5     0.5    -0.5   -    ||2 SigmaU+

Zero point energy (ZPE) =           0.000000

 Memory Report:
 Active Arrays                                        size (GiB)
 GRID                    0.000004                                                
 r-field                 0.000003                                                
 contrfunc               0.000075                                                
 Total memory   =                                     0.00008338 GiB
 Maximal memory =                                     0.00020479 GiB (        100000.0)

 (       25 arrays contributing less than 1% are not shown)

              Timing data at 2019/04/23 15:38:11

                                            Total time (seconds)  Self time (seconds)
 Timer                               Calls  --------------------  -------------------
 -----                               -----       Real       CPU        Real       CPU

 Duo                                    1.        0.0       0.0         0.0       0.0
 Map on grid                            1.        0.0       0.0         0.0       0.0
 Solve vibrational part                 1.        0.0       0.0         0.0       0.0

 (  7 timers contributing less than 1% are not shown)
Trovemaster commented 5 years ago

Dear Martin,

This feature is only working for analytical representations. It has not being implemented for the grid-type functions, sorry. It is on my todo list.

Cheers, Sergey

MartinBeseda commented 5 years ago

Dear Sergey,

so, if I understand you well, you don't support piecewise-defined functions currently, do you? I was thinking about interpolating the dataset manually and providing the resulting function to DUO.

Cheers, Martin

Trovemaster commented 5 years ago

We do support piecewise-defined functions all right, just not for this feature of estimating the spectroscopic constants. The latter have never being a priority as the main purpose of Duo is a direct solution of the Schroedinger equation. Besides, spectroscopic constants, at least of the higher order, are heavily empirical. There is no much sense, from my experience, deriving them from the potential energy curves.

Cheers, Sergey

MartinBeseda commented 5 years ago

Well, I see the sense of deriving them from PESs in two things: 1) Checking the PES comparing the constants with already known experients 2) Providing them to experimentalists, if they have no previous results.

And just off-topic question - which software did you use for the diatomic spectroscopic constants computations then?

Trovemaster commented 5 years ago

Yes, I am aware of these applications. (1) Personally I prefer comparing energies not constants. (2) I do see sometimes (not very often) experimentalists trying to use ab initio spectroscopic constants, but not with much success.

RE: spectroscopic constants, I don't compute them at all. Spectroscopic constants only make sense as part of the perturbative/empirical methodology based on the effective rotational Hamiltonians for computing energies or line positions. We tend to use the variational methodology, based on the PECs (spin-orbit curves etc) directly. If I need to do a effective rotational Hamoltonian treatment I use PGOPHER.

MartinBeseda commented 5 years ago

Ok, I understand. Still, I'd like to compute them at least for this one case, but I understand, that it isn't your priority. Thank you for the prompt responses!

Martin

Trovemaster commented 5 years ago

I will surely introduce this feature at some point. As a quick fix I can suggest to represent your potential energy curve as, e.g. a Morse potential, or even as a polynomial, or something like this poten 1 name "X2Sigmau+" symmetry u + lambda 1 mult 2 type EMO values V0 0.00000000000000E+00
RE 1.08435981870765E+00 DE 4.44521643363054E+04 RREF -1.00000000000000E+00 PL 8.00000000000000E+00 PR 8.00000000000000E+00 NL 1.00000000000000E+00 NR 3.00000000000000E+00 B0 3.53434543259561E+00 fit B1 3.41588567189263E-01 fit B2 2.36501323107290E+00 fit B3 -3.81206035003490E+00 fit end

The spectroscopic constants only depend on the firs few derivatives, i.e. it is only the very bottom of the potential you need to fit accurately. However, looking at you potential, you do want to have more points around the equilibrium (below ~5000 cm-1) in order to have spectroscopic constants accurate. Otherwise the uncertainly will be too high. You can try using Duo to fit to ab initio data, see example below, but it will only work with a better coverage around equilibrium and accurate (smooth) ab initio data. A side note, you should use accurate masses or switch to "atoms" (see below).

Cheers, Sergey

(N2+ 2SigmaU+)

atoms N N

nstates 1

jrot -2 2

symmetry C2v

grid npoints 100 range 0.85 5 type 0 end

vibrationalbasis vmax 1 end

eigensolver nroots 50 end

poten 1 name "X2Sigmau+" symmetry u + lambda 1 mult 2 type EMO values V0 0.00000000000000E+00
RE 1.08435981870765E+00 DE 4.44521643363054E+04 RREF -1.00000000000000E+00 PL 8.00000000000000E+00 PR 8.00000000000000E+00 NL 1.00000000000000E+00 NR 3.00000000000000E+00 B0 3.53434543259561E+00 fit B1 3.41588567189263E-01 fit B2 2.36501323107290E+00 fit B3 -3.81206035003490E+00 fit end

DO_NOT_SHIFT_PECS

abinitio potential 1 name "2 SigmaU+" symmetry u + lambda 1 mult 2 spin 0.5 type grid interpolationtype Cubicsplines units angstrom eH weighting ps1997 0.005 50000.0 values 0.8 0.4663431999999972 0.85 0.2716072000000054 0.9 0.14446329999999818 0.95 0.06631509999999707 1 0.023365799999993442 1.05 0 1.1 -0.00121289999999874 1.15 0.008618799999993598 1.2 0.0252145999999982 1.25 0.04443120000000533 1.3 0.06253599999999437 1.35 0.0777470999999963 1.4 0.08978539999999668 1.45 0.099500699999993 1.5 0.1078013999999996 1.55 0.11524179999999262 1.6 0.12209389999999587 1.65 0.1284830000000028 1.7 0.13446980000000508 1.75 0.14007789999999432 1.8 0.14535119999999324 1.85 0.1503034999999926 1.9 0.15495169999999803 1.95 0.1593098000000026 2 0.16338980000000447 2.05 0.16720200000000318 2.1 0.17075570000000084 2.15 0.17405920000000208 2.2 0.1771202000000045 2.25 0.1799457999999987 2.3 0.18254299999999546 2.35 0.18491869999999722 2.4 0.18708049999999332 2.45 0.1890367000000026 2.5 0.19079659999999876 2.55 0.19237090000000023 2.6 0.1937710999999922 2.65 0.1950096000000059 2.7 0.1960983999999968 2.75 0.1970500000000044 2.8 0.19787619999999606 2.85 0.198588799999996 2.9 0.1991991000000013 2.95 0.19971769999999367 3 0.20015519999999754 3.05 0.20054960000000221 3.1 0.20084439999999404 3.15 0.2010855000000049 3.2 0.20128060000000403 3.25 0.20143659999999386 3.3 0.2015597999999983 3.35 0.20165559999999516 3.4 0.20172870000000387 3.45 0.20178300000000604 3.5 0.20182219999999518 3.55 0.20184899999999573 3.6 0.20186619999999778 3.65 0.20167109999999866 3.7 0.20167419999999936 3.75 0.20167290000000548 3.8 0.20166840000000263 3.85 0.2016617000000025 3.9 0.20165350000000615 3.95 0.20164450000000045 4 0.20163519999999835 4.05 0.20162589999999625 4.1 0.20161690000000476 4.15 0.20160850000000607 4.2 0.20160079999999425 4.25 0.20159400000000005 4.3 0.201588000000001 4.35 0.2015829999999994 4.4 0.20157890000000123 4.45 0.20157580000000053 4.55 0.20157229999999515 4.6 0.20157179999999641 4.65 0.20157209999999282 4.7 0.2015731000000045 4.75 0.20157480000000305 4.8 0.20157720000000268 4.85 0.20157999999999276 4.9 0.20158340000000408 4.95 0.2015872999999999 5 0.20159160000000043 30 0.20215919999999699 end

FITTING JLIST 0 itmax 20 fit_factor 1e-12 robust 0.0 output ESigma_Elander energies (J parity NN energy ) (e-state v ilambda isigma omega weight) 0 + 1 0.00 1 0 1 0 0 1.0 end