SebWouters
commented
5 years ago Hi @hungpham2017
Be in STO-3G has 4 orbitals, which is a rather small test case, better to just use FCI then and no DMRG.
Once DMRG gets stuck in a local minimum, it's hard to get out. Normally, the test case is big. In that case:
- The MPS state underlying DMRG is large and contains sufficient noise, due to which it's (comparatively) unlikely to get stuck in the wrong state.
- The utilized virtual dimension in DMRG is smaller than the one needed for representing the FCI state via an MPS. In that case, there is a discarded weight in DMRG optimizations near the middle of the MPS. Noise is introduced in the algorithm in an amount 0.5 noise_factor weight_discarded (see e.g. https://arxiv.org/pdf/1405.1225.pdf, section 3.4.3). If there never is a discarded weight, because the virtual dimension is sufficiently large, no noise will be introduced.
In your case Be in STO-3G is so small (4 orbitals), that it is very likely to get stuck in the wrong state and no noise is added as the virtual dimension for representing the FCI solution would be (not taking SU(2) into consideration) 16.
Given the above explanation, please test with a little bit larger system, for example O (oxygen) and 6-31G:
import numpy as np
from functools import reduce
from pyscf import gto, scf, fci, ao2mo
mol = gto.M(atom={'O 0 0 0'}, basis='6-31G')
m = scf.RHF(mol).run()
norb = m.mo_coeff.shape[1]
nelec = mol.nelec
h1e = reduce(np.dot, (m.mo_coeff.T, m.get_hcore(), m.mo_coeff))
g2e = ao2mo.incore.general(m._eri, (m.mo_coeff,)*4, compact=False).reshape(norb,norb,norb,norb)
fs = fci.addons.fix_spin_(fci.FCI(mol, m.mo_coeff), .5)
fs.nroots = 16
e, fcivec = fs.kernel(verbose=0)
import PyCheMPS2
import ctypes, os, sys
Initializer = PyCheMPS2.PyInitialize()
Initializer.Init()
Group = 0
orbirreps = np.zeros([norb], dtype=ctypes.c_int)
HamCheMPS2 = PyCheMPS2.PyHamiltonian(norb, Group, orbirreps)
#Feed the 1e and 2e integral (T and V)
for orb1 in range(norb):
for orb2 in range(norb):
HamCheMPS2.setTmat(orb1, orb2, h1e[orb1, orb2])
for orb3 in range(norb):
for orb4 in range(norb):
HamCheMPS2.setVmat(orb1, orb2, orb3, orb4, g2e[orb1, orb3, orb2, orb4]) #From chemist to physics notation
TwoS = mol.spin
Nel_up = (mol.nelectron+ TwoS) // 2
Nel_down = mol.nelectron - Nel_up
Irrep = 0
maxMemWorkMB = m.max_memory
Prob = PyCheMPS2.PyProblem(HamCheMPS2, TwoS, mol.nelectron, Irrep)
OptScheme = PyCheMPS2.PyConvergenceScheme(4)
OptScheme.setInstruction(0, 200, 1e-4, 5, 0.03)
OptScheme.setInstruction(1, 500, 1e-5, 5, 0.03)
OptScheme.setInstruction(2, 1000, 1e-6, 5, 0.001)
OptScheme.setInstruction(3, 1000, 1e-8, 100, 0.00)
theDMRG = PyCheMPS2.PyDMRG(Prob, OptScheme)
EDMRG0 = theDMRG.Solve()
EDMRG = []
theDMRG.activateExcitations(15)
for i in range(15):
theDMRG.newExcitation(20.0)
EDMRG.append(theDMRG.Solve())
print("PySCF : Root0: %15.8f, Root1: %15.8f, Root2: %15.8f, Root3: %15.8f, Root4: %15.8f, Root5: %15.8f, Root6: %15.8f, Root7: %15.8f, Root8: %15.8f, Root9: %15.8f, Root10: %15.8f, Root11: %15.8f, Root12: %15.8f, Root13: %15.8f, Root14: %15.8f, Root15: %15.8f" % (e[0],e[1],e[2],e[3],e[4],e[5],e[6],e[7],e[8],e[9],e[10],e[11],e[12],e[13],e[14],e[15]))
print("CheMPS2: Root0: %15.8f, Root1: %15.8f, Root2: %15.8f, Root3: %15.8f, Root4: %15.8f, Root5: %15.8f, Root6: %15.8f, Root7: %15.8f, Root8: %15.8f, Root9: %15.8f, Root10: %15.8f, Root11: %15.8f, Root12: %15.8f, Root13: %15.8f, Root14: %15.8f, Root15: %15.8f" % (EDMRG0, EDMRG[0], EDMRG[1], EDMRG[2], EDMRG[3], EDMRG[4], EDMRG[5], EDMRG[6], EDMRG[7], EDMRG[8], EDMRG[9], EDMRG[10], EDMRG[11], EDMRG[12], EDMRG[13], EDMRG[14]))The results are given here:
Run 1
PySCF : Root0: -74.75708817, Root1: -74.75708817, Root2: -74.75708817, Root3: -74.75708817, Root4: -74.75708817, Root5: -74.69661726, Root6: -73.91627244, Root7: -73.91627239, Root8: -73.91627239, Root9: -73.77153553, Root10: -73.77134495, Root11: -73.76997568, Root12: -73.64050889, Root13: -73.63932580, Root14: -73.63835551, Root15: -73.63829142
CheMPS2: Root0: -74.75708817, Root1: -74.75708817, Root2: -74.75708817, Root3: -74.75708817, Root4: -74.75708817, Root5: -74.69661726, Root6: -73.91627281, Root7: -73.91627281, Root8: -73.91627281, Root9: -73.77134991, Root10: -73.77134991, Root11: -73.77134991, Root12: -73.64080297, Root13: -73.64080297, Root14: -73.64080297, Root15: -73.64080297
Run 2
PySCF : Root0: -78.84326124, Root1: -74.75708817, Root2: -74.75708817, Root3: -74.75708817, Root4: -74.75708817, Root5: -74.75708817, Root6: -74.69661726, Root7: -73.91627281, Root8: -73.91627281, Root9: -73.91627278, Root10: -73.77611243, Root11: -73.77133784, Root12: -73.77125551, Root13: -73.64057296, Root14: -73.64043002, Root15: -73.64035129
CheMPS2: Root0: -74.75708817, Root1: -74.75708817, Root2: -74.75708817, Root3: -74.75708817, Root4: -74.75708817, Root5: -74.69661726, Root6: -73.91627281, Root7: -73.91627281, Root8: -73.91627281, Root9: -73.77134991, Root10: -73.77134991, Root11: -73.77134991, Root12: -73.64080297, Root13: -73.64080297, Root14: -73.64080297, Root15: -73.64080297
Run 3
PySCF : Root0: -74.75709041, Root1: -74.75708820, Root2: -74.75708817, Root3: -74.75708817, Root4: -74.75708816, Root5: -74.69661726, Root6: -73.91627299, Root7: -73.91627283, Root8: -73.91627255, Root9: -73.77134943, Root10: -73.77110371, Root11: -73.77094994, Root12: -73.64067206, Root13: -73.64029826, Root14: -73.63934603, Root15: -73.63610232
CheMPS2: Root0: -74.75708817, Root1: -74.75708817, Root2: -74.75708817, Root3: -74.75708817, Root4: -74.75708817, Root5: -74.69661726, Root6: -73.91627281, Root7: -73.91627281, Root8: -73.91627281, Root9: -73.77134991, Root10: -73.77134991, Root11: -73.77134991, Root12: -73.64080297, Root13: -73.64080297, Root14: -73.64080297, Root15: -73.64080297As you can observe, it is now PySCF which gives varying and wrong solutions in different runs, while CheMPS2 is consistent over the consecutive runs.
@sunqm: Any idea?
Best regards, Sebastian
P.S.: The magic number "20.0" in the script above should be about 10 times the difference between the largest energy and the smallest energy you target. Better to make the number too large than too small. In chemistry, it makes sense to set it at fabs(energy_ground_state).
sunqm
hungpham2017
loriab
Hello Seb,
I am trying to compare the PySCF/FCI and CheMPS2/DMRG solver for excited state calculations. When I executed the script over and over, the FCI is stable while the DMRG solver gave random orders of excited state energies, and even sometimes give an excited state energy for the ground state (one root calculation). What would be the possible reason? Maybe there were anything wrong in the way I used it. Thank you, Hung
Here is the script: