In [1]:
import adcc
from pyscf import gto, scf
In [2]:
# Run SCF in pyscf
mol = gto.M(
    atom='O 0 0 0;'
         'H 0 0 1.795239827225189;'
         'H 1.693194615993441 0 -0.599043184453037',
    basis='cc-pvtz',
    unit="Bohr"
)
scfres = scf.RHF(mol)
scfres.conv_tol = 1e-12
scfres.conv_tol_grad = 1e-9
scfres.kernel()
refstate = adcc.ReferenceState(scfres)
converged SCF energy = -76.0571904154805
In [3]:
# Run ADC (2 and 2x for comparison)
state_2 = adcc.adc2(refstate, n_singlets=7)
print()
state_2x = adcc.adc2x(state_2.ground_state, n_singlets=7, guesses=state_2.excitation_vectors)
print()
state_3 = adcc.adc3(state_2x.ground_state, n_singlets=7, guesses=state_2x.excitation_vectors)
Starting adc2  singlet Jacobi-Davidson ...
Niter n_ss  max_residual  time  Ritz values
  1    14       0.39135   1.5s  [0.42513831 0.50590019 0.5059818  0.58483322 0.65566604 0.74712915
 0.76426908]
  2    28      0.013426   1.4s  [0.29512418 0.36981007 0.37674447 0.4528981  0.54751869 0.63327043
 0.63888939]
  3    42    0.00037195   1.6s  [0.2944929  0.36753152 0.37534001 0.45005037 0.54562195 0.63209821
 0.63769807]
  4    56    3.4634e-05   1.8s  [0.29448216 0.36746578 0.37521223 0.44984749 0.54547344 0.63193585
 0.6375989 ]
  5    70    9.4546e-07   2.1s  [0.29448182 0.36746378 0.3752094  0.44983771 0.54546393 0.63193085
 0.63759499]
=== Converged ===
    Number of matrix applies:    140
    Total solver time:             8s 428ms

Starting adc2x  singlet Jacobi-Davidson ...
Niter n_ss  max_residual  time  Ritz values
  1     7      0.014401   4.4s  [0.28228808 0.35677177 0.36353882 0.43995009 0.53537904 0.61680728
 0.62868375]
  2    14     0.0040087   4.5s  [0.27868225 0.3533497  0.3595459  0.43633214 0.53216524 0.61157475
 0.6255099 ]
  3    21    0.00035123   4.5s  [0.27795349 0.35277147 0.35881959 0.43571915 0.53151265 0.610253
 0.62495499]
  4    28    6.5036e-05   4.3s  [0.27786543 0.35270143 0.35871258 0.4356038  0.53142691 0.61014039
 0.62485118]
  5    35    1.6608e-05   4.6s  [0.27785588 0.35269586 0.35870035 0.43559008 0.53141433 0.6101313
 0.62482164]
  6    42    3.9341e-06   4.4s  [0.27785506 0.35269559 0.35869932 0.43558885 0.53141285 0.61013052
 0.62481454]
=== Restart ===
  7    14    7.0219e-07   4.5s  [0.27785499 0.35269557 0.35869921 0.43558868 0.53141263 0.61013044
 0.62481315]
=== Converged ===
    Number of matrix applies:    98
    Total solver time:            31s  40ms

Starting adc3  singlet Jacobi-Davidson ...
Niter n_ss  max_residual  time  Ritz values
  1     7     0.0065045   8.8s  [0.30766171 0.37989246 0.3899041  0.46360558 0.55137138 0.63738607
 0.64655262]
  2    14     0.0011243   8.3s  [0.30496765 0.37733362 0.3872384  0.46098356 0.54962477 0.63490646
 0.64453702]
  3    21    0.00030729   8.0s  [0.30438559 0.3768414  0.386693   0.46030024 0.54925434 0.63444585
 0.64402302]
  4    28    6.5818e-05   8.0s  [0.3043287  0.37680166 0.3866366  0.46021292 0.54922069 0.63441703
 0.64391612]
  5    35    6.5588e-06   8.3s  [0.30432329 0.37679956 0.38663105 0.46020415 0.54921711 0.63441506
 0.6438919 ]
  6    42    1.1205e-06   8.2s  [0.30432283 0.37679946 0.38663055 0.46020324 0.5492167  0.6344148
 0.64388927]
=== Restart ===
  7    14    2.3924e-07   8.0s  [0.3043228  0.37679946 0.38663051 0.46020317 0.54921665 0.63441476
 0.64388878]
=== Converged ===
    Number of matrix applies:    98
    Total solver time:            57s 565ms
In [4]:
# Show results
print(state_2.describe())
print(state_2x.describe())
print(state_3.describe())
+--------------------------------------------------------------+
| adc2                                    singlet ,  converged |
+--------------------------------------------------------------+
|  #        excitation energy     osc str    |v1|^2    |v2|^2  |
|          (au)           (eV)                                 |
|  0     0.2944818      8.013258   0.0385    0.9438   0.05623  |
|  1     0.3674638      9.999199   0.0000    0.9449   0.05514  |
|  2     0.3752094      10.20997   0.0992    0.9429   0.05707  |
|  3     0.4498377      12.24071   0.0524    0.9447   0.05534  |
|  4     0.5454639      14.84283   0.2498    0.9563   0.04366  |
|  5     0.6319309      17.19571   0.0000    0.9476   0.05241  |
|  6      0.637595      17.34984   0.0704    0.9524   0.04765  |
+--------------------------------------------------------------+

+--------------------------------------------------------------+
| adc2x                                   singlet ,  converged |
+--------------------------------------------------------------+
|  #        excitation energy     osc str    |v1|^2    |v2|^2  |
|          (au)           (eV)                                 |
|  0      0.277855      7.560819   0.0356    0.9268   0.07321  |
|  1     0.3526956      9.597335   0.0000    0.9295   0.07049  |
|  2     0.3586992      9.760703   0.0936    0.9258   0.07418  |
|  3     0.4355887      11.85297   0.0492    0.9292   0.07076  |
|  4     0.5314126      14.46047   0.2444    0.9421   0.05792  |
|  5     0.6101304      16.60249   0.0000    0.9235   0.07649  |
|  6     0.6248131      17.00203   0.0708    0.9386   0.06139  |
+--------------------------------------------------------------+

+--------------------------------------------------------------+
| adc3 (adc2)                             singlet ,  converged |
+--------------------------------------------------------------+
|  #        excitation energy     osc str    |v1|^2    |v2|^2  |
|          (au)           (eV)                                 |
|  0     0.3043228      8.281045   0.0378    0.9428    0.0572  |
|  1     0.3767995      10.25324   0.0000    0.9449    0.0551  |
|  2     0.3866305      10.52075   0.1000    0.9418   0.05822  |
|  3     0.4602032      12.52277   0.0570    0.9438   0.05623  |
|  4     0.5492167      14.94495   0.2533    0.9526   0.04741  |
|  5     0.6344148       17.2633   0.0000     0.938   0.06197  |
|  6     0.6438888      17.52111   0.0776    0.9497   0.05033  |
+--------------------------------------------------------------+

In [5]:
from matplotlib import pyplot as plt
plt.figure(figsize=(12, 7))
state_2.plot_spectrum(label="2")
state_2x.plot_spectrum(label="2x")
state_3.plot_spectrum(label="3")
plt.legend()
Out[5]:
<matplotlib.legend.Legend at 0x7fe7db7de898>