Charge localization and Jahn-Teller distortions in the benzene dimer cation Piotr A. Pieniazek, Stephen E. Bradforth* and Anna I. Krylov* Department of Chemistry, University of Southern California, Los Angeles, CA 90089-0482 (Dated: July 3, 2008) I. EOM-IP ENERGY EXPRESSIONS In the matrix form, the EOM-IP/EA-CCSD equations assume the following form: ��� ��� ��� �� H SS - ECC �� H SD �� H DS �� H DD - ECC ������ ������ ������ R1(n) R2(n) ��� ��� ��� = ��n ��� ��� ��� R1(n) R2(n) ��� ��� ��� (1) L1(n) L2(n) ��� ��� ��� �� H SS - ECC �� H SD �� H DS �� H DD - ECC ��� ��� ��� = ��n L1(n) L2(n) (2) where the IP/EA superscript is dropped. These equations are usually solved using the David- son iterative diagonalization procedure [1���3], which requires the computation of the ��-vectors, the products of the Hamiltonian and trial vectors. For the EOM-IP-CC(2,3) model the matrix equation for the �� vectors assumes the following form: ��� ��� ��� ��� ��� ��� �� H SS - ECC �� H SD �� H ST �� H DS �� H DD - ECC �� H DT �� H T S �� H T D �� H T T - ECC ������ ������ ������ ������ ������ ������ R1 R2 R3 ��� ��� ��� ��� ��� ��� = ��� ��� ��� ��� ��� ��� ��1 ��2 ��3 ��� ��� ��� ��� ��� ��� (3) (4) The EOM-IP-CCSD model is recovered by setting R3 = 0. The left eigenvalue problem for EOM-IP-CCSD has the form: L1 L2 ��� ��� ��� ��SS H - ECC ��SD H ��DS H ��DD H - ECC ��� ��� ��� = ��1 �� ��2 �� (5) Programmable expressions for the left and right EOM-IP-CCSD ��-vectors, as well as the right EOM-IP-CC(2,3) ��-vectors are given in Sec. III.

2 II. EOM-IP-CC GRADIENTS The Hellmann-Feynman theorem states that once the wave function of a system has been fully variationally optimized, the forces on the system in response to a perturbation �� can be determined by computing the corresponding expectation value using the unperturbed wave function |�� : ���E ����� = ��� �� H �� ����� = �� ���H ����� �� (6) Using the second-quantized Hamiltonian, the energy of the system can be expressed as: E = ��L H ��R = pq hpq��pq + 1 4 pqrs pq||rs ��pqrs (7) ��pq = 1 2 ��L p+q + q+p ��R (8) ��pqrs = 1 2 ��L p+q+sr + s+r+pq ��R (9) where �� pq and �� pqrs are the one and two-particle reduced density matrices that do not depend on the perturbation. The derivative is simply: ���E ����� = pq hpq��pq �� + 1 4 pqrs pq||rs �� ��pqrs (10) where the superscript �� denotes the derivative of the respective integrals. However, this holds true only for wave functions optimized variationally w.r.t all their parameters. The wave functions obtained by electronic structure methods (including EOM-CC) are subject to constraints, and their response to perturbation �� must be accounted for through the additional terms: ���E ����� = �� ���H ����� �� + ����� ����� H �� + �� ���H ����� ����� ����� (11) The direct determination of these terms is inefficient, and the Z-vector [4] and the Lagrangian [5���10] techniques have been developed to compute additional wave function response terms. Below we summarize the Lagrangian-based derivation of the EOM-IP-CCSD gradients follow- ing the presentation in Ref. [10]. The Lagrangian is constructed to incorporate the constraints on the wave function as variational parameters through the undetermined multipliers. In the case of EOM-IP-CCSD, the constraints are: (i) the CC equations are satisfied (ii) the orbitals are eigenfunctions of