Sorting of CCSD amplitudes

[bsenjean]

Hi,

I am currently checking if the CCSD amplitudes extracted from Psi4 are sorted correctly in OpenFermion, and I have the impression that something is wrong (but maybe I'm doing something wrong). Here is an example for H4 (I enclose the whole script below if you want to run it).

I find that in Psi4, the CCSD amplitudes are the following:

    Largest TIjAb Amplitudes:
      1   1   1   1        -0.2129966367
      1   1   0   0        -0.1959962733
      0   1   1   0        -0.1937183230
      1   0   0   1        -0.1937183230
      0   1   0   1        -0.1876273767
      1   0   1   0        -0.1876273767
      0   0   1   1        -0.1861373089
      0   0   0   0        -0.1697814772

Where the indices are: "I" --> occupied alpha spin-orbital (labelled from 0 to 1) "j" --> occupied beta spin-orbital (labelled from 0 to 1) "A" --> virtual alpha spin-orbital (labelled from 0 to 1) "b" --> virtual beta spin-orbital (labelled from 0 to 1).

So if I'm not mistaken, those spin-orbitals have the following indices in OpenFermion: I = 0 --> 0 I = 1 --> 2 j = 0 --> 1 j = 1 --> 3 A = 0 --> 4 A = 1 --> 6 b = 0 --> 5 b = 1 --> 7

So the amplitude -0.2129966367 should correspond to `2 3 6^ 7^'.

However, if I look at the (normal ordered) ucc amplitudes generated by OpenFermion, I get the following:

-0.1861373089 [1^ 0^ 5 4] +
-0.1697814772 [1^ 0^ 7 6] +
-0.193718323 [2^ 1^ 6 5] +
0.1876273767 [2^ 1^ 7 4] +
0.1876273767 [3^ 0^ 6 5] +
-0.193718323 [3^ 0^ 7 4] +
-0.2129966367 [3^ 2^ 5 4] +
-0.1959962733 [3^ 2^ 7 6] +
0.1861373089 [5^ 4^ 1 0] +
0.2129966367 [5^ 4^ 3 2] +
0.193718323 [6^ 5^ 2 1] +
-0.1876273767 [6^ 5^ 3 0] +
-0.1876273767 [7^ 4^ 2 1] +
0.193718323 [7^ 4^ 3 0] +
0.1697814772 [7^ 6^ 1 0] +
0.1959962733 [7^ 6^ 3 2]

And so at the end it's like the indices for the virtual orbitals are swapped, i.e. 4 <--> 6 and 5 <--> 7. I checked in more details and it seems the indices are already swapped (according to me) in molecule.ccsd_double_amplitudes.

I enclose a script to test it below. Am I missing something obvious or not ? Thanks !

import psi4
import openfermion
import openfermionpsi4
import numpy as np
import os

data_dir = os.getcwd()
theta = 60
radius = 1.5 #Angstrom
phi = theta*0.5*np.pi/180. #phi = 1/2*theta. Convert to radians
x = radius*np.cos(phi)
y = radius*np.sin(phi)
basis = "sto-3g"
multiplicity = 1
description = str(theta)
geometry = [('H', (x,y, 0.)), ('H', (x, -y, 0.)), ('H', (-x,y, 0.)), ('H', (-x, -y, 0.))]

# TEST PSI4
psi4.geometry('''
H  {0}  {1}  0.
H  {0} -{1}  0.
H -{0}  {1}  0.
H -{0} -{1}  0.
'''.format(x,y))
psi4.core.set_output_file("ccsd_H4.txt", False)
psi4.set_options({"basis": basis,
                  "mp2_amps_print":True,
                  "num_amps_print":26})
cc_energy, cc_wfn = psi4.energy('ccsd',return_wfn = True)

# TEST OPENFERMION
molecule = openfermion.MolecularData(
    geometry,
    basis,
    multiplicity,
    description=description,
    data_directory=data_dir)

molecule = openfermionpsi4.run_psi4(molecule,
                    run_ccsd=1)

hamiltonian, one_body, two_body = molecule.get_molecular_hamiltonian()

ucc_ccsd = openfermion.utils.uccsd_generator(molecule.ccsd_single_amps,molecule.ccsd_double_amps)
ucc_ccsd_normalordered = openfermion.normal_ordered(ucc_ccsd)

ind1, ind2, ind3, ind4 = np.nonzero(molecule.ccsd_double_amps)
individual_t2_values_ccsd = []
for indice in range(len(ind1)):
 t2_value = abs(molecule.ccsd_double_amps[ind1[indice],ind2[indice],ind3[indice],ind4[indice]])
 individual_t2_values_ccsd.append(t2_value) if abs(t2_value) not in individual_t2_values_ccsd else individual_t2_values_ccsd

with open("ccsd_H4.txt","a") as f:
 f.write("From OpenFermion, t2 amps:\n")
 for indice in range(len(ind1)):
    # double_amplitudes is stored as t[i,j,k,l] * (a_i^\dagger a_j a_k^\dagger a_l - H.C.)
    f.write(str(ind1[indice]) + " " + str(ind2[indice]) + " " + str(ind3[indice]) + " " + str(ind4[indice]) + " " + str(2*molecule.ccsd_double_amps[ind1[indice], ind2[indice], ind3[indice], ind4[indice]]) + "\n")
 f.write("  double uccsd amps:\n" + str(ucc_ccsd) + "\n")
 f.write("  double uccsd amps normal ordered:\n" + str(ucc_ccsd_normalordered) + "\n")

You can check the ccsd_H4.txt file to see the amplitudes and corresponding indices.

[bsenjean]

Eeeeh, I just found the origin of my problem.

If I take the same script as in the previous mail, but setting symmetry c1 in Psi4, then I get the correct ordering (meaning there is no swapping between 4 <--> 6 and 5 <--> 7 anymore). So that's due to the IRREP.

And indeed in OpenFermion-Psi4, in the _psi4_template the keyword symmetry c1 is set.