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 20import numpy as np 21import sys 22import PyCheMPS2 23import ctypes 24 25# Set the seed of the random number generator and cout.precision 26Initializer = PyCheMPS2.PyInitialize() 27Initializer.Init() 28 29# Read in the FCIDUMP 30psi4group = 5 # c2v: see chemps2/Irreps.h 31filename = b'../../tests/matrixelements/CH4.STO3G.FCIDUMP' 32orbirreps = np.array([-1, -1], dtype=ctypes.c_int) # CheMPS2 reads it in from FCIDUMP 33Ham = PyCheMPS2.PyHamiltonian( -1, psi4group, orbirreps, filename ) 34 35# Define the symmetry sector 36TwoS = 0 # Two times the targeted spin 37Nelec = 10 # The number of electrons 38Irrep = 0 # The targeted irrep 39 40# Setting up the Problem 41Prob = PyCheMPS2.PyProblem(Ham, TwoS, Nelec, Irrep) 42 43# Setting up the ConvergenceScheme 44# setInstruction(instruction, D, Econst, maxSweeps, noisePrefactor) 45OptScheme = PyCheMPS2.PyConvergenceScheme(2) # 2 instructions 46OptScheme.setInstruction(0, 30, 1e-10, 3, 0.1) 47OptScheme.setInstruction(1, 1000, 1e-10, 10, 0.0) 48 49# Do DMRG calculation 50theDMRG = PyCheMPS2.PyDMRG(Prob, OptScheme) 51EnergyDMRG = theDMRG.Solve() 52theDMRG.calc2DMandCorrelations() 53 54# Do FCI calculation 55Nel_up = ( Nelec + TwoS ) / 2 56Nel_down = ( Nelec - TwoS ) / 2 57maxMemWorkMB = 10.0 58FCIverbose = 1 59theFCI = PyCheMPS2.PyFCI(Ham, Nel_up, Nel_down, Irrep, maxMemWorkMB, FCIverbose) 60GSvector = np.zeros([ theFCI.getVecLength() ], dtype=ctypes.c_double) 61GSvector[ theFCI.LowestEnergyDeterminant() ] = 1.0 62EnergyFCI = theFCI.GSDavidson(GSvector) 63theFCI.CalcSpinSquared(GSvector) 64TwoRDM = np.zeros([ Ham.getL()**4 ], dtype=ctypes.c_double) 65theFCI.Fill2RDM(GSvector, TwoRDM) 66RMSerror2DM = 0.0 67for orb1 in range(0, Ham.getL()): 68 for orb2 in range(0, Ham.getL()): 69 for orb3 in range(0, Ham.getL()): 70 for orb4 in range(0, Ham.getL()): 71 temp = TwoRDM[orb1 + Ham.getL()*(orb2 + Ham.getL()*(orb3 + Ham.getL()*orb4))] - theDMRG.get2DMA(orb1,orb2,orb3,orb4) 72 RMSerror2DM += temp*temp 73RMSerror2DM = np.sqrt(RMSerror2DM) 74print("Frobenius norm of the difference of the DMRG and FCI 2-RDMs =", RMSerror2DM) 75 76# Clean-up 77# theDMRG.deleteStoredMPS() 78theDMRG.deleteStoredOperators() 79del theFCI 80del theDMRG 81del OptScheme 82del Prob 83del Ham 84del Initializer 85 86# Check whether the test succeeded 87if ((np.fabs(EnergyDMRG - EnergyFCI) < 1e-8) and (RMSerror2DM < 1e-3)): 88 print("================> Did test 3 succeed : yes") 89else: 90 print("================> Did test 3 succeed : no") 91 92