/*
   CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
   Copyright (C) 2013-2018 Sebastian Wouters

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License along
   with this program; if not, write to the Free Software Foundation, Inc.,
   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/

#include <iostream>
#include <math.h>
#include <string.h>

#include "Initialize.h"
#include "CASSCF.h"
#include "DMRGSCFoptions.h"
#include "MPIchemps2.h"

using namespace std;

int main(void){

   #ifdef CHEMPS2_MPI_COMPILATION
   CheMPS2::MPIchemps2::mpi_init();
   #endif

   CheMPS2::Initialize::Init();

   // Setup the Hamiltonian
   string matrixelements = "${CMAKE_SOURCE_DIR}/tests/matrixelements/N2.CCPVDZ.FCIDUMP";
   const int psi4groupnumber = 7; // d2h -- see Irreps.h and N2.ccpvdz.out
   CheMPS2::Hamiltonian * Ham = new CheMPS2::Hamiltonian( matrixelements, psi4groupnumber );

   // Setup CASSCF --> number of irreps = 8
   int DOCC[]  = { 3, 0, 0, 0, 0, 2, 1, 1 }; // see N2.ccpvdz.out
   int SOCC[]  = { 0, 0, 0, 0, 0, 0, 0, 0 };
   int NOCC[]  = { 1, 0, 0, 0, 0, 1, 0, 0 };
   int NDMRG[] = { 2, 0, 1, 1, 0, 2, 1, 1 };
   int NVIRT[] = { 4, 1, 2, 2, 1, 4, 2, 2 };
   CheMPS2::CASSCF koekoek( Ham, DOCC, SOCC, NOCC, NDMRG, NVIRT );

   // Setup symmetry sector
   int Nelec = 14;
   int TwoS  = 0;
   int Irrep = 0;

   // Run CASSCF
   const int root_num = 1; //Ground state only
   CheMPS2::DMRGSCFoptions * scf_options = new CheMPS2::DMRGSCFoptions();
   scf_options->setDoDIIS( true );
   const double IPEA = 0.0;
   const double IMAG = 0.0;
   const bool PSEUDOCANONICAL = false;
   double Energy1 = koekoek.solve( Nelec, TwoS, Irrep, NULL, root_num, scf_options);
   double Energy2 = koekoek.caspt2(Nelec, TwoS, Irrep, NULL, root_num, scf_options, IPEA, IMAG, PSEUDOCANONICAL);

   // Clean up
   if (scf_options->getStoreUnitary()){ koekoek.deleteStoredUnitary( scf_options->getUnitaryStorageName() ); }
   if (scf_options->getStoreDIIS()){ koekoek.deleteStoredDIIS( scf_options->getDIISStorageName() ); }
   delete scf_options;
   delete Ham;

   // Check succes
   const bool success = (( fabs( Energy1 + 109.103502335253 ) < 1e-8 ) && ( fabs( Energy2 + 0.159997813112638 ) < 1e-8 )) ? true : false;

   #ifdef CHEMPS2_MPI_COMPILATION
   CheMPS2::MPIchemps2::mpi_finalize();
   #endif

   cout << "================> Did test 13 succeed : ";
   if (success){
      cout << "yes" << endl;
      return 0; //Success
   }
   cout << "no" << endl;
   return 7; //Fail

}