CISD Quadrupoles

I am trying to use psi4.core.variable("CISD QUADRUPOLE") to get the expectation values of quadrupole momments for the CISD ground state. But, I see the following error:

KeyError: 'psi4.core.variable: Requested variable CISD QUADRUPOLE was not set!\n'

Note: I call psi4.core.variable() only after running psi4.energy('cisd')

You never asked Psi to compute the quadrupole. Try psi4.properties('cisd').

1 Like

Thank you so much, I tried looking for this in the documentation and couldn’t find anything. But, I just realized I was looking in the wrong place as psi uses CC for CISD calculation and the CC page does mention about properties.

But, doesn’t psi4 compute the properties when asked to return the wavefunction?

No. Returning the wavefunction and computing properties are completely unrelated.

psi4.core.variable('CISD QUADRUPOLE') still returns variable not set error, when called after psi4.properties('cisd'). Strangely enough, when I checked the psi output file, it did print the Quadrupoles there.

Psi saves it as CI QUADRUPOLE.

It really should save as both, and I’ll put up an issue about it.

By the way, this is an unrelaxed property, not a relaxed property. There are two formulas for properties that are inequivalent for approximate methods. You’re using the cheaper but less accurate one (unrelaxed). The disagreement between the two is small when you include singles, but it’s good to be aware of this.

Actually, there seems to be a significant difference between the “relaxed” and “unrelaxed” properties.

Apart from that, while calculating the properties DETCI oddly seems to start from a very high electronic energy, for a Helium atom with a few ghost basis functions around it.
(See below the DETCI and CC output snippets)

         ---------------------------------------------------------
                          Configuration Interaction
                            (a 'D E T C I' module)

                 C. David Sherrill, Daniel G. A. Smith, and
                              Matt L. Leininger
         ---------------------------------------------------------


   ==> Parameters <==

    EX LEVEL       =        2      H0 BLOCKSIZE  =     1000
    VAL EX LEVEL   =        0      H0 GUESS SIZE =     1000
    H0COUPLINGSIZE =        0      H0 COUPLING   =       NO
    MAXITER        =       50      NUM PRINT     =       20
    NUM ROOTS      =        1      ICORE         =        1
    PRINT LVL      =        1      FCI           =       NO
    R CONV         = 1.00e-04      MIXED         =      YES
    E CONV         = 1.00e-12      MIXED4        =      YES
    R4S            =       NO      REPL OTF      =       NO
    DIAG METHOD    =      SEM      FOLLOW ROOT   =        0
    PRECONDITIONER = DAVIDSON      UPDATE        = DAVIDSON
    S              =   0.0000      Ms0           =      YES
    GUESS VECTOR   =  H0BLOCK      OPENTYPE      =     NONE
    COLLAPSE SIZE  =        1      HD AVG        = EVANGELISTI
    MAX NUM VECS   =       51      REF SYM       =     AUTO
    IOPEN        =       NO

    EX ALLOW       =  1  1 
    STATE AVERAGE  =  0(1.00) 

   ==> CI Orbital and Space information <==

    RAS1 LVL      =        0      A RAS3 MAX    =        1
    RAS1 MIN      =        0      B RAS3 MAX    =        1
    A RAS1 LVL    =        0      RAS4 LVL      =      126
    A RAS1 MIN    =        0      A RAS4 MAX    =        0
    A RAS1 MAX    =        1      B RAS4 MAX    =        0
    B RAS1 LVL    =        0      RAS4 MAX      =        0
    B RAS1 MIN    =        0      A RAS34 MAX   =        1
    B RAS1 MAX    =        1      B RAS34 MAX   =        1
    RAS3 LVL      =        1      RAS34 MAX     =        2
    RAS3 MAX      =        2
   ------------------------------------
               Space    Total     A
   ------------------------------------
                 Nso      142   142
                 Nmo      126   126
               Ndocc        1     1
               Nsocc        0     0
   ------------------------------------
              CI Spaces
   ------------------------------------
        Dropped DOCC        0     0
                RAS1        1     1
                RAS2        0     0
                RAS3      125   125
                RAS4        0     0
      Active (total)      126   126
        Dropped UOCC        0     0
   ------------------------------------

   ==> Setting up CI strings <==

    There are 126 alpha and 126 beta strings
    The CI space requires 15876 (1.59E+04) determinants and 4 blocks

   ==> Transforming CI integrals <==

	Presorting SO-basis two-electron integrals.
	Sorting File: SO Ints (nn|nn) nbuckets = 1
	Constructing frozen core operators
	Starting first half-transformation.
	Sorting half-transformed integrals.
	First half integral transformation complete.
	Starting second half-transformation.
	Two-electron integral transformation complete.

   ==> Starting CI iterations <==

    H0 Block Eigenvalue = -447005.72378967

    Simultaneous Expansion Method (Block Davidson Method)
    Using 1 initial trial vectors

     Iter   Root       Total Energy       Delta E      C RMS

   @CI  0:     0  -447005.723826216068   -4.4701E+05   1.6537E+03  
   @CI  1:     0  -447014.629862402566   -8.9060E+00   1.7494E+03  
   @CI  2:     0  -894016.679259274621   -4.4700E+05   4.3103E+02  
   @CI  3:     0  -894017.069101640605   -3.8984E-01   1.0213E+02  
   @CI  4:     0  -894017.081594379270   -1.2493E-02   3.5157E-01  
   @CI  5:     0  -894017.081594656222   -2.7695E-07   8.6757E-03  
    Warning: Norm of correction (root 0) is < 1.0E-13
   @CI  6:     0  -894017.081594655989   2.3283E-10   1.9923E-04  
    Warning: Norm of correction (root 0) is < 1.0E-13
   @CI  7:     0  -894017.081594655407   5.8208E-10   3.9912E-07  
    Warning: Norm of correction (root 0) is < 1.0E-13
   @CI  8:     0  -894017.081594656804   -1.3970E-09   1.6994E-07  
    Warning: Norm of correction (root 0) is < 1.0E-13
   @CI  9:     0  -894017.081594655989   8.1491E-10   1.4811E-07  
    Warning: Norm of correction (root 0) is < 1.0E-13
   @CI 10:     0  -894017.081594656571   -5.8208E-10   1.1319E-06  
    Warning: Norm of correction (root 0) is < 1.0E-13
   @CI 11:     0  -894017.081594656571   0.0000E+00   3.3058E-06 c

   ==> Energetics <==

    SCF energy =           -2.861532850861122
    Total CI energy =    -894017.081594656570815

   ==> CISD root 0 information <==

    CISD Root 0 energy =  -894017.081594656570815

   The 20 most important determinants:

    *   1    0.560280  (   84,   84)  85AX 
    *   2    0.434146  (   83,   84)  84AA 85AB 
    *   3    0.434146  (   84,   83)  84AB 85AA 
    *   4    0.336409  (   83,   83)  84AX 
    *   5    0.202332  (   82,   84)  83AA 85AB 
    *   6    0.202332  (   84,   82)  83AB 85AA 
    *   7    0.156782  (   82,   83)  83AA 84AB 
    *   8    0.156782  (   83,   82)  83AB 84AA 
    *   9    0.084703  (   75,   84)  76AA 85AB 
    *  10    0.084703  (   84,   75)  76AB 85AA 
    *  11    0.073068  (   82,   82)  83AX 
    *  12    0.071591  (   70,   84)  71AA 85AB 
    *  13    0.071591  (   84,   70)  71AB 85AA 
    *  14    0.065634  (   75,   83)  76AA 84AB 
    *  15    0.065634  (   83,   75)  76AB 84AA 
    *  16    0.055474  (   70,   83)  71AA 84AB 
    *  17    0.055474  (   83,   70)  71AB 84AA 
    *  18    0.052377  (   71,   84)  72AA 85AB 
    *  19    0.052377  (   84,   71)  72AB 85AA 
    *  20    0.040586  (   71,   83)  72AA 84AB 

        *******************************************************
        *                                                     *
        *                        CISD                         *
        *      Singles Doubles Configuration Interaction      *
        *                                                     *
        *                   Eugene DePrince                   *
        *                                                     *
        *******************************************************



  ==> Input parameters <==

        Freeze core orbitals?                  no
        Use frozen natural orbitals?           no
        r_convergence:                  1.000e-07
        e_convergence:                  1.000e-12
        Number of DIIS vectors:                 8
        Number of frozen core orbitals:         0
        Number of active occupied orbitals:     1
        Number of active virtual orbitals:    125
        Number of frozen virtual orbitals:      0

  ==> Memory <==

        available memory =                           1907.35 mb
        minimum memory requirements for QCISD =         0.60 mb
        (T) algorithm:                                  0.24 mb (low-memory)
        memory requirements for QCISD(T) =              0.24 mb

  ==> Define tiling <==

        v(ab,cd) diagrams will be evaluated in   1 blocks.
        v(ab,ci) diagrams will be evaluated in   1 blocks over ov2.
        v(ab,ci) diagrams will be evaluated in   1 blocks over ov.

  Allocate cpu memory (   473.74 mb).....done.
  Initialize cpu memory..................done.

  Begin CISD iterations

   Iter  DIIS          Energy       d(Energy)          |d(T)|     time
      0   0 1   -0.0357337652   -0.0357337652    0.0747914141        0
      1   1 1   -0.0398195123   -0.0040857471    0.0129344286        0
      2   2 1   -0.0409419885   -0.0011224762    0.0033528833        0
      3   3 1   -0.0410169652   -0.0000749766    0.0004093105        1
      4   4 1   -0.0410178110   -0.0000008458    0.0000798056        0
      5   5 1   -0.0410178359   -0.0000000249    0.0000152596        0
      6   6 1   -0.0410178368   -0.0000000009    0.0000032928        1
      7   7 1   -0.0410178368   -0.0000000000    0.0000005766        0
      8   8 1   -0.0410178368   -0.0000000000    0.0000000839        0
      9   8 2   -0.0410178368   -0.0000000000    0.0000000162        1

  CISD iterations converged!

  CISD variational energy:      -0.041017836827
  CISD transition energy:       -0.041017836836

        OS SCS-MP2 correlation energy:          -0.042880518290
        SS SCS-MP2 correlation energy:           0.000000000000
        SCS-MP2 correlation energy:             -0.042880518290
      * SCS-MP2 total energy:                   -2.904413369151

        OS MP2 correlation energy:              -0.035733765242
        SS MP2 correlation energy:               0.000000000000
        MP2 correlation energy:                 -0.035733765242
      * MP2 total energy:                       -2.897266616103

        OS CISD correlation energy:             -0.041017836836
        SS CISD correlation energy:              0.000000000000
        CISD correlation energy:                -0.041017836827
      * CISD total energy:                      -2.902550687688

The detci energy is clearly incorrect, but this is unrelated to unrelaxed vs. relaxed dipoles. Two codes that are supposed to give the same answer do not give the same answer, seemingly because one is buggy.

Please post a complete example of your input to the detci computation, and I can refer this to the specialists.

Here’s my psi4 input in python. See the attachments provided below, for the .xyz file, where Gh(He1)'s are ghost atoms with a diffuse gaussian. All these functions are placed at points on a spherically grid as the one shown in the image.

He.xyz (7.0 KB)

grid

import psi4
import numpy as np
# psi4 settings
psi4.core.set_num_threads(8)
psi4.set_memory('16 Gb')
options_dict = {'scf_type': 'pk',
                'guess': 'sad',
                'e_convergence': 1e-12,
                'd_convergence': 1e-8,
                'ints_tolerance': 1e-12,
                'maxiter' : 200,
                'ci_maxiter':50,
                'ms0': True}
psi4.set_options(options_dict)
psi4.set_output_file('cisd_properties_test.out', False)
# Setting custom diffuse gaussians on ghost atoms
def basisspec_psi4_yo__mybasis(mol, role):
    basstrings = {}
    mol.set_basis_all_atoms('aug-cc-pvqz', role=role)
    mol.set_basis_by_label("He1", "gh-bs", role=role)
    basstrings['gh-bs']= """
    cartessian
    ****
    He     0
    S    1   1.00
        0.08665577       1.00000000
    ****
    """
    return basstrings
psi4.qcdb.libmintsbasisset.basishorde['MYBASIS'] = basisspec_psi4_yo__mybasis
psi4.core.set_global_option('basis', 'mybasis')
# Reading mol file and running calculations
with open('He.xyz', 'r') as file:
    molstr = file.read()
mol = psi4.geometry(molstr)
e_cisd, wfn_cisd = psi4.energy('cisd', return_wfn=True)
psi4.properties('cisd')

This input file segfaults when I try it. It looks like your psi4.energy call is Garbage In, Garbage Out, where the psi4.properties call is fine. I’ll file an issue, but I can’t afford the time to fix this one.