import psi4
import pathlib
import os

psi4.set_output_file('fci_h2o.dat', False)
psi4.set_memory('120GB')
psi4.set_num_threads(24)

# Reference equilibrium bond length and angle taken from
# Table 1 in JCTC, 2013, 9, 316-323
# Specifically, the CCSDT/jun-cc-pvdz data is used
# Frozen core (1s orbital of oxygen) geometry optimizations
# with convergence parameters (in Psi4 notation)
#  'R_CONVERGENCE': 1e-09
#  'MAX_ENERGY_G_CONVERGENCE': 1.0e-08
#  'MAX_FORCE_G_CONVERGENCE': 1.0e-06
#  'RMS_FORCE_G_CONVERGENCE': 1.0e-06
#  'MAX_DISP_G_CONVERGENCE': 1.0e-06
#  'RMS_DISP_G_CONVERGENCE': 1.0e-06

R_ref = 0.96735
theta_ref = 102.799

mol_string = """
symmetry c2v
0 1
O
H 1 {R}
H 1 {R} 2 {theta}
units angstrom
"""

# Obtain FCIDUMP and FCI energy
parent = pathlib.Path.cwd()
basis = 'CC-PVDZ'

mol = psi4.geometry(mol_string.format(R=R_ref, theta=theta_ref), name='h2o')
molname = mol.name()
# Information string
info = 'Running FCI/{basis:s} for geometry {geometry:s} without freezing the 1s orbital of oxygen.\n'.format(basis=basis, geometry=mol.name())
psi4.core.print_out(info)
psi4.set_options({
  'SCF_TYPE': 'PK',
  'BASIS': basis,
  'E_CONVERGENCE': 1e-09,
  'D_CONVERGENCE': 1e-09,
  'R_CONVERGENCE': 1e-06,
  'REFERENCE': 'RHF',
  'FCI': 'TRUE'
})
scf_e, scf_wfn = psi4.energy('scf', molecule=mol, return_wfn=True)
basename = '.'.join([mol.name(), basis, mol.point_group().symbol()])
fname_fcidump = basename + '.INTDUMP'
psi4.fcidump(scf_wfn, fname=fname_fcidump, oe_ints=['EIGENVALUES'])
e = psi4.energy('detci', ref_wfn=scf_wfn)
psi4.core.clean_options()
psi4.core.clean()