1/* 2 CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry 3 Copyright (C) 2013-2018 Sebastian Wouters 4 5 This program is free software; you can redistribute it and/or modify 6 it under the terms of the GNU General Public License as published by 7 the Free Software Foundation; either version 2 of the License, or 8 (at your option) any later version. 9 10 This program is distributed in the hope that it will be useful, 11 but WITHOUT ANY WARRANTY; without even the implied warranty of 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 GNU General Public License for more details. 14 15 You should have received a copy of the GNU General Public License along 16 with this program; if not, write to the Free Software Foundation, Inc., 17 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 18*/ 19 20#include <iostream> 21#include <math.h> 22#include <string.h> 23 24#include "Initialize.h" 25#include "DMRG.h" 26#include "FCI.h" 27#include "MPIchemps2.h" 28 29using namespace std; 30 31int main(void){ 32 33 #ifdef CHEMPS2_MPI_COMPILATION 34 CheMPS2::MPIchemps2::mpi_init(); 35 #endif 36 37 CheMPS2::Initialize::Init(); 38 39 //The path to the matrix elements 40 string matrixelements = "${CMAKE_SOURCE_DIR}/tests/matrixelements/CH4.STO3G.FCIDUMP"; 41 42 //The Hamiltonian 43 const int psi4groupnumber = 5; // c2v -- see Irreps.h and CH4.sto3g.out 44 CheMPS2::Hamiltonian * Ham = new CheMPS2::Hamiltonian( matrixelements, psi4groupnumber ); 45 cout << "The group was found to be " << CheMPS2::Irreps::getGroupName(Ham->getNGroup()) << endl; 46 47 //The targeted state 48 int TwoS = 0; 49 int N = 10; 50 int Irrep = 0; 51 CheMPS2::Problem * Prob = new CheMPS2::Problem(Ham, TwoS, N, Irrep); 52 53 //The convergence scheme 54 CheMPS2::ConvergenceScheme * OptScheme = new CheMPS2::ConvergenceScheme(2); 55 //OptScheme->setInstruction(instruction, DSU(2), Econvergence, maxSweeps, noisePrefactor); 56 OptScheme->setInstruction(0, 30, 1e-10, 3, 0.1); 57 OptScheme->setInstruction(1, 1000, 1e-10, 10, 0.0); 58 59 //Run ground state calculation 60 CheMPS2::DMRG * theDMRG = new CheMPS2::DMRG(Prob, OptScheme); 61 const double EnergyDMRG = theDMRG->Solve(); 62 theDMRG->calc2DMandCorrelations(); 63 64 //Calculate FCI reference energy and compare the DMRG and FCI 2-RDMs 65 double EnergyFCI = 0.0; 66 double RMSerror2DM = 0.0; 67 #ifdef CHEMPS2_MPI_COMPILATION 68 if ( CheMPS2::MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER ) 69 #endif 70 { 71 const int Nel_up = ( N + TwoS ) / 2; 72 const int Nel_down = ( N - TwoS ) / 2; 73 const double maxMemWorkMB = 10.0; 74 const int FCIverbose = 1; 75 CheMPS2::FCI * theFCI = new CheMPS2::FCI(Ham, Nel_up, Nel_down, Irrep, maxMemWorkMB, FCIverbose); 76 double * inoutput = new double[theFCI->getVecLength(0)]; 77 theFCI->ClearVector(theFCI->getVecLength(0), inoutput); 78 inoutput[ theFCI->LowestEnergyDeterminant() ] = 1.0; 79 EnergyFCI = theFCI->GSDavidson(inoutput); 80 theFCI->CalcSpinSquared(inoutput); 81 const int L = Ham->getL(); 82 double * TwoDMspace = new double[ L*L*L*L ]; 83 theFCI->Fill2RDM(inoutput, TwoDMspace); 84 for (int orb1=0; orb1<L; orb1++){ 85 for (int orb2=0; orb2<L; orb2++){ 86 for (int orb3=0; orb3<L; orb3++){ 87 for (int orb4=0; orb4<L; orb4++){ 88 const double difference = TwoDMspace[orb1 + L*(orb2 + L*(orb3 + L*orb4))] 89 - theDMRG->get2DM()->getTwoDMA_HAM(orb1, orb2, orb3, orb4); 90 RMSerror2DM += difference * difference; 91 } 92 } 93 } 94 } 95 delete [] TwoDMspace; 96 delete [] inoutput; 97 delete theFCI; 98 RMSerror2DM = sqrt(RMSerror2DM); 99 cout << "Frobenius norm of the difference of the DMRG and FCI 2-RDMs = " << RMSerror2DM << endl; 100 } 101 #ifdef CHEMPS2_MPI_COMPILATION 102 CheMPS2::MPIchemps2::broadcast_array_double( &EnergyFCI, 1, MPI_CHEMPS2_MASTER ); 103 CheMPS2::MPIchemps2::broadcast_array_double( &RMSerror2DM, 1, MPI_CHEMPS2_MASTER ); 104 #endif 105 106 //Clean up DMRG 107 if (CheMPS2::DMRG_storeMpsOnDisk){ theDMRG->deleteStoredMPS(); } 108 if (CheMPS2::DMRG_storeRenormOptrOnDisk){ theDMRG->deleteStoredOperators(); } 109 delete theDMRG; 110 delete OptScheme; 111 delete Prob; 112 delete Ham; 113 114 //Check success 115 const bool success = (( fabs( EnergyDMRG - EnergyFCI ) < 1e-8 ) && ( RMSerror2DM < 1e-3 )) ? true : false; 116 117 #ifdef CHEMPS2_MPI_COMPILATION 118 CheMPS2::MPIchemps2::mpi_finalize(); 119 #endif 120 121 cout << "================> Did test 3 succeed : "; 122 if (success){ 123 cout << "yes" << endl; 124 return 0; //Success 125 } 126 cout << "no" << endl; 127 return 7; //Fail 128 129} 130 131 132